SFH of massive quiescent galaxies

Pathways to quiescence: SHARDS view on the star formation histories of massive quiescent galaxies at 1.0 1.5

Helena Domínguez Sánchez , Pablo G. Pérez-González, Pilar Esquej, E-mail: helenads@ucm.es    M. Carmen Eliche-Moral, Guillermo Barro, Antonio Cava, Anton M. Koekemoer,    Belén Alcalde Pampliega, Almudena Alonso Herrero, Gustavo Bruzual, Nicolás Cardiel,    Javier Cenarro, Daniel Ceverino, Stéphane Charlot, Antonio Hernán Caballero
Departamento de Astrofísica, Facultad de CC. Físicas, Universidad Complutense de Madrid, E-28040 Madrid, Spain
University of California, 1156 High Street, Santa Cruz, CA 95064
Observatoire de Genève, Université de Genève, 51 Ch. des Maillettes, 1290 Versoix, Switzerland
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore MD 21218, USA
Instituto de Física de Cantabria, CSIC-UC, 39005 Santander, Spain
Instituto de Radioastronomía y Astrofísica, UNAM, Campus Morelia, México
Centro de Estudios de Física del Cosmos de Aragón, Plaza San Juan 1, 44001, Teruel, Spain
Centro de Astrobiología (CSIC-INTA), Ctra de Torrejón a Ajalvir, km 4, E-28850 Torrejón de Ardoz, Madrid, Spain
Astro-UAM, Universidad Autónoma de Madrid, Unidad Asociada CSIC, E-28049 Madrid, Spain
UPMC-CNRS, UMR7095, Institut d´Astrophysique de Paris, F-75014 Paris, France
Accepted. Received ; in original form

We present Star Formation Histories (SFHs) for a sample of 104 massive (stellar mass M  10 M) quiescent galaxies (MQGs) at  = 1.0–1.5 from the analysis of spectro-photometric data from the SHARDS and HST/WFC3 G102 and G141 surveys of the GOODS-N field, jointly with broad-band observations from ultraviolet (UV) to far-infrared (Far-IR). The sample is constructed on the basis of rest-frame UVJ colours and specific star formation rates (sSFR=SFR/Mass). The Spectral Energy Distributions (SEDs) of each galaxy are compared to models assuming a delayed exponentially declining SFH. A Monte Carlo algorithm characterizes the degeneracies, which we are able to break taking advantage of the SHARDS data resolution, by measuring indices such as MgUV and D4000. The population of MQGs shows a duality in their properties. The sample is dominated (85%) by galaxies with young mass-weighted ages,   2 Gyr, short star formation timescales,   60-200 Myr, and masses log(M/M 10.5. There is an older population (15%) with 2 – 4 Gyr, longer star formation timescales,   400 Myr, and larger masses, log(M/M 10.7. The SFHs of our MQGs are consistent with the slope and the location of the Main Sequence (MS) of star-forming galaxies at   1.0, when our galaxies were 0.5-1.0 Gyr old. According to these SFHs, all the MQGs experienced a Luminous Infrared Galaxy (LIRG) phase that lasts for  500 Myr, and half of them an Ultra Luminous Infrared Galaxy (ULIRG) phase for  100 Myr. We find that the MQG population is almost assembled at   1, and continues evolving passively with few additions to the population.

galaxies: formation; galaxies: evolution; galaxies: high-redshift; galaxies: stellar content
pagerange: Pathways to quiescence: SHARDS view on the star formation histories of massive quiescent galaxies at 1.0 1.5Apubyear: 2015

1 Introduction

In the current paradigm of cosmic evolution, galaxies grow their mass by accreting gas from the cosmic web (e.g. Tacconi et al. 2010) and transforming it into stars. A tight relation between mass and SFR exists for normal star-forming galaxies, known as the MS (e.g. Noeske et al. 2007; Elbaz et al. 2007; Rodighiero et al. 2010). Galaxies grow in mass within or above this MS until, eventually, the exhaustion of gas or some feedback mechanism halt the star formation. At this time, galaxies reach a quiescence state. The mass evolution of the quiescent population is then dominated by (dry) mergers of already assembled (smaller or similar in mass) galaxies. Lambda Cold Dark Matter (CDM) theory requires more massive halos to generally assemble later than less massive ones. However, over the past few decades there have been a series of observational results indicating that many processes, such as star formation, occur earlier in the most massive and luminous galaxies than in less massive galaxies, in a scenario called downsizing (e.g., Cowie et al. 1996; Madau et al. 1996; Cimatti et al. 2006; Bundy et al. 2007; Pérez-González et al. 2008; Ilbert et al. 2010).

Focusing on the evolution of the most massive galaxies along the lifetime of the Universe, we now know that many of such systems were already assembled at high redshifts (Pérez-González et al., 2008; Marchesini et al., 2009; Ilbert et al., 2013; Muzzin et al., 2013; Tomczak et al., 2014; Grazian et al., 2015). Moreover, a significant fraction of high-z massive galaxies were not actively forming stars and were evolving passively (Cimatti et al. 2004; Daddi et al. 2005; Papovich et al. 2006; Fontana et al. 2009; Santini et al. 2009; Domínguez Sánchez et al. 2011; Ilbert et al. 2013), i.e., they were MQGs. The number densities of MQGs at intermediate redshifts are in disagreement with semi-analytical models (e.g. Pozzetti et al. 2010; Ilbert et al. 2013; Muzzin et al. 2013). Besides, the observed MQGs at high- are found to be much more compact than their local analogues, implying a strong mass-size evolution with cosmic time (e.g. Trujillo et al. 2006; Buitrago et al. 2008; van Dokkum et al. 2008; Barro et al. 2014).

Understanding the formation mechanisms of this population of MQGs is fundamental to improve our picture of galaxy evolution. In particular, having an accurate description of the SFHs of these galaxies is crucial to have a good estimation of the time needed by physical processes to ignite the star formation and then quench it in massive galaxies. However, up to date, there are few results on the individual properties and SFHs of MQGs at high-, mostly because estimating SFHs is very hard with the data we typically have at hand.

Quiescent galaxies at high- are difficult to observe in the optical bands as their emission in the rest-frame UV is weak (due to the absence of star formation processes), and dominated by absorption features (e.g., Daddi et al. 2005; Cimatti et al. 2008; van Dokkum et al. 2011). Indeed, these galaxies show very weak or no emission lines (Kriek et al., 2009; Belli et al., 2014) and absorption features related to an old and passively evolved stellar population. Spectroscopic observations are time consuming and are therefore limited to a small number of galaxies or to the analysis of stacked spectra (Kriek et al., 2009; Onodera et al., 2012; Toft et al., 2012; Bedregal et al., 2013; Whitaker et al., 2013; Onodera et al., 2015; Belli et al., 2015).

Thanks to the arrival of deep and wide multi-wavelength photometric surveys such as COSMOS (Scoville et al., 2007), ULTRA-VISTA (McCracken et al., 2012) or CANDELS (Grogin et al., 2011; Koekemoer et al., 2011), the mass functions and number densities of quiescent galaxies have been studied up to   3-4. But the spectral resolution of photometric data is not enough to study in detail the stellar population properties. The main reason is the presence of strong degeneracies in the analysis of their SEDs using stellar population synthesis models. These degeneracies are mainly related to the similar effect of different levels of dust attenuation, age, and metallicity in very low spectral resolution data (R  7, typical of broad-band studies). This indeed complicates the accurate determination of high- galaxy properties based on SED-fitting; being the stellar mass the only reliable parameter obtained in these kind of studies (Elsner et al., 2008; Santini et al., 2015).

In order to break the typical degeneracies inherent to any study of stellar populations in distant galaxies, we need data with higher spectral resolution than broad-band photometry (Pacifici et al., 2013). There are spectral features which help to break these degeneracies. For example, the  index probes several absorption lines (e.g., MgI 2852, MgII 2796, 2804, FeII 2600, 2606) and has been shown to be extremely reliable to identify galaxies dominated by evolved stars. Moreover, these absorption lines can be used to easily distinguish the SED of a MQG from the featureless spectrum of a dusty starburst (Daddi et al., 2005). The break in the stellar continuum at 4000 , D4000, is also a good age indicator (Bruzual 1983; Balogh et al. 1999; Kauffmann et al. 2003). It arises because of the accumulation of a large number of spectral lines in a narrow wavelength region. Given the dependence of D4000 on stellar atmospheric parameters (Gorgas et al., 1999), it is very prominent for galaxies older than  1 Gyr. The strength of the 4000 Å  break for a single stellar population increases with its age and depends weakly on the metallicity at low ages ( 1 Gyr, Hernán-Caballero et al. 2013). The  and D4000 indices have been successfully used in the past to obtain redshifts and ages of stellar populations in massive galaxies at high- (Saracco et al., 2005; Kriek et al., 2011; Ferreras et al., 2012).

In this work, we take advantage of the spectro-photometric resolution of the Survey for High- Absorption Red and Dead Sources survey (SHARDS, Pérez-González et al. 2013) to obtain robust estimations of the SFH of MQGs at high-. SHARDS is an ultra-deep optical survey of the GOODS-N field covering the wavelength range between 500 and 950 nm with 25 contiguous medium-band filters, providing a spectral resolution R  50. We combine the SHARDS observations with spectroscopic data from the HST WFC3 G102 and G141 grisms covering 900-1600 nm, as well as with multi-wavelength ancillary data from UV to Far-IR extracted from the Rainbow database (Pérez-González et al., 2008; Barro et al., 2011). We, therefore, use data with a spectral resolution R  50 or better from 500 nm to 1600 nm (jointly with broad-band photometry), a significant improvement from previous works on the subject. We perform a SED-fitting to the whole wavelength range (up to the IRAC bands) using delayed exponentially declining stellar population models. Thanks to the unique photometric dataset and spectral coverage of this work, we are able to account for the degeneracies and study in detail and individually the SFHs (masses, SFRs, ages, star formation timescales, dust attenuations and metallicities) of a sample of MQGs at 1.0 – 1.5, discussing the implications for the early mass assembly of galaxies.

The structure of this paper is as follows: in Sect. 2, we present the data and the sample selection. In Sect. 3, we explain the SED-fitting method and the procedure to characterize the degeneracies inherent to the stellar population analysis. In Sect. 4, we analyse the derived properties and the time-evolution of our sample of MQGs on the basis of their SFHs. Finally, in Sect. 5, we summarize our conclusions.

Throughout this work, we standardize to a (, , )(0.7, 0.3, 0.7) Wilkinson Microwave Anisotropy Probe (WMAP) concordance cosmology (Spergel et al., 2003), AB magnitudes, (Oke & Gunn, 1983), a Kroupa 2001 IMF (integrated from 0.1 – 100 M) and Bruzual & Charlot (2003) (BC03, hereafter) stellar population synthesis models.

2 Data and Sample selection

In this work, we analyse a sample of MQGs (M  10 M) at selected with two different criteria: their rest-frame colours and their sSFR. In the following subsections, we describe the datasets gathered for this work, as well as the details about the selection criteria.

2.1 Data

Our sample was selected in the 130 arcmin covered by SHARDS (Pérez-González et al., 2013) in the GOODS-N region. SHARDS is an ESO/GTC Large Program carried out with the OSIRIS instrument on the 10.4 m Gran Telescopio Canarias (GTC). It consists of an ultra-deep optical spectro-photometric survey of the GOODS-N field at wavelengths between 500 and 950 nm using 25 contiguous medium-band filters which provide a spectral resolution R  50. The data reach an AB magnitude of 26.5 (at least at a 3 level) with sub-arcsec seeing in each one of the 25 bands. More details about the reduction and calibration of the SHARDS data can be found in Pérez-González et al. (2013).

For this paper, we complement the SHARDS data with the G102 and G141 grism observations of the GOODS-N field carried out with the HST/WFC3 instrument. The G102 spectroscopic program (PI: Barro) covers the spectral region between 800 and 1150 nm with R  210 and a 5 magnitude limit of 21.5 mag in the F140W band. The G141 data (PI: Weiner) covers from 1100 to 1700 nm with R  130 and a 5 limit of 21.5 mag in the F160W band. The WFC3 grism spectroscopic data were reduced using the aXe software version 2.3. Based on the 2D spectra provided by aXe, we extracted 1D spectra for each galaxy in our sample using the effective radius as the extraction width. We visually inspected each spectrum, adjusting the extraction parameters (aperture width and spectral range), to avoid contamination from nearby sources. We refer the reader to Esquej et al. (in preparation) for a detailed description of the reduction and extraction of WFC3 grism data.

Jointly with the SHARDS and WFC3 grism data, we also use in our analysis the ancillary multi-wavelength catalogue and advanced products in the GOODS-N field presented in Pérez-González et al. (2008) and compiled in the Rainbow Cosmological Surveys Database111http://rainbowx.fis.ucm.es (see also Pérez-González et al., 2008; Barro et al., 2011). This dataset includes observations from X-rays to the Far-IR and radio bands, as well as spectroscopic data in the GOODS-N field from the literature. In Pérez-González et al. (2008), merged photometric catalogues were presented, including broad-band data for a stellar mass selected sample based on ultra-deep IRAC observations. For this work, we have merged this catalogue with the SHARDS and WFC3 grism datasets. Using the full SED, we carried out a stellar population modeling which provided accurate photometric redshifts, stellar masses, SFRs, and rest-frame synthetic colours for  26,000 stellar mass selected galaxies in the region of the GOODS-North field covered by SHARDS.

Concerning photometric redshifts, the ultra-deep medium band SHARDS data allowed us to obtain high quality photometric redshifts for all sources, which will be presented in Barro et al. (in preparation). The median /(1+z) is 0.0067 for the 2650 sources with I  25 and spectroscopic redshifts (see also Ferreras et al., 2014).

As described in Pérez-González et al. (2008), the parent sample used in this paper, selected with ultra-deep IRAC imaging as a proxy for stellar mass, is complete for galaxies with M  10 M up to and a maximally old stellar population. The IRAC selection is biased against younger galaxies with large attenuations and masses M  10 M, so we impose this mass limit in the definition of our final sample of MQGs.

Synthetic rest-frame colours were estimated for all galaxies in the parent sample by convolving the best-fitting stellar population models with transmission curves for standard filters. In this paper, we will use the Johnson and filters, as well as the 2MASS -band filter to construct a diagram and select quiescent galaxies. We note that the actual transmission curves for which we estimated rest-frame absolute magnitudes were taken from EAZY (Brammer et al., 2008) and FAST (Kriek et al., 2009) 222The central wavelengths (widths) in nm of these filters are 359.84 (58.36), 549.02 (85.79), 1237.59 (169.48) for the bands, respectively., in order to match the colour distribution in the diagrams used in Whitaker et al. (2011), from which we extracted the quiescence definition.

Figure 1: SFR versus redshift (left panel) and IRX- relation (right panel) for the sample of massive (M  10 M) galaxies at detected by MIPS at 24 m. The blue line is the corresponding SFR for the 24 m MIPS detection limit (30 Jy at 5). The red points (23 galaxies) are IR-faint emitters, i.e., galaxies with an IR emission lower than the 5 limit, selected to derive the IRX- relation. To avoid biasing the fits in the IRX- plot, we randomly chose a subsample of faint IR emitters homogeneously covering the whole range of UV slopes . We fit the IRX- (log linear) relation for the faint IR emitters (red line) to obtain a -based dust attenuation for the population of galaxies with low levels of star formation. The green line is the IRX- fit from Meurer et al. (1999), based on the analysis of regular star-forming galaxies and starbursts. Note that the dust attenuation correction at high values is much smaller for the IR-faint sample compared to IR-bright galaxies. Typical uncertainties are shown as error bars in the right corner of each panel.

Finally, in order to select quiescent galaxies for this paper, SFRs were estimated for all galaxies in the parent sample based on either mid- and far-IR data from Spitzer and Herschel, or from UV luminosities. In the case of IR emitters, GOODS-N has been observed with the deepest MIPS data in the sky, with a 5 limit of 30 Jy (Pérez-González et al., 2005), and also very deep PACS and SPIRE images (Elbaz et al., 2011), with the following 5 limits: 1.7 and 3.6 mJy for PACS 100 and 160 m, and 9, 12, and 11 mJy for SPIRE 250, 350, and 500 m (these limits include the effect of confusion in the SPIRE bands).

We looked for counterparts for the galaxies in our mass selected sample in the Spitzer MIPS bands using a search radius of 2″  for 24 m, as described in Pérez-González et al. (2005, 2008). For Herschel PACS and SPIRE bands, we built merged photometric catalogues using position priors from the MIPS data, as described in Pérez-González et al. (2010) and Rawle et al. (2015). The Herschel catalogues are then linked to the MIPS catalog, so we used the same search radius. We only considered detections above the 5 flux limit for each band.

Infrared Luminosities, L, were estimated from the Spitzer and Herschel data by fitting dust emission models from Chary & Elbaz (2001) to all available photometric data points with rest-frame wavelength longer than 6 m. We also checked that similar results (with typical differences smaller than 0.2 dex) were obtained using the templates from Rieke et al. (2009) and Dale & Helou (2002). When we only had a single photometric data point (in almost all cases, 24 m) we scaled the models from Chary & Elbaz (2001) to the monochromatic luminosity probed by that observation. We compared our IR-based SFRs with those obtained by applying the method described in Wuyts et al. (2011); Rujopakarn et al. (2012) and Elbaz et al. (2011), obtaining a good ( 0.2 dex) agreement.

The SFR for IR detected galaxies was derived following Kennicutt (1998, see also ) normalized to a Kroupa (2001) IMF:


where is the luminosity at 280 nm rest-frame in erg s.

In order to estimate SFRs for galaxies not detected in the mid- or far-IR, we used UV luminosities at 280 nm rest-frame alone. These luminosities were converted into SFRs by applying Kennicutt (1998) equation normalized to a Kroupa (2001) IMF:


The UV-based SFRs were corrected for dust attenuation following a recipe based on the relation between the UV slope, , and the UV/IR ratio, IRX, know as the IRX- relation. The UV slope for each galaxy is calculated using a linear interpolation between 150 and 280 nm in the templates fitting the SED of each galaxy (from which a photometric redshift and stellar mass estimate were obtained). The typical uncertainty in the values is  20%. In order to convert these slopes to a dust attenuation, we used the comparison between the observed UV and IR-based SFRs for galaxies detected by MIPS (Figure 1). Typical attenuation recipes based on the UV slope are derived for galaxies with high levels of star formation (Meurer et al., 1999; Takeuchi et al., 2012), which are dustier than relatively quiescent galaxies. Quiescent galaxies have systematically lower ratio of total far-IR to UV luminosity than starburst galaxies (Kong et al., 2004). Given that in this paper we are interested in galaxies with very low levels of star formation, and in order to avoid an over-correction for dust attenuation, we derived a IRX- relation for faint-IR emitters, i.e., those galaxies for which their IR detection is below the 5 IR detection limit at that redshift. In comparison with bright-IR galaxies or the local starbursts, this sub-sample should present dust attenuation properties which are closer to the MQGs that we want to study. The procedure is outlined in Figure 1, where we show SFR versus redshift, and the IRX- relation for the galaxies at detected by MIPS. We highlight the galaxies which are faint-IR emitters, for which we derive the following IRX- relation:


We apply the Meurer et al. (1999) IRX- relation for values lower than the cross-point (=-0.97) and Eq. 3 for higher values.

2.2 Selection of quiescent galaxies

To construct a complete and uncontaminated sample of MQGs at z1.0 – 1.5 we used two complementary methods: a diagram and sSFRs. We only consider galaxies with M  10 M, for which our survey is complete within the considered redshift range.

Figure 2: Plots showing the criteria to select the MQG sample. Left panel: diagram for galaxies at z1.0 – 1.5 (grey dots). The final sample of MQGs (104 objects) are marked with black empty circles. The continuous line delimits the region where quiescent galaxies should be located, according to Whitaker et al. (2011). The dashed line splits the previous region in two, with post-starburst galaxies located on the left. The -selected galaxies are coloured in green, and the sSFR-selected in red. The large blue circles mark IR detected galaxies (only for galaxies with sSFR  0.2 Gyr). Right panel: sSFR versus stellar mass, colour coded as in the left panel. The dashed line represents our limit to consider a galaxy as quiescent (sSFR  0.2 Gyr). Uncertainties in the colours calculated by propagating the redshift uncertainties are negligible, given the superb quality of our SHARDS-based photometric redshifts. Therefore, we assumed for our rest-frame absolute magnitudes the average uncertainty of the two observed filters which lie closer to the central wavelength probed by the rest-frame filters (shown as error bars in the right corner of the left panel).
Figure 3: 5″5″postage stamps made with HST ACS/WFC3 for the 104 MQGs selected in this work (there are 2 sSFR-selected galaxies missing because they are outside the WFC3 region). Galaxies in the upper panel are -selected, while galaxies in the lower panel are sSFR-selected. In each panel, they are ranked (from left to right, top to bottom) by increasing mass-weighted age  (see Sect. 4). Postage stamps for the 3 sub-samples of galaxies defined in the text based on  are framed in blue (mature galaxies), green (intermediate), and red (senior galaxies).

2.2.1 Quiescence criterion based on the colour-colour diagram

There are 410 galaxies with M  10 M at z1.0 – 1.5 in the Rainbow catalogue of the GOODS-N field. The median stellar mass and sSFR of this sample are M  10 M and sSFR  Gyr. In Figure 2, we plot a diagram including all these galaxies. Whitaker et al. (2011) defines the following as the region in the diagram where quiescent galaxies are located:


We find 87 galaxies with M  10 M in the quiescent region of the diagram. Out of these, 25% are detected by MIPS, implying either some level of (obscured) star formation or nuclear activity. We eliminate them from this -selected sample of quiescent galaxies. The sample of MQGs selected with the diagram and no IR detection is composed of 65 galaxies. The median and quartile stellar mass for this sample is log(M/M, they lie at , and have sSFR Gyr. We will refer to the galaxies selected in this way as “UVJ-selected” galaxies and we will discuss their properties in Sect. 4.2.

2.2.2 Quiescence criterion based on sSFR values

We complement the selection based on the diagram with a cut in sSFR. We arbitrarily choose sSFR  0.2 Gyr as our limit for quiescence, given that by imposing this value we are able to virtually recover all the sources identified as dead by the criterion (see right panel of Figure 2).

We find 102 galaxies with M  10 M and sSFR  0.2 Gyr. All galaxies identified as quiescent in the diagram have sSFR  0.2 Gyr, except 2. These 2 galaxies have relatively young stellar populations (t  0.8 Gyr) and short star formation timescales (  60 Myr), as we discuss in Sect. 3. Note that among the sSFR-based sample we have 22 IR emitters (21% of the sSFR-selected sample). These galaxies must have some (residual) star formation or nuclear activity, but their mass is large enough (log(M/M 10.4) to present sSFRs values comparable to dead galaxies undetected in the IR. The median and quartiles for the stellar mass distribution of this sample is log(M/M M, they lie at , and have sSFR Gyr. With the sSFR criterion, we recover 8 galaxies located in the quiescent region which were discarded in the colour-colour selection because of their IR detection. However, they have low enough sSFR values ( 0.18 Gyr) to be selected in our MQG sample. We will refer to galaxies selected in this way as “sSFR-selected” galaxies, and we will discuss their properties in Section 4.2. Note that this is a complementary sample to the -selected, i.e., they are galaxies with sSFR  0.2 Gyr but excluding the 65 previously selected with the colour criteria, which results in 39 additional galaxies.

Combining the and sSFR criteria, we arrive to a final sample of MQGs with M  10 M composed of 104 sources. In Figure 3 we show the postage stamps made with HST ACS/WFC3 for this final sample of MQGs.

2.2.3 Statistical properties of the sample

The combined - and sSFR-selected sample of MQGs was refined by carrying out a detailed stellar population synthesis analysis of each galaxy. We constructed the most detailed SED possible for each galaxy combining the ultra-deep SHARDS data with the G102 and G141 spectro-photometric observations. To increase the S/N of the grism spectra, we binned them in order to have 10 nm per pixel. This corresponds to one and two resolution element for G102 and G141, respectively. This provides SEDs with up to 150 photometric points at a photo-spectral resolution from 500 nm to 1700 nm. Galaxies without G141 or G102 spectra have at least 30 photometric data points. This unique photometric dataset encompasses, within the whole observed redshift range, significant spectral features related to the age of the galaxies, such as D4000 or . Our final sample of MQGs at   consists of 104 galaxies, 65 -selected plus 39 sSFR-selected. There are 54 (52%) galaxies for which spectroscopic redshifts are available, and the photometric redshift quality for them is characterized by a median 0.0047. They are detected at 3 in at least 13 SHARDS bands. Concerning the availability of grism data, 60 and 70% of the final sample have usable G102 and G141 spectra with at least S/N  3 and median S/N  10 per pixel (i.e., 10 nm). The rest have either severe contamination problems or are too faint for the grism observations. MQGs represent  25 of the population of massive galaxies within the same redshift interval. They have median values log(M/M) M, , and sSFR Gyr.

The fraction of quiescent galaxies in a mass-selected sample according to our analysis is in good agreement with the 28 reported by Muzzin et al. (2013) for a purely selected sample of quiescent galaxies with log(M/M 9.5 at  . The fraction of MQGs that we find is larger than the 13 of quiescent galaxies with log(M/M 9.6 at   selected on the basis of their NUV, r and J colours found by Ilbert et al. (2013). The fraction of very massive (log(M/M 10.85 M) quiescent galaxies (selected on the basis of IR colours and sSFR   0.01 Gyr) derived in Fontana et al. (2009) at is %, larger than the  25% that we obtain with the same mass cut. We will further discuss our results about the properties of our sample of MQGs and compare them with the literature in Sect. 4.3.

Figure 4: Example of SED-fitting for one of our MQGs, SHARDSJ123727.9+622034.7, a galaxy at 1.148. In the left panel we show the SED-fitting to the whole dataset gathered for our work (194 data points). The red line is the best-fitting model. The white circles represent the SHARDS spectro-photometric data and the diamonds are broad-band data. The G102 and G141 spectra are plotted as black and dark grey lines and zoomed in the central and bottom right panels. We depict 3 errors in all plots (and the average S/N in the grism spectrum plots). The vertical dotted lines show the location of typical emission lines (H, H and [OII]3727). We also show a 5″5″postage stamp made with HST ACS/WFC3 data. In the upper right panel, we show a zoom in the SHARDS region. The coloured areas represent the bands used to determine the  and D4000 indices, whose values are given in the legend. For this galaxy, only one stellar population model was compatible with the data after applying the method described in Sect. 3.2. The best-fitting parameters are shown in the legend, including metallicity, star formation timescale, age, mass-weighted age, dust attenuation, stellar mass, predicted 24 m flux (see Sect. 3.3) and statistical significance of the solution.

3 Methodology to determine the Star Formation Histories of Massive Quiescent Galaxies

3.1 SED-fitting

Our main goal in this paper is to study in detail the SFHs of MQGs at . For this purpose, we fit the observed photometry to stellar population synthesis models using the synthesizer fitting code described in Pérez-González et al. (2003, 2008). We use models with one burst of star formation characterized by a delayed exponentially declining SFH:


The most common SFH parameterization used in the literature is an exponentially declining function. However, we have chosen a more realistic parameterization with an initial increase in the star formation activity followed by an exponential decay, avoiding the nonphysical infinite derivative at time equals zero obtained in the pure exponential. Our parameterization is also closer to the SFHs predicted by galaxy evolution models (e.g. Pacifici et al. 2013; Barro et al. 2014).

We compare our SEDs with the BC03 stellar population library, assuming a Kroupa (2001) IMF and the Calzetti et al. (2000) attenuation law. We allow 3 different metallicity values, Z/Z = 0.4, 1.0 and 2.5, which correspond to values of sub-solar, solar and super-solar metallicity. Although it is commonly assumed in the literature that galaxies have solar metallicity, it is well known that galaxies may have different metal contents. For example, the mass-metallicity relation (Tremonti et al., 2004) points out that, at low redshift, more massive galaxies are more metal-rich than less massive ones. The metallicity also affects the shape of the SED, making galaxies look redder as we move to higher metallicities. Although the metallicities may take on more values than assumed, using only the 3 discrete values around solar metallicity given by the BC03 models is a sufficient approximation to see the effect of this parameter in our results.

Our stellar population synthesis code, synthesizer, performs a minimization and returns the model which best fits the data. The age (time since the formation of the stellar population that dominates the SED of the galaxy, t), star formation timescale (), dust attenuation (A) and metallicity (Z) are set as free parameters. The parameter space allowed in the fitting procedure is given in Table 1. The stellar mass (M) of each galaxy is derived as the normalization of the best-fitting template to the median observed flux.

In the following sections, apart from the stellar population parameters mentioned above, we will also discuss mass-weighted ages (). The  is defined as:


where t is the best-fitting age, and M is the best-fitting mass. While t corresponds to the beginning of the star-formation period that dominates the SFH of the galaxy,  is a better approximation to the average age of the stellar population and is a more useful parameter to understand the evolution of galaxies taking into account the duration of the star formation processes.

An example of the SED-fitting is presented in Figure 4, where we show the fit to the whole wavelength interval, as well as zooms in the spectral regions covered by SHARDS, and the G102 and G141 grisms. Our fitting code also includes an algorithm to estimate uncertainties in the derived parameters and to study (and break when possible) the degeneracies typically present in SED-fitting studies. This algorithm involves a Monte Carlo technique, and the usage of direct measurements of spectral indices (especially D4000 and ) with the SHARDS and WFC3/grism spectro-photometric data, as well as the analysis of the mid- and far-IR fluxes and upper limits. We describe this algorithm in detail in the following subsection.

Parameter range/values units step
age () 0.04 – 6.3 Gyr 0.1 dex
timescale () 3 – 10000 Myr 0.1 dex
dust attenuation (A) 0 – 1.5 mag 0.1 mag
metallicity (Z) 0.4, 1.0, 2.5 Z/Z discrete
Table 1: Free parameters (age, star formation timescale, dust attenuation and metallicity) and their allowed ranges used in the SED-fitting procedure.

3.2 Estimating uncertainties and analysing degeneracies with a Monte Carlo algorithm

We used a Monte Carlo approach to estimate uncertainties in the stellar population properties and to take into account the possible degeneracies of the SED-fitting technique. For each galaxy, we constructed 1000 modified SEDs by allowing the photometric data points to randomly vary following a Gaussian distribution, with a width given by the photometric errors. We performed the SED-fitting to the modified photometric data and obtained 1000 different solutions with their corresponding set of parameters for every galaxy in our sample. In the SED-fitting procedure, the stellar population modeling code looks for best-fitting ages, timescales, metallicities, and dust attenuations in a grid of discrete values.

Once we have 1000 SED-fitting solutions for a given galaxy, we look for clusters of solutions in the -t parameter space. In order to account for the discrete distribution of the fitting parameters probed by the minimization algorithm, we introduce Gaussian noise to the output parameters of each galaxy using a width equal to the step used for each fitted property (see Table 1). We identify clusters in the -t plane with a k-means method and a minimal separation of 0.2 dex between different solutions (i.e., the difference between the median cluster properties must be at least 0.2 dex in age and ). Solutions which provide similar results are grouped as a single solution identified by a median value and a scatter in the multi-dimensional t---Z space. Using the full set of solutions for a given cluster, we calculate the values enclosing 68% of the data around the median. These are assumed to be the uncertainties of our estimations for each cluster of solutions. The typical relative uncertainties of the parameters in our analysis are t=12%, =16%, A=0.06 mag and M=0.05 dex. The metallicity uncertainties cannot be determined because the allowed Z values are discrete.

Each cluster is assigned a statistical significance given by the fraction of solutions belonging to that cluster. The clusters represent the significant solutions of a galaxy taking into account the degeneracies in the determination of the stellar properties from the SED-fitting. Although we use the t- plane to look for different solutions, we tested whether looking for clusters in the whole t---Z multi-dimensional space changed the results. Given that all the parameters are highly correlated, we found no difference between both approaches, i.e., a cluster analysis in one plane was able to robustly recover clusters in the 4-dimensional space.

Figure 5: Example of our SED-fitting method, including the analysis of degeneracies, for SHARDSJ123620.3+6200844.3, a galaxy at . Colours and symbols as in Figure 4. For this galaxy, three clusters of solutions were obtained after applying the method described in Sect. 3.2. The main parameters of each solution are shown in the legend.

In Figure 5, we show the SED-fitting result for a galaxy for which we find three clusters of solutions. The most significant (74% of solutions belong to this cluster) is a relatively young stellar population (  1.4 Gyr) with star formation timescale   160 Myr, moderate dust attenuation ( 0.3 mag), and super-solar metallicity. Another solution consistent with the data is an older population (  2.6 Gyr) with a longer star formation timescale (  310 Myr), higher dust attenuation ( 0.6 mag), and sub-solar metallicity. And finally, another possible solution is characterized by an intermediate age population (  1.9 Gyr) with very short star formation timescale (  24 Myr), higher dust attenuation ( 0.7 mag), and also sub-solar metallicity. We remark that all 3 solutions fit the data equally well, although we must warn the reader that most of the spectro-photometric data points are bluer than 1 m rest-frame, so the fits are biased towards the bluer part of the SED. In Figure 6, we show the procedure used to identify the three clusters of solutions. We plot in the t- space the solutions obtained for the 1000 Monte Carlo simulations. The 3 clusters mentioned above can be identified by colours, and we also depict the median values of each solution.

Figure 6: Example of the procedure used to break degeneracies among best-fitting solutions using measurements of absorption indices for SHARDSJ123620.3+6200844.3. Upper panel: we show a zoom of the SED, centered in the SHARDS region, for the same galaxy depicted in Figure 5. The shaded bands mark the regions used to measure the  and D4000 indices. The coloured circles represent the data points used to measure the indices from the SHARDS photometry. Lower left panel: Distribution of the 1000 best-fitting solutions in the age/timescale plane. Note that this is the direct age from the SED-fitting (t) and not . The coloured dots correspond to each one of the identified clusters, with the same colour as in the SED plot. For each one of the three clusters, the median values are represented by large green circles and the black contours enclose 68% of the solutions. Lower right panel: we show the evolutionary tracks in the -D4000 plane for the three best-fitting models found for SHARDSJ123620.3+6200844.3. The empty circles show the expected values of the indices at different ages, t (from smaller to larger symbols: 0.5, 1.0, 2.0, 3.0 and 5.0 Gyr). The large coloured diamonds show the location of the indices at the age of the best-fitting solution. The black circle represents the indices measured from the SHARDS photometry. In this case, unlike in the following plots where we will use , t is the appropriate parameter to separate clusters of solutions and to identify the evolutionary tracks consistent with the indices, as t is the actual best-fitting parameter of the stellar population models. The indices predicted by the track of the first solution are incompatible with the values measured directly from photometry and therefore, this solution is discarded.

We investigate how our results are affected by the degeneracies by analysing the clusters of solutions. We obtain 176 different solutions for the 104 galaxies in our sample. This means that each galaxy has, on average, 1.7 solutions. There are 48 galaxies ( 46 of the sample) for which we find only one cluster of solutions, i.e., only one set of properties fits the data, after taking into account observational uncertainties and the degeneracies linked to them. There are 41 galaxies (39) with 2 clusters of solutions and less than 15 of the sample has 3 or more clusters (with the maximum number of clusters being 4, which happens only for two galaxies). The fact that almost half of the sample presents no degenerate solutions reveals the power of combining SHARDS and WFC3 data to constrain the properties of MQGs at high-z. When using broad-band data alone, the number of galaxies with degenerate solutions increases up to 76%.

We derived the maximal difference of the galaxy parameters obtained from the different clusters of solutions for each galaxy. The mean relative differences are t  40%,   30%,   80%, A=0.1 mag, M   0.1 dex. The differences between solutions for the stellar masses and the dust attenuations are of the order of the typical uncertainties of each cluster of solutions, meaning that these parameters are not strongly affected by the degeneracies. In fact, only two galaxies with degenerate solutions present A values larger than 0.5 mag. The degeneracies in t are smaller than 50% in relative values (and even smaller (30%) for the ). The largest differences are found for the star formation timescale, . Again, the metallicity uncertainties cannot be determined because the allowed Z values are discrete. The metallicity values are unique (i.e., Z =0) for 55 of the galaxies with more than one cluster of solutions.

3.3 Using spectral indices and IR detections to break degeneracies

The SED-fitting technique described in the previous section does not make use of the full power of the dataset gathered for this work. Including all photometric data points from the UV to the mid-IR gives us the global shape of the stellar emission. But this global shape may easily wash out the higher spectral information given by the spectro-photometric data from SHARDS and the grism observations. Thus, the stellar population properties can be even better constrained when taking advantage of the ultra-deep spectro-photometric SHARDS and grism data. The high resolution (R  50) photometry from SHARDS and the grisms allows us to measure spectral indices such as  and D4000, which are well correlated to stellar population properties (Kauffmann et al., 2003; Daddi et al., 2005).

We measure the  index presented in Daddi et al. (2005) as a strong age-dependent feature in the UV, detectable in relatively low resolution spectra in the region 260-290 nm. This feature is an absorption band formed by a combination of several strong Mg and Fe lines. The  index is defined as:


where the integration ranges are defined in nm. Note that to measure the  index we need 10 nm windows at rest-frame, which at redshift  1 translates to 20 nm or more (observed-frame). This means that the SHARDS filters, with a typical width around 16 nm, present a sufficient spectral resolution to measure the   absorption index.

We also measured the D4000 index, introduced by Bruzual (1983), which gives an estimate of the strength of the 4000  break. This index is defined as:


where the integration ranges are again defined in nm.

We prefer to use this definition instead of the narrow index D4000 (Balogh et al., 1999) to reduce the uncertainty measurement of the index using photometric data alone, thanks to the broader index bands. See Hernán-Caballero et al. (2013) for more details about measuring D4000 and D4000 with SHARDS data.

To measure the spectral indices based on our spectro-photometric data, we first select the data points for each galaxy that fall into the bands used in the index definition. These bands depend on the redshift of the galaxy, but we note that for SHARDS we also have to consider the position of the galaxy in the FOV, given that the central wavelength of the pass-band seen for each galaxy depends on the position in the OSIRIS focal plane (see Pérez-González et al. 2013). We directly measured the ratios between the corresponding fluxes to get a first estimation of the indices. We refer to these indices measurements as * and D4000*. To be able to compare these values with the typical index definitions, we correct them by using the ratio between the index measured in the best-fitting model of each galaxy at the central wavelength of each filter convolved to the SHARDS resolution and the standard index measured directly in the models. The typical value of this correction is small, the average is 1.11 for , and 0.99 for D4000. We also checked that the correction is rather insensitive to the use of only one stellar population model to measure the correction, i.e., using a larger set of models with different timescales and ages has a very small effect ( 10%). We remark that the indices calculated in this way are completely independent of the SED-fitting procedure and are based on the observed photometric data alone.

In Figure 6, we plot a zoom of the SED in the spectral range covered by SHARDS for the galaxy presented in Figure 5. We also plot the corrected   and D4000 indices for that particular galaxy, as well as the tracks expected for the evolution of these indices according to the models of the 3 best-fitting (degenerate) solutions obtained for that galaxy (see Sect. 3.2). The indices measurements allow us to eliminate solutions which are incompatible with the spectro-photometric data. For this particular galaxy, the measured indices are more compatible with the second and third solutions, and we discard the most significant solution (the one with t  1.7 Gyr).

In addition to the use of spectral indices, we evaluated the robustness of the different solutions obtained for each galaxy by making use of the (un-)detection in the mid- and far-IR and energy balance arguments. For each best-fitting solution, given the dust attenuation and mass, we derived the expected flux at an observed wavelength of 24 m (F(24), for the galaxy shown in Figure 4, the values are written in the legend). The procedure starts calculating the luminosity absorbed by dust in the UV/optical, which has to be re-emitted in the IR. We assumed that the stellar emission absorbed by dust must be equal to the IR luminosity integrated from 8 to 1000 m, L. We then used the Rujopakarn et al. (2013) relation to transform the L into a 24 m observed flux (a relation that depends on redshift). We note that the Rujopakarn et al. (2013) templates are based on star-forming galaxies, which may not be comparable to the quiescent galaxies from this sample. However, we consider it a good approximation, since the contribution to the dust heating of the older stellar population is larger at wavelengths   250 m (see Bendo et al. 2012) and does not significantly affect the 24 m emission.333Since for the same L, the predicted 24 m flux for a quiescent galaxy should be lower than for a star-forming galaxy, our predicted 24 m fluxes calculated with the method described above may be overstated. This reduces the significance of our rejection of non IR detections, although we only do so for 4% of the degenerate solutions. For galaxies undetected in the IR, we were able to eliminate solutions which present expected 24 m fluxes larger than the detection limit ( 50 Jy at 50% completeness, as estimated in Pérez-González et al. 2005 to build the IR luminosity function at ). For example, for the SHARDSJ123620.3+6200844.3, the third solution predicts a 24 m flux of  55 Jy. As this galaxy is not detected in the IR, we rejected the third solution and chose as the best solution the second one (  2.6 Gyr, 300 Myr). The second solution has F(24)  51, which is still larger than the 50% completeness limit, but only by 1%. With this method, we were able to reject only  4 of the degenerate solutions, but it was a strong argument to limit the dust attenuation range probed by our stellar population analysis to values within  mag. Larger values would imply IR detections for galaxies as massive as ours.

We are able to break the degeneracies making use of the spectral indices for   32% of the galaxies with more than one cluster of solutions. For the remaining objects, the indices uncertainties are too large (due to photometric errors) or the indices predicted by the tracks of the different solutions are consistent with the measured values. The energy balance argument helps breaking the degeneracies in 4% of the cases. When all the solutions of a galaxy are compatible with the measured indices and do not violate the energy balance argument, we choose the most significant one, i.e., the most populated cluster ( 32% of the cases). For 22% of the galaxies, the solutions were very similar in , with only significant differences in . We also discarded 10% of the solutions for being unrealistic ( 10 Myr and log(M/M) 11.0, which would imply SFR  10000 Myr or larger).

In summary, out of the complete sample of 104 galaxies, 46% of them had only one possible solution (i.e., no degeneracies). Out of the 54% of galaxies with degenerate solutions, we were able to break the degeneracies either by measuring indices or by using the Mid-IR/Far-IR data for 20% of the sample. For the remaining  34%, we used the most significant solution. In the following sections, we will only consider one solution for each galaxy and we will refer to them as primary solutions. In only 12% of the cases the primary solutions are not the most significant ones. The properties of the primary solutions for each galaxy are given in Table 2.

ID sSFR log(M/M) t A Z/ Z SFR Sign.
[Gyr] [Gyr] [Gyr] [Myr] [mag] [Myr] %
SHARDS123737.94+621309.0 1.2410 (s) 1.34 0.09 1.95 0.09 0.10 0.06 10.78 0.05 2.1 1.7 199 0.00 0.02 2.5 0.08 100
SHARDS123657.46+621451.2 1.2534 (s) 1.20 0.08 2.00 0.07 0.02 0.01 11.09 0.05 2.4 1.8 316 0.66 0.07 0.4 1.40 100
SHARDS123723.91+621520.7 1.39 (p) 1.0 0.2 1.7 0.2 0.1 0.1 10.0 0.1 0.8 0.8 14 0.2 0.3 2.5   10 100
Table 2: Galaxy properties for the sample of 104 galaxies here presented: ID, redshift (; for spectroscopic, for photometric), UVJ colours, sSFR used in the sample selection (see Sect. 2.2), mass (M), best-fitting age (t), mass-weighted age (), star formation timescale (), dust attenuation (A), metallicity (Z), SFR from SED-fitting (SFR) and significance of the primary solution. The full version of the table is available online in the supplementary files.

4 Analysis of The Star Formation Histories of Massive Quiescent Galaxies

In this section, we analyse the stellar populations properties of MQGs at . For the discussion, we only consider the primary solutions identified with the methodology described in Section 3.

4.1 Statistical properties of the stellar populations of MQGs

Figure 7: Mass-weighted ages versus stellar mass (left panel) and star formation timescale (right panel). Galaxies are colour coded by star formation timescale and mass, respectively. Only the primary solutions (those selected with the different methods described in Sect.s 3.2 and 3.3) are plotted. Circles represent galaxies with only one solution, while triangles represent galaxies with degenerate solutions. The light red areas represent the region between the 1 and 3 quartiles for each parameter. Diamonds represent the median values for the whole sample (yellow) and for the subsamples (colour coded according to the legend of each panel). The coloured error bars show the region encompassing 68% of the data for each subsample. Typical individual uncertainties are plotted in the lower right corner of each plot.

In Figure 7, we show the location of our galaxies in the mass-weighted age versus mass (–M) plane, colour-coded using 3 bins in star formation timescale. The bins have been chosen at the values were the relation changes trend, according to the right panel of the same figure, where we plot the relationship between  and , using a colour-code to distinguish 3 bins of mass. The mass bins limits are approximately the median and 3 quartile values. The population of MQGs at are dominated by “new arrivals”, i.e., galaxies with relatively young stellar populations:   2 Gyr. These galaxies younger than 2 Gyr account for 85 of the sample. We also identify a tail of old galaxies (  2 Gyr) summing up   15 of the total population. Hereafter we divide the MQGs in three sub-samples depending on their  values:  mature (  1.0 Gyr, 38 galaxies), intermediate (=1.0-2.0 Gyr, 50 galaxies), and senior ( 2.0 Gyr, 16 galaxies). The statistical properties of the galaxies divided in age and mass bins are listed in Table 3 and 4, respectively.

The population of mature galaxies presents an average mass-weighted age =0.8 Gyr and relatively short timescales =60 Myr. The typical mass is log(M/M)=10.4. On the other hand, the senior population presents average values of =2.6 Gyr, longer star formation timescales, =400 Myr, and larger masses log(M/M)=10.7. The intermediate population has transitional parameters: =1.4 Gyr, =200 Myr and log(M/M)=10.5.

Note that the fact that the mature population has short values is a direct consequence of our sample selection. As we are selecting quiescent galaxies (galaxies with low levels of sSFR), galaxies with longer ( 100 Myr) and young ages ( 1 Gyr) would have too high sSFR values to enter in our quiescent selection criteria. In principle, the senior galaxies might present short or long star formation timescales, i.e., we do not have any selection bias against any of those types of galaxies. The bias in our results due to our sample selection are further discussed in Sect. 4.3, together with the average SFHs of galaxies divided in mass and age bins.

With respect to the dust attenuation, the senior galaxies are less dusty, =0.4 mag, than the mature ones, =0.8 mag. The attenuation values are in agreement with those found by other authors (e.g., Belli et al. 2015). The galaxies with the highest dust attenuations are among the youngest ( 90 of the galaxies with A  0.8 mag are younger than 2 Gyr), suggesting that they could be recently quenched galaxies.

With respect to the metallicities, 41 of the galaxies are best-fitted with solar metallicity, 36 with super-solar metallicity and 23 with sub-solar metallicity. No clear metallicity trends are found with respect to mass or age.

Figure 8: Mass-weighted ages versus star formation timescale. Galaxies are colour coded by their sSFR. Circles represent -selected galaxies, while the diamonds are sSFR-selected. The sizes of the symbols represent the percentage of mass assembled in the last Gyr. The grey areas represent the region between the 1 and 3 quartiles. Typical uncertainties are plotted in the lower right corner.

At this point, we want to mention the difference between old and quiescence when using -models. In single stellar population (SSP) models, all the star formation takes place at the same time (t), and then the population passively evolves. Therefore, an old galaxy is always quiescent and, in fact, a young galaxy would also be literally quiescent (although it may be not selected with our criteria). When using models, the duration of the star formation depends on . Galaxies with long star formation timescales may have been formed more than 1-2 Gyr ago but could still be forming stars. What actually indicates the quiescence of a galaxy (and consequently the star formation activity) is the sSFR or the relation.

In Figure 8, we show again the  versus plot, but now the galaxies are colour-coded according to their sSFR (averaged over the last 100 Myr from the best-fitting models, i.e., these sSFRs are different from those used in the selection, which were based on the observed UV luminosity and dust attenuation estimates – see Section 2.2). The galaxies with similar sSFR values are distributed along diagonal lines of constant / values. The senior population of galaxies with long values show higher sSFR values than mature galaxies with very short . The symbol size represents the percentage of mass assembled in the last Gyr. The mature galaxies ( 1 Gyr) assembled most of their mass during that time (by definition of mature objects). There is a population of intermediate galaxies with    100 Myr that has assembled between 10-50% during the last Gyr, while the senior population has assembled less than 10% of its mass. We also differentiate in Figure 8 between -selected and sSFR-selected galaxies. sSFR-selected objects, i.e., galaxies outside the quiescent -region (or inside the -quiescent region but detected in the IR, see Sect. 2.2), are preferentially located in the lower right corner of the - plane. This indicates that the sSFR-selected galaxies have larger / values and are less quiescent than the pure galaxies. In Table 5 we show the median properties of each sub-sample and we discuss the distribution of the derived galaxy properties in the -plot in Sect. 4.2 and Figure 11.

Parameter Sample Number Q1 Median Q3
Mature N=38 1.12 1.23 1.33
Interm. N=50 1.02 1.17 1.24
Senior N=16 1.02 1.06 1.19
t [Gyr] Mature 0.8 0.9 1.0
Interm. 1.3 1.8 2.2
Senior 3.0 3.2 4.2
[Myr] Mature 20 60 100
Interm. 30 200 300
Senior 300 400 500
A [mag] Mature 0.5 0.8 1.1
Interm. 0.1 0.5 0.8
Senior 0.1 0.4 0.6
Z/ Z Mature 1.0 1.0 2.5
Interm. 0.4 1.0 2.5
Senior 0.4 1.0 1.0
log(M/M) Mature 10.3 10.4 10.7
Interm. 10.3 10.5 10.8
Senior 10.6 10.7 11.0
 [Gyr] Mature 0.7 0.8 0.9
Interm. 1.2 1.4 1.7
Senior 2.2 2.6 3.4
sSFR [Gyr] Mature   10   10 10
Interm.   10 10 10
Senior 10 10 10
t [Gyr] Mature 0.2 0.5 0.7
Interm. 0.7 0.9 1.1
Senior 1.5 1.9 2.6
Table 3: Galaxy properties of our sample in three age bins: mature (  1.0 Gyr), intermediate (1.0     2.0 Gyr) and senior (  2.0 Gyr). We show the 1, median and 3 quartiles for redshift (), the best-fitting ages (t), star formation timescales (), dust attenuations (A), metallicity (Z), masses(M), mass-weighted ages (), sSFR and time since quenching (t, explained in Sect. 4.2)


Parameter Sample Number Q1 Median Q3
Low M N=45 1.13 1.24 1.33
Interm. M N=36 1.02 1.14 1.23
High M N=23 1.01 1.12 1.24
t [Gyr] Low M 0.9 1.1 1.6
Interm. M 1.0 2.0 2.3
High M 1.0 2.0 3.0
[Myr] Low M 20 80 200
Interm. M 30 200 400
High M 30 250 300
A [mag] Low Mass 0.2 0.6 1.0
Interm. M 0.3 0.6 0.9
High M 0.3 0.5 0.8
Z/ Z Low M 1.0 1.0 2.5
Interm. M 1.0 1.0 2.5
High M 1.0 1.0 1.0
log(M/M) Low M 10.1 10.3 10.4
Interm. M 10.6 10.7 10.8
High M 10.8 10.9 11.0
 [Gyr] Low M 0.8 1.1 1.3
Interm. M 1.0 1.4 1.8
High M 1.0 1.6 2.4
sSFR [Gyr] Low M   10 10 10
Interm. M   10 10 10
High M   10 10 10
t [Gyr] Low M 0.4 0.6 0.9
Interm. M 0.7 0.9 1.2
High M 0.8 1.0 1.7
Table 4: Galaxy properties of our sample in 3 mass bins: log(M/M)=[10.0,10.5], log(M/M)=[10.5,10.8] and log(M/M)=[10.8,11.4]. The parameters are the same as those in Table 3.


Figure 9: Differences in the average stellar population models for the mature and senior galaxy sub-samples (see text for details). Upper panel: stellar population models with the average properties of the mature and senior populations (as shown in the legend) in the spectral region covered by SHARDS. The models have been normalized to the fluxes in the 230–450 nm range. We show both the original resolution of the models (thin lines) and the models convolved to the SHARDS resolution (R=50, thick lines). Lower panel: ratio of the two models shown above at the original model resolution (thin grey line) and at R50 (thick black line). We remark the significant (and measurable) differences between the two models at the resolution of our spectro-photometric data, especially for the the   and D4000 indices (marked in both plots as coloured areas).

Our result about a duality (mature versus senior systems) in the population of MQGs at and their differences not only in ages but also in timescales and dust attenuation was tested. The test was planned to check that the population of old (2 Gyr) MQGs at , which is just a 15% and a tail in the age distribution, is real and its existence is not an effect of the degeneracies (in the - plane). For this purpose, we constructed spectro-photometric stacks for the mature (  1.0 Gyr), intermediate (=1.0–2.0 Gyr) and senior (  2.0 Gyr) populations.

The broad-band colours of the three sub-samples of galaxies are very similar (they were selected with the same method based on colours) in the whole wavelength range. However, at the resolution achieved by the spectro-photometric data from SHARDS and the WFC3 grisms, especially in the   and D4000 spectral regions, the sub-samples and the average fitting models show measurable differences. In Figure 9, we show the stellar population models with the average properties of the mature and senior populations, concentrating in the spectral region covered by SHARDS. There are significant differences in the relative fluxes for the two models. In the D4000 region there is an excess of flux for the mature population model and the D4000 break is stronger for the senior population model. This translates to a   20% difference in the relative fluxes of the two models in the blue band of the D4000 index. In the  region the relative flux of the mature model is even higher, up to 70%, revealing the signature of a younger population, which is less significant in the senior model. We should be able to distinguish between the two solutions by measuring the differences in the indices, and this is possible in the SHARDS and grism data (and not with broad-band observations).

Figure 10: Left panel: Stacks for the senior (thick red line) and mature (thick blue line) populations normalized in the   (upper panel) and D4000 regions (lower panel). The average error of the stacks is  10%. We used the dispersion as uncertainties (actual errors calculated by propagating the observed flux uncertainties were negligible). We also plot the normalized models with the average properties of each sub-population (thin red and blue lines). Right panel:  versus D4000 plane, with the indices measured in the stacks for the mature, intermediate, and senior galaxies (blue, green and red diamonds with error bars). For comparison, we also show the evolutionary tracks for 2 models with different star formation timescales (a SSP and a delayed exponential model with 1000 Myr) and two different dust attenuations (A=0.0 mag and A=1.0 mag). The empty circles represent the expected index values for each track at different ages (given in Gyr in the plot). The three coloured lines show the tracks predicted by the models with the average properties of each galaxy population (given in the legend). The filled coloured circles mark the indices values at different ages (0.5, 1.0, 2.0, 3.0, 5.0 Gyr, from smaller to larger symbols). The indices measured in the stacks are compatible with the average properties derived from the SED-fitting for each sub-population.

In Figure 10, we present the stacks built for the mature and senior galaxy populations together with the models characterized by their average properties. We also show the indices measured in the stacks of the mature, intermediate and senior galaxy populations, and the tracks predicted by the model characterized by the average timescales and dust attenuations quoted in Table 3 for each sub-population. The indices measured directly from the stacks are consistent with those predicted by the tracks of the models with the average properties of each sub-population. We demonstrate that the spectro-photometric data directly show the differences in the stellar population properties for the 3 sub-samples. Note that the origin of the average stellar population properties is the full-SED spectral fits, and here we are comparing those models with finer resolution data in a limited wavelength range. The two analyses are not completely independent, but the information encoded in the data at the resolution which makes index measurements possible could easily be erased or degraded by fitting the whole SED from the UV to the mid-IR. Our test demonstrates that we are seeing differences in the galaxy populations in both the global SED and the spectral indices. We, thus, conclude that our results about the ages of MQGs at are robust and not an artefact linked to the SED-fitting degeneracies.

The longer timescales for the senior galaxies, which are also relatively massive, is not directly consistent with a typical conception in the downsizing scenario which states that the most massive galaxies formed their stars early but also very rapidly. A short timescale for the formation of massive galaxies has also been claimed to be necessary to explain the -element enhancement seen in early-type massive galaxies (such as ellipticals) in the nearby Universe (see, e.g., Thomas et al. 2005 and Renzini 2006). Our estimations of the timescales for the most massive galaxies are not extremely long (typically 400 Myr) and would be rightly consistent with the values needed to match the chemical abundances (less than 1 Gyr; see Thomas et al. 1999; Worthey et al. 1992; Worthey 1998 and references therein). In addition, these massive galaxies may be the product of mergers (or clumps) involving progenitors with shorter timescales, but with different ages, maybe linked to small offsets (of the order of tens or a few hundred Myr) in the ignition of the star formation. The SFHs obtained for the merged systems (i.e., our galaxies) would then present longer timescales, but do not violate any constrain linked to -element enhancement. In Sect. 4.3, we further discuss the SFH of our galaxies as a function of the mass and age and the importance of selection effects.

The variety in the ages of quenched galaxies at is in agreement with previous works such as Bedregal et al. (2013) and Belli et al. (2015). Both papers present a wide range of stellar population properties for similar samples as ours. Bedregal et al. (2013), analysing the G102-G141 WFC3 grisms of 41 massive (log(M/M 10.65 M) quiescent galaxies at , with exponentially declining SFH found that they had short star formation timescales ( 100 Myr) and a wide distribution in stellar ages (1–4 Gyr). Belli et al. (2015) analysed Keck LRIS spectra of quiescent galaxies at with log (M/M 10.6 M and concluded that there are two different quenching routes. The youngest galaxies (t  1 Gyr) show a fast quenching (  100 Myr), while the older galaxies (up to 4 Gyr old) show slowly-declining SFHs ( 200 Myr). Note, however, that our values and theirs are not directly comparable, since we used different parameterizations for the SFH (i.e., the parameter for an exponential is not exactly the same as the value for a delayed exponential).

4.2 Dissecting the diagram: distribution of stellar population properties

Figure 11: diagrams with galaxies colour-coded on the basis of different properties, from left to right, top to bottom: mass-weighted age (), star formation timescale (), dust attenuation (A), metallicity (Z), SED-based sSFR (sSFR), and time since quenching (t). The circles represent galaxies with no degenerate solutions, while the triangles are galaxies with more than one solution in the Monte Carlo simulations (see Sect. 3.2). -selected galaxies are located inside the quiescent region of each diagram, while sSFR-selected galaxies fall outside by definition (see sect. 2.2, Table 5).
Parameter Sample Number Q1 Median Q3
N=65 1.05 1.17 1.25
sSFR N=39 1.03 1.15 1.24
t [Gyr] 0.9 1.1 2.0
sSFR 1.0 1.7 2.5
[Myr] 20 60 200
sSFR 160 250 400
A [mag] 0.1 0.5 0.8
sSFR 0.4 0.7 1.1
Z/ Z 1.0 1.0 2.5
sSFR 0.4 1.0 2.5
log(M/M) 10.3 10.5 10.8
sSFR 10.3 10.7 10.8
 [Gyr] 0.9 1.1 1.7
sSFR 0.8 1.2 1.7
sSFR [Gyr]   10   10 10
sSFR 10 10 10
t [Gyr] 0.7 0.9 1.2
sSFR 0.3 0.6 1.0
Table 5: Galaxy properties of our sample according to their selection criteria (/sSFR, see Sect. 2.2). The parameters are the same as those in Table 3.

In Figure 11, we plot the colour diagram, which was the starting point of our sample selection, to study how the colours do actually correlate with the derived galaxy properties. The main properties of the - and sSFR-selected sub-samples are shown in Table 5. We recall that we refer to galaxies in the quiescent region of the diagram and without IR detection as -selected, while the sSFR-selected have sSFR  0.2 Gyr but are located outside the passive region (or inside but detected in the IR, see Sect. 2.2).

Concerning ages, there is not a significant difference between those derived for the -selected and sSFR-selected samples. The average ages are =1.1 and 1.2 Gyr for each sub-sample, respectively. However, the senior galaxies seem to have redder colours, with average values V-J=1.46, than the mature population, V-J=1.31.

The timescales are more clearly segregated in the plot. The of the -selected galaxies are shorter than those of the sSFR-selected galaxies: typical values are 60 Myr and 250 Myr, respectively. Comparing our results with those in Belli et al. (2015), who also presented stellar population properties within the diagrams, we find a similar trend: they found that galaxies with   100 Myr were distributed in a narrow region parallel to the diagonal line defined in the passive region.

Belli et al. (2015) also found that galaxies with higher dust attenuations (A  1 mag) were outside the passive region. We find typical A=0.5 and A=0.7 mag for the -selected and sSFR-selected galaxies, respectively, in good agreement with Belli et al. (2015). The percentage of -selected galaxies with A  1 mag is only 7, while this percentage increases up to 28 for the sSFR-selected galaxies. All of the galaxies falling in the post-starburst region defined by Whitaker et al. (2011) (see Fig. 2) have A  0.5, i.e., they are relatively dust-free. However, we do not see such a clear gradient in A with the distance to the division line as in Belli et al. (2015). A possible explanation for this discrepancy is the difference in the metallicity values used in each work. While Belli et al. (2015) used a very narrow range in metallicity (a normal prior centered on the solar value Z=0.02 with a width of 0.005), we use 3 discrete values (Z=0.02, 0.05 and 0.008). Due to the dust attenuation-metallicity degeneracy, the dust attenuation distribution is obviously affected by the assumed metallicity values.

With respect to the metallicity, no clear trends are present for any of the two selection criteria. The fraction of -selected galaxies with sub-solar, solar and super-solar metallicity is 20, 40 and 40%, respectively. For the sSFR-selected, the percentages are 25, 46 and 30%. We note that, if the metallicity was a fixed parameter (e.g., solar metallicity), the galaxy properties would be better segregated within the diagram, as the best-fitting model would be constrained by only 3 parameters instead of 4.

In Figure 11, we also present how SED-based sSFRs correlate with the position in the diagram. The most quiescent galaxies are located in a similar region to the galaxies with low values. In fact,  60 of the -selected galaxies have sSFR  10 Gyr, while this only happens for  8% of the sSFR-selected galaxies.

We have also derived the time since quenching, defined as the time since the galaxy became quiescent using our definition from Sect. 2.2, i.e., how much time has passed since the galaxy had sSFR  0.2 Gyr. The -selected galaxies have been dead, on average, for almost 1 Gyr, =0.9 Gyr, while the sSFR-selected galaxies have been more recently quenched, =0.6 Gyr.

4.3 Tracing back the SFHs of MQGs: clues about their past

Figure 12: Evolutionary tracks in the SFR-Mass plane for our sample of MQGs. Left panel: the small coloured dots mark the location of galaxies at different times after their formation, as indicated in the legend. Note that these are direct ages from the SED-fitting and that we have assumed the same origin for all galaxies. We mark as thin grey lines the past path for each galaxy. The thick grey lines show the location of the MS at   1.2 (dashed-dotted line) and   2 (solid line) according to Speagle et al. (2014), and at   1 according to Elbaz et al. (2007) (dashed line). The green shadow is the 1 dispersion of the MS at z=2.0. Most of the galaxies are   0.5 – 1.0 Gyr old when they are on the MS region at   1.0. In fact, 75% of our sample is 1 below the MS at z=1.2 when they are 1 Gyr old and all of our galaxies are at least 1 below the MS, considering the MS evolution from  2 to  1, after 2 Gyr. Right panel: Location of our galaxies in the SFR-Mass plane at the epoch of observation (large filled circles), colour-coded by their best-fitting ages (t). We note that both the ages and SFRs used in this plot are the output of the SED-fitting (and not SFR neither ). To better visualize the location of the galaxies in the SFR-Mass plane, we plot galaxies with very low SFR ( 10M yr) around SFR  10M yr.

In this Section we analyse the SFHs for MQGs at in the observed redshift range and at earlier lookback times. One of the advantages of determining SFHs with stellar population models is that they provide us with a time-dependent evolution of parameters such as the SFR and the galaxy stellar mass. Therefore, assuming closed-box evolution (i.e., no merging events or re-activation of the star formation activity), we do not only characterize the properties of galaxies at the epoch of observation, but we can also trace back their properties and see how they assembled their mass at earlier cosmic times. At this point, we want to stress that the results hereafter presented are strongly affected by the assumption of the SFH as a single burst, delayed exponentially declining model. This SFH parameterization can account for the star formation taking place after one gas-rich major merger, but neglects the stellar population of the mass previously assembled (see Figure 1 from Hopkins et al. 2008). CDM models and observed merger rates (e.g. Bluck et al. 2009; Newman 2013; Man et al. 2014) predict that massive galaxies undergo at least one major merger since =3. However, assuming multiple bursts of star formation, would imply the analysis of two stellar population models. This would significantly increase the degeneracies in the SED-fitting procedure and would severely complicate the interpretation of the results. We discuss the impact in our results of a more complicated SFH in Appendix A, but we delay a comprehensive analysis of more complicated SFH parameterizations to future papers.

4.3.1 The SFR-Mass relation

One of the most fundamental relations between galaxy properties is the SFR-Mass relation, the so-called MS. In Figure 12 we show the evolutionary paths of the galaxies in our sample in the SFR-Mass plane as a function of time. For this plot, we assumed that all galaxies started their formation at the same time (i.e., all SFHs share the t0 point). The shape of the track depends mainly on the final mass and the star formation timescale, : galaxies with short timescales form the bulk of their mass very quickly, they arrive relatively early to the MS and then their star formation decreases while their stellar mass remains almost unchanged (falling vertically). Galaxies with very long values continuously increase their mass without changing their SFR (almost horizontal tracks) and may stay in the MS up to   2 Gyr. According to our best-fitting SFHs (note that the choice of the parameterization is also important), after 100-500 Myr of evolution our galaxies would form a MS with a very similar slope to that observed directly at with samples of star-forming galaxies (see, e.g., Speagle et al. 2014). We remark that, for this exercise, we offset the SFHs of our sample to make them match at t0. The effect of galaxies starting their formation at different epochs would result in a widening of the MS shown in Figure 12 for different ages. According to our results, galaxies would come out of the MS (1 below) after approximately 1.5 Gyr (96% of our MQGs are located below the MS at that age), and quenching would then proceed almost vertically in this plot. We remark that this statement considers the main sequence evolution from  2 to  1. Therefore, we find that MQGs would pass a fair fraction of their lifetimes in the MS, with the possibility to live above the MS (in the starburst locus) for short periods of time of the order of 100 Myr. We conclude that the SFHs determined for the MQGs at are consistent with the slope and even the location of the MS at and that the existence of the MS for our sample of MQGs is mainly an age effect, in the sense that the MS is formed by galaxies with similar ages ( 0.5-1 Gyr).

We also want to remark that the current SFR for many of our galaxies derived from the SED modeling (averaged over the last 100 Myr of their history) is much smaller than the SFR obtained from observables. The current SFR are not fully consistent with the SFR and SFR measurements explained in Section 2.1 and used in the sample selection. Indeed, 45% of the sample have SFR  10M yr but the lowest SFR estimated from the typical tracers is 10M yr. This could be due to an overestimation of the dust attenuation from the IRX- relation, but the choice of a delayed exponentially decreasing SFH is most probably responsible for this effect. Typically, this parameterization produces very low SFRs for dead galaxies (whose emission is dominated by an evolved stellar population), but there might exist a second stellar population with some (very low and negligible in terms of mass and emission) on-going star formation which would not be possible to recover with the assumed SFH parameterization. However, the fits of the model to our data are very good for all the sample and the inclusion of a second population would multiply by 2 the uncertainties and degeneracies. We, therefore, assume that our galaxies are dominated by the older stellar population and that one (composite) stellar population model explains the main features of our galaxies. We caution the reader, however, that the current SFR should be taken as lower limits.

Following Kennicutt (1998), we can derive the SFR corresponding to LIRGs/ULIRGs: L=10 L corresponds to 12 M yr, and L=10 L to 121 M yr (transformed to Kroupa 2001 IMF). All of the galaxies in our sample had SFR peaks larger than the LIRG limit (except one which has log(M/M)= 10.0). The typical fraction of their lifetime spent in the LIRG phase is  32 ( 500 Myr). A significant fraction of the sample (46) had SFR peaks larger than the ULIRG limit, but the typical fraction of their lifetime in this phase is much smaller,  8 ( 100 Myr). Our results favor LIRGs at as the most likely progenitors of MQGs at , and that most of their stars were formed in star-forming events with SFRs around 100 M yr. The ULIRGs seem to be the progenitors of the most massive galaxies. The fraction of galaxies which have undergone a ULIRG phase increases up to 65% (75%) when considering galaxies with masses larger than 10 M (10 M). We discuss the connection between MQGs and ULIRGs on the basis of their number densities in Sect. 4.3.2.

Figure 13: Evolution of the MQGs in the SFR-Mass plane at different redshifts. The observed properties (those measured directly from the data) are plotted as black filled circles in the redshift bin . The expected location of the galaxies 1 Gyr () and 2 Gyr () before the observations are plotted as black filled circles in the upper right and upper left panels, respectively. The lower right panel represents the expected location of the galaxies 1 Gyr after the observations (). We use the age, timescale, mass and SFR given by our SED fits to move forward and backwards in time in the SFR-Mass plane, assuming passive evolution. For comparison, we plot as a green line the MS from Speagle et al. (2014). The dark and light green areas mark the 1 and 2 scatter, 0.2 and 0.4 dex. The blue line in the panel is the MS from Elbaz et al. (2007). The orange/red lines show the SFR limit for LIRGs/ULIRGs converted into SFRs using Kennicutt (1998) relation. The light blue dotted line is the value of constant sSFR used in the sample selection (see Sect. 2.2; sSFR=0.2 Gyr). Galaxies with SFR  10M yr are plotted as black empty diamonds at SFR  10M yr.
1.7 – 3.0 N=2 N=7 N=8 N=15 N=12
7% 23% 27% 50% 40%
1.3 – 2.0 N=8 N=10 N=5 N=52 N=46
12% 15% 7% 78% 69%
1.0 – 1.5 N=104 N=101
100% 97%
0.8 – 1.2 N=104 N=104
100% 100%
Table 6: Number and percentage of galaxies above (1 or 2), in and below (1 or 2) the MS from Speagle et al. (2014) at each redshift bin, as show in Figure 13. Note that the sum of percentages at each redshift does not equal 100 as the  1 bins also include the galaxies at  2.

In Figure 13, we show where our galaxies would be placed in the SFR-Mass plane at different redshifts according to their past and future evolutionary tracks. We study 4 redshift bins corresponding to 2 and 1 Gyr before the observations, the actual observed redshift and 1 Gyr after the observations. In Table 6, we show the percentages of galaxies above, within, and below the MS at different redshifts. Figure 13 shows that at , the epoch of observation, all of our galaxies lie 1 below the MS (this happens also when using the SFRs based on classical tracers), as expected given that we started with a selection of passive galaxies. If we move 1 Gyr backwards in time,  22% of the galaxies are located in or above the MS, but a significant number fraction ( 78%) of the galaxies at are located 1 below the MS, meaning that they were already quiescent or in the process of quenching. At , the bulk of the galaxies cannot be plotted because  70 have best-fitting ages smaller than 2 Gyr, so at that epoch their masses were much smaller than our limits in the plot (or they were not even formed yet). At that redshift bin, half of the galaxies (50%) are in the MS or above but the other half (15 galaxies) are 1 below the MS as early as   2. In the opposite time direction, if we move 1 Gyr after the observations, all of our galaxies have very low SFR and are completely dead and well below the MS from Speagle et al. (2014) and Elbaz et al. (2007), as expected from purely passive evolution.

Comparing with the parent population of galaxies more massive than 10 M at each redshift bin, we find fractions of quiescent galaxies of 25, 12 and 3 at , and , respectively. Different studies have shown that the fraction of quiescent galaxies constantly increases with time. For example, Domínguez Sánchez et al. (2011) derived a fraction of quiescent galaxies (log (M/M 10.6, sSFR  0.01 Gyr) of  20 and 10 at and respectively. The percentage found in Ilbert et al. (2013) for a sample of quiescent galaxies selected by their NUV-r vs r-J colours with log(M/M 9.6 increases from 6 to 13 from to . The evolution of the fraction of quiescent galaxies ( selected) in Muzzin et al. (2013) is not so strong (28 at and 24 at ), although the mass limits used in this calculation change for each redshift bin, log(M/M 9.48 at , log(M/M 10.54 at . The fractions of MQGs that we derive at z  1.5 by studying the past evolution of MQGs at z1.0-1.5 predicted from their SFHs are consistent with observational results at higher redshifts.

4.3.2 Number densities

Figure 14: Number density of MQGs galaxies at  1.0 – 1.5 as a function of their mass-weighted age  (red stars). The large black star is the total number density of MQGs from our work. We also show the number density of quiescent galaxies at the same redshift from Muzzin et al. (2013) (grey area). The small empty orange triangle depicts the number density and mean age of MQGs at = 1.4 – 2.2, as measured by Whitaker et al. (2013). Assuming passive evolution for these galaxies, they would move to the orange filled triangle. The number density uncertainties (shown as 2) have been calculated assuming Poisson statistics, while the error bars in age represent the standard deviation in each age bin.

In the next paragraphs, we compare the number density of quiescent galaxies in our work with a compilation of estimations from the literature. The number density of quiescent galaxies in our work for are derived considering our results for the SFHs of MQGs and the position of the galaxies in a SFR-Mass plot at different epochs (Figure 13). In this case, we assume that galaxies evolve passively once their star formation is quenched and they do not experience any merger event. We warn the reader that a fair comparison of number densities requires taking into account the differences in quiescent fraction linked to the stellar mass cut or the definition of quiescence, which vary from paper to paper and are difficult to consider in their full extent. Here we consider a galaxy as quiescent when is located 1 below the MS at each redshift. With this method, the number densities of MQGs from our work are =(7.0 0.7)10 Mpc at , =(2.3 0.3)10 Mpc , and =(0.31 0.08)10 Mpc . This is in good agreement with the number densities of quiescent galaxies reported in Muzzin et al. (2013) at (=7.610 Mpc), but smaller than their number densities at higher redshifts =(3.3, 1.2, 0.65)10 Mpc at , and , respectively. There are several factors that could be affecting this comparison. In the first place, the redshift bins are not exactly the same. Besides, Muzzin et al. (2013) selection is only based on colours, while we are considering as quiescent galaxies 1 below the MS at each redshift, as derived from the backwards evolution of the MQGs at . We recall that we have eliminated from our MQGs at galaxies detected in the IR in the quiescent region, when they have sSFR  0.2 Gyr. We also note that galaxies in the quiescent region of the diagram at can have much larger sSFR (up to 1.0 Gyr) than our quiescent limit. We recall that the values derived for the past evolution of MQGs must be taken with care, as we are assuming pure passive evolution and no merger events (which could help to decrease the number density of quiescent sources).

We also derive the number densities of our observed galaxies (=1.0 – 1.5) for the different populations defined in Sect. 4.1 (mature, intermediate, senior). This is shown in Figure 14, together with the number density of Muzzin et al. (2013) at =1.0 – 1.5 and Whitaker et al. (2013) at =1.4 – 2.2. The existence of a population of old (t  1.3 Gyr) galaxies at z  2 has already been found by Whitaker et al. (2013) after analysing the stellar populations of 171 massive -selected galaxies using stacked 3D-HST spectra. Given that in this work they use SSP models and we use more general (and realistic) delayed exponential models, the comparison with our sample is not straight forward, but we compare their results with our mass-weighted ages. Taking into account that there is a difference of  1.6 Gyr between the median redshift of Whitaker et al. (2013) sample,  = 1.64, and the median redshift of the senior population,  = 1.06, the galaxies observed by Whitaker et al. (2013) would be  3 Gyr old at the redshifts probed in our work, should they stay passive and not re-ignite. In our sample, we find a number density of galaxies with   2.0 of (1.1 0.2) 10 Mpc for log(M/M 10.0 and (0.9 0.2) 10 Mpc for log(M/M 10.5. These results are consistent within errors with the number density of galaxies from Whitaker et al. (2013): (0.8 0.1) 10 Mpc at and with log(M/M 10.50. The similar number densities from Whitaker et al. (2013) for MQGs with ages around 1.3 Gyr at z  2 and those obtained in this work with ages  2 Gyr at z   1.1 is a good consistency check on the existence of old galaxies at . However, this value is smaller than the number density of Muzzin et al. (2013) at for galaxies with log(M/M 10.50 (  1.910 Mpc). This difference could be due to the more restrictive selection from Whitaker et al. (2013), where the authors require G141 WFC3/HST grism detection for their sample. Given the discrepancies on previous results and the difficulty in making a fair comparison between number densities, it is difficult to reach firm conclusions regarding the assembly of the red sequence. However, the number densities (and fractions of quiescent galaxies) that we find when considering passive galaxies at z1.0 – 1.5 and moving back in time seem to be consistent or smaller than those reported in studies based on galaxies at those actual redshifts above z1.5. This could suggest that mergers play an important role or that the SFH of galaxies could be more complicated than the delayed exponentially declining assumed in this work (i.e., with the re-ignition of the star-formation for some galaxies).

To check the possibility that ULIRGs are the progenitors of MQGs, as mentioned in Section 4.3.1, we compare them in terms of number densities. Magnelli et al. (2013) found a number density   1.010 Mpc for ULIRGS at  = 2. Taking into account that 46% of our galaxies may have undergone a ULIRG phase, this would mean a number density of   3.210 Mpc. However, the duty cycle of the ULIRG phase is very short ( 100 Myr, see Sect. 4.3.1), which would significantly reduce the observed number density of ULIRGs at z=2. We conclude that the the possibility that ULIRGs  = 2 are the progenitors of the MQGs  = 1.2 is consistent in terms of number densities and timescales in rough terms.

4.3.3 Age evolution

In Figure 15 we show the evolution of the ages of MQGs over the last 10 Gyr. We compare our results with previous works based on stacked spectra (Mendel et al., 2015; Whitaker et al., 2013; Onodera et al., 2015; Schiavon et al., 2006; Choi et al., 2014) and individual spectral measurements (Toft et al., 2012; Krogager et al., 2013; van de Sande et al., 2013; Bedregal et al., 2013; Belli et al., 2015; Barro et al., 2015). We compare the data with the predictions from pure passive evolution of a delayed -models with different formation redshifts. The median ages derived for our sample (  1.2 Gyr) are consistent with those from Belli et al. (2015) at the same redshift. The ages derived for the senior population could be explained in terms of passive evolution of the galaxies studied at higher redshifts in Toft et al. (2012); Krogager et al. (2013); van de Sande et al. (2013); Whitaker et al. (2013); Onodera et al. (2015); Mendel et al. (2015) and Barro et al. (2015), suggesting that at least part of the quiescent population at   1.5 does not restart the star formation once they are quenched. However, the number density of old (age  2 Gyr) MQGs is small and the average properties of MQGs at   1.2 are dominated by galaxies with age  2 Gyr.

The observed properties of galaxies are not fully consistent with a pure passive evolution of our sample of MQGs at   1.2. If the mature population evolved passively, they would have typical ages consistent with those derived for the most massive galaxies (log (M/M 11.3) of the Choi et al. (2014) sample at   1. However, the median mass of the mature population is  10.4, meaning that they should grow by almost 0.9 dex in mass in 3-5 Gyr to be consistent with Choi et al. (2014) results. Recent works suggest that galaxies can increase their mass by a factor of  2 due to residual SFR (e.g. Pérez-González et al. 2008; Fumagalli et al. 2014) or via major mergers (Shankar et al., 2015), which is not enough to account for this mass discrepancy. The masses of the senior population (log(M/M)=10.7) would be consistent with the less massive or intermediate mass sample of Choi et al. (2014), suggesting that they could be the progenitors of MQGs. In this case, some level of rejuvenation should take place since 1.2, in the form of residual star formation or via mergers of senior and mature (or even star forming) galaxies, in order to reconcile the ages of our sample with those from Choi et al. (2014). For example,  50% of the total mass should be formed in a burst of  0.5 Gyr to make the ages of our senior galaxies (=2.6 at =1.06) consistent with the ages observed by Choi et al. (2014) (=2.5 Gyr at =0.6) for the intermediate mass population. The resulting total mass would be log (M/M)=11.0, also consistent with the values observed by Choi et al. (2014). SFH with longer star formation episodes (  400 Myr) could also help alleviating the discrepancies. The properties of the less massive (log (M/M 10.7) quiescent galaxies at from Choi et al. (2014) cannot be explained in terms of purely passive evolution of quiescent galaxies at higher-, suggesting that the less massive population is still forming stars at . However, the lack of number density considerations complicates the comparison with Choi et al. (2014), and may explain some of the discrepancies.

The relatively uniform ages (around 1-2 Gyr) measured at  1 suggest that the quiescent galaxy population is being kept young by the constant addition of recently quenched galaxies, i.e., “new arrivals”, as we mentioned in Sect. 4.1 and was already stated in previous studies (e.g. Mendel et al. 2015). We conclude that the formation of the red sequence of quiescent galaxies is actually occurring at  1.0 – 2.0 (no results are available beyond that redshift). At these redshifts, the number density of the oldest population is small and the population of dead galaxies is dominated by new arrivals with ages around 1-2 Gyr (or mature galaxies, as defined in this work). Only at redshifts below   1, the MQG population is totally assembled and evolves passively with no significant new additions. This is in agreement with previous studies supporting a significant evolution of the number density of MQGs at lower redshifts: by a factor of  2 from 1.5 down to 0.7, epoch at which the bulk of this population seems to be definitively assembled (see Daddi et al. 2005; Eliche-Moral et al. 2010; Davidzon et al. 2013; Prieto et al. 2013; Mendel et al. 2015; Prieto & Eliche-Moral 2015).

Figure 15: Evolution of the median stellar age of MQGs over the last  10 Gyr. The results from this work (red stars) are plotted for three  bins at their median redshift, the size of the symbols being proportional to the number density. The largest symbol is at   1.2 Gyr, the median age of our sample. Our results are compared with other ages from the literature at different redshifts measured from stacks using SSP models (Mendel et al. 2015; Whitaker et al. 2013; Onodera et al. 2015; Schiavon et al. 2006) or from individual spectra using -models (Toft et al., 2012; Krogager et al., 2013; van de Sande et al., 2013; Bedregal et al., 2013; Belli et al., 2015). We also plot the mean luminosity-weighted age for 3 quiescent galaxies and the luminosity-weighted age for a galaxy in the process of quenching derived in Barro et al. (2015). The results from Choi et al. (2014) at lower redshifts are divided in 3 mass bins, log(M/M 10.7, 11.0, 11.3 from lighter to darker blue. The grey lines show the age evolution of models with SFHs of the form SFR(t) t/e, with  = 100 Myr and different formation epochs () and =400 Myr and . The yellow area indicates the redshift range studied in this work.

4.3.4 Mass and age dependence of the SFH

To study the dependence of the SFH with mass and age, we divide our sample in 3 mass bins in log(M/M) units: 10.0 – 10.5, 10.5 – 10.8, and 10.8 – 11.5; and also in 3 bins in  1 Gyr, 1 – 2 Gyr, 2 Gyr (corresponding to the mature, intermediate and senior populations introduced in Sect. 4). We then construct the typical SFH for each sub-sample using the median and SFR values. The SFHs are normalized so that, by integrating them, we recover the median mass of each subsample. We also derive the mean age of the Universe when the galaxies were formed (t). We plot the results of the SFR evolution for the three galaxy sub-samples divided by age and mass in Figure 16, and we give the obtained values in Table 7.

If we study the SFH of our sample divided in age bins (left panel of Figure 16), we find that the senior galaxies must have been formed, on average, when the Universe was only 2 Gyr old, while the intermediate and mature galaxies are formed when the Universe was 3.4 and 4.1 Gyr old, respectively. The age of the Universe when the galaxies had their maximum SFR is t  2.3, 3.6, 4.2 Gyr with median SFR60, 110, 230 M yr and median 400, 200, 60 Myr for the senior, intermediate and mature galaxies, respectively. This is a direct implication of the galaxy properties derived for each population and summarized in Table 3. According to this scenario, the senior MQGs would have been formed in a relatively young Universe in longer and not very intense star formation episodes, while, on the other hand, the recently quenched MQGs at 1.0-1.5 must have been formed at later times in shorter and more intense bursts of star formation. We must notice that there are important selection effects in these results. The mature population must have been formed relatively quickly to be able to quench their star formation in less than 1 Gyr (their maximal age, by definition), meaning that they must present short values to satisfy our quiescent selection criteria. Instead, the senior population does not present any selection bias, and could, in principle, have short or long timescales.

Figure 16: Average SFHs (SFR versus time) for our sample of MQGs divided in 3 age bins (left panel) and 3 mass bins (right panel). The yellow shaded area indicates the redshift studied in this work.

With respect to the SFH divided in mass bins, the age of the Universe when they were formed was 3.6, 3.4 and 3.1 Gyr for the less massive, intermediate and most massive galaxies, respectively. The less massive galaxies are formed  0.5 Gyr later than the most massive ones. This difference in time is larger than our age uncertainties ( 0.2 Gyr), supporting that the most massive galaxies were formed first in time, in agreement with the downsizing scenario (Cowie et al., 1996). Recent cosmological simulations also predict that the most massive galaxies form most of the stars before and quench earlier with time (Zolotov et al., 2015). We derive that the peak of the SFR occurred at t  3.7, 3.6, 3.4 Gyr with median SFR90, 120, 210 M yr and median 80, 200, 250 Myr for the low, intermediate and high mass sub-samples, respectively.

We compare these results with the scenario proposed in Thomas et al. (2005, 2010) for the formation of massive galaxies based on data in the nearby Universe (), which established that the most massive galaxies (log(M/M 12) formed 2 Gyr after the Big Bang in relatively short formation timescales (  200 Myr), while lower mass galaxies assembled later and had longer star formation episodes (   1000 Myr). The galaxies in our sample are formed  2 Gyr earlier than the results from Thomas et al. (2010) suggest. For example, the age of the Universe at the formation epoch of galaxies with log(M/M 11 proposed by Thomas et al. (2010) is 5 Gyr (  1.2), while we find 3 Gyr (  2). However, a direct comparison with Thomas et al. (2010) is not straight forward. We recall that we are limiting our study to =1.0-1.5, which restricts the latest formation epoch to t  6.0 Gyr, while Thomas et al. (2010) sample considers local galaxies which could have been formed at . The mass-dependence of the star formation timescale found by Thomas et al. (2010) is not seen in our sample. In contrast, the more massive galaxies in our sample have typically longer than the less massive galaxies. Our uncertainties are quite large ( 20%) and the mass range analysed in this work is limited to 1.5 orders of magnitude, which may be reducing the mass-dependence of the star formation timescales. But, more importantly, the median SFHs that we plot in Figure 16 are strongly affected by our selection criteria as we are only considering quiescent galaxies at 1.0    1.5. This biases our sample against less massive galaxies with long star formation episodes (long ) which were formed at   1.2 but are still forming stars at the redshift studied in this work. These galaxies did not have time to become quiescent at the redshift studied and therefore are not included in our sample. This biases the average towards lower values. The number density of quiescent galaxies with log(M/M)=10.0 – 10.5 is 10 times larger at than at , meaning that this is an important selection effect.

Age bins [Gyr] [Myr] [M yr]
  1.0 Gyr 4.1 60 230
 = 1.0 - 2.0 Gyr 3.4 200 110
  2.0 Gyr 1.9 400 60
Mass bins
log(M/M) = 10.0 - 10.5 3.6 80 90
log(M/M) = 10.5 - 10.8 3.4 200 120
log(M/M) = 10.8 - 11.5 3.1 250 210
Table 7: Median SFH of MQGs at 1.0    1.5 divided in 3   bins and 3 mass bins. We show the age of the Universe when the galaxies were formed (t), the star formation timescale () and the maximal SFR (SFR) for each subsample.

5 Summary and conclusions

In this work, we have analysed the SFH of 104 MQGs (log(M/M 10.0) at . They were selected on the basis of their location in the quiescent region of the diagram (and no IR detection) or imposing a cut in the specific star formation rate, sSFR  0.2 Gyr. We constructed the best possible SEDs by combining SHARDS spectro-photometric data, HST/WFC3 G102 and G141 grisms, and multi wavelength ancillary data form the Rainbow database. The SEDs were compared to delayed exponentially declining stellar population models and a Monte Carlo algorithm was used to characterize in detail the uncertainties and the inherent degeneracies in this kind of study.

Our main conclusions are:

  • The combined use of SHARDS and WFC3/HST grism data represents a significant improvement in the study of MQGs at =1.0-1.5. These data allowed to characterize and even break degeneracies in some cases. For  50% of our sample, we obtained a single solution in our analysis of the stellar population properties.

  • Two spectral features, the  and D4000 indices, were found to be very useful to disen