Possible identification of massive and evolved galaxies at z ≳ 5

Abstract

1 Introduction

When and how did massive quiescent or evolved galaxies like present-day elliptical galaxies form? According to the current Λ cold dark matter paradigm, dark matter halos grow by assembling material from their neighborhood, or by merging with other halos. This view is supported by cosmological simulations (White & Rees 1978; Davis et al. 1985; Springel et al. 2005; Mo et al. 2010). On the other hand, the growth of galaxies, the baryonic components of the halos, seems not to proceed in a simple bottom-up manner. Observations have shown that many of the most massive galaxies (M_star > 10¹¹ M_⊙) in the local universe assembled quickly and at early epochs (z > 2) while less massive galaxies grew more gradually (Kodama et al. 1998; Marchesini et al. 2009, 2010; Ilbert et al. 2010; Brammer et al. 2011; Muzzin et al. 2013). It is clear that the baryonic physics of galaxy formation involves more than just cooling and accretion of gas at a rate dictated by gravity.

In order to understand the balance among several processes in galaxy formation [major/minor mergers, harassment, feedback by active galactic nuclei (AGNs) or supernovae, cold gas accretion, and so on], it is important to unravel the relationship among different galaxy populations such as actively star forming galaxies, quiescent galaxies, and AGNs. The first quiescent galaxies to emerge are especially interesting because their redshifts and masses will have much to say about the history of galaxy formation and star formation activity in the very early universe. However, our present knowledge of galaxies at these epochs are based primarily on star-forming galaxies and AGNs, due to the limited sensitivity and volume of surveys carried out to-date, as well as an intrinsic rarity of quiescent galaxies at high redshift (Kajisawa et al. 2011; Domínguez Sánchez et al. 2011; Muzzin et al. 2013). A new, larger sample of quiescent galaxies at higher redshift is needed to obtain a more complete picture of galaxy formation and evolution at early epochs.

In contrast with star-forming galaxies, which typically exhibit numerous emission lines or strong Lyman breaks, quiescent galaxies are relatively difficult to identify at high redshift because of their nearly featureless spectral energy distributions (SEDs). The Balmer/4000-Å break is the only strong feature in the SEDs of quiescent galaxies, and it has long been used to identify them at z = 1–3 (Elston et al. 1988; Dickinson 1995; Yamada et al. 1997, 2005; Franx et al. 2003; Hatch et al. 2011). For example, objects with R − K > 5 or I − K > 4 in the Vega magnitude system select elliptical galaxies with a strong 4000-Å break at z ∼ 1 (Extremely Red Objects: EROs; Thompson et al. 1999; McCarthy et al. 2001; Cimatti et al. 2002; Caputi et al. 2004). This ERO selection was extended to higher redshift, z ∼ 2.4, by adopting the red z − [3.6] color criterion [f_ν(3.6 μm)/f_ν(z) > 20, IRAC-selected extremely red object (IERO): Yan et al. 2004]. The near-infrared (NIR) color criterion J − K > 2.3 (Vega) is often used to select quiescent galaxies at z = 2–3 (Distant Red Galaxies: DRGs; van Dokkum et al. 2003; Kajisawa et al. 2006a; Kubo et al. 2013). While dusty star-forming galaxies also satisfy these simple R − K and J − K color criteria, a substantial fraction of EROs and DRGs are spectroscopically confirmed to be passively evolved galaxies (Cimatti et al. 2002, 2004; Doherty et al. 2005; Kriek et al. 2006a, 2006b). As spectroscopy is time-consuming, many authors have instead separated passive and dusty star-forming galaxies photometrically using two colors so that one color brackets the Balmer break and the other color represents the shape of spectrum at longer wavelengths. The JHK color selection [J − K > 2(H − K) + 0.5 and J − K > 1.5 in the Vega magnitude system] is similar to the DRG selection, but more sensitive to the evolved galaxy population at z ∼ 2.5 (Kajisawa et al. 2006b; Doherty et al. 2010). It is also known that passive and dust-obscured galaxies (DOGs) at z > 2 separate well in the I − K vs. K − [4.5] two-color diagram (Labbé et al. 2005; Papovich et al. 2006; Kubo et al. 2013). Recently, rest-frame U − V and V − J (rest-UVJ) color selection is favored, which allows us to distinguish passive galaxies from star-forming galaxies independently of their redshifts as long as SED fitting works well (Wuyts et al. 2007; Williams et al. 2009; Whitaker et al. 2011; Brammer et al. 2011; Morishita et al. 2014).

At z > 3, the Balmer/4000 Å break is shifted to λ > 1.6 μm, and detecting quiescent galaxies at those redshifts is a challenging task. However, some candidates at z = 4–6 have been found in extremely deep NIR and IR surveys (Rodighiero et al. 2007; Wiklind et al. 2008; Mancini et al. 2009; Huang et al. 2011; Caputi et al. 2012; Nayyeri et al. 2014). Huang et al. (2011) analyzed very deep Hubble Space Telescope (HST) and Spitzer imaging in the Great Observatories Origins Deep Survey South (GOODS-South) field and identified four galaxies with H − [3.6] > 4.5. Huang et al. (2011) argued that the four EROs are likely distant passive galaxies rather than nearby dusty galaxies, because their Multiband Imaging Photometer (MIPS) 24 μm fluxes are weaker than expected if their rest-infrared SEDs are similar to those of the central region of M 82 (a local dusty galaxy with 5 < A_V < 51) or those of the DOGs at z ∼ 2. Employing a similar strategy, Caputi et al. (2012) identified 25 objects with H − [4.5] > 4 over ∼180 arcmin² of the HST Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS: Grogin et al. 2011; Koekemoer et al. 2011) Ultra Deep Survey (UDS) field. Their sample includes one z ∼ 6 massive galaxy candidate whose SED is fitted well with a 0.5-Gyr-old stellar population template with a stellar mass of 3 × 10¹¹ M_⊙. Many other red H − [4.5] objects also prefer high-redshift photometric solutions (3 < z < 5), but their red SEDs are consistent with significant dust extinction: 1.1 ≤ A_V ≤ 4.2 (Caputi et al. 2012). Wiklind et al. (2008) proposed a color selection that is more sensitive to galaxies with strong Balmer breaks (Balmer Break Galaxies: BBGs) at z = 4.9–6.5, and which is based on the three colors, K − [3.6], J − K, and H − [3.6]. They adopted their selection to the multi-band imaging available in the GOODS-South field, and identified 11 candidates of passive galaxies with ages 0.2–1 Gyr at z ≳ 5. They mentioned that the number of candidates decreases to four if they allow only objects undetected in the MIPS 24 μm.

In this work, we applied a similar color selection of high-redshift evolved galaxies to the 3.6 μm and 4.5 μm imaging from the Spitzer Extended Deep Survey (SEDS: Ashby et al. 2013) and the K-band imaging from the United Kingdom Infrared Telescope (UKIRT) Infrared Deep Sky Survey (UKIDSS: Lawrence et al. 2007). The K − [3.6] color is expected to be more sensitive to the Balmer break at z = 5–6 than H − [3.6] or H − [4.5] because it spans a relatively short wavelength interval. The wide area of the UKIDSS UDS (O. Almaini et al. in preparation ) enables us to search for rare galaxies over the entire SEDS UDS field (0.34 deg²), whereas the survey area of the previous studies using HST are relatively small. The aim of this work is to identify evolved or old stellar population galaxies at z ≳ 5.

This paper is structured as follows. In section 2 we summarize data used in this work. We describe the K − [3.6] versus [3.6] − [4.5] two-color selection scheme and further selection using multi-wavelength imaging in section 3. Stacking analysis and results for the most likely evolved galaxy candidates are described in section 4. In section 5, we discuss the implications about the high-z evolved galaxies. Our conclusions are summarized in section 6. Unless specified otherwise, we use the AB magnitude system (Oke & Gunn 1983) and adopt a cosmology with H₀ = 70.4 km s⁻¹ Mpc⁻¹, Ω_M = 0.272, and Ω_Λ = 0.728 (Komatsu et al. 2011).

2 Data

This research is primarily based on data from SEDS (Ashby et al. 2013), a very deep infrared survey using Spitzer's Infrared Array Camera (IRAC: Fazio et al. 2004) to cover five well-known extragalactic science fields. SEDS reaches a depth of 25.4 mag (0.25 μJy; 5σ) in both of the warm IRAC bands at 3.6 and 4.5 μm. In this work we analyze the ∼0.34 deg² of the SEDS UDS field, where rich multi-wavelength imaging is available from a diverse array of ground- and space-based facilities. The surveys used in this study are summarized in table 1 and figure 1.

Fig. 1.

Coverages of the survey data used in this study: SEDS, UKIDSS, SXDS, New z′-band, and SpUDS images. HerMES covers all other FoVs. We investigated objects in the overlapped region of SEDS and UKIDSS (gray shaded region). (Color online.)

Open in new tabDownload slide

The visible-wavelength photometry comes from the Subaru/XMM-Newton Deep Survey (SXDS) which covers the entire SEDS UDS area, reaching 5σ limits (2″ diameter aperture) of B = 27.5, V = 27.1, R = 26.9, and i′ = 26.8 (Furusawa et al. 2008). For z′-band imaging, we used the new data obtained in 2008–2009 (P.I.: H. Furusawa; Furusawa et al. (2016)) with Suprime-Cam (Miyazaki et al. 2002) on the Subaru Telescope. This new z′-band imaging reaches the 5σ limits of z′ = 26.3 (2″ diameter aperture), which is ∼0.5 mag deeper than that of the SXDS z′-band imaging.

The depth and coverage of the near-infrared imaging are the most important factors in this study. Deep K- or H-band imaging is required to detect or constrain the flux at shorter wavelengths than the Balmer break of quiescent galaxies at z ≳ 5. We obtained the deep J-, H-, and K-band imaging from the UKIDSS UDS (DR10¹), which have ∼0.8 deg² coverage containing nearly all of the SEDS UDS coverage. The UKIDSS project is defined in Lawrence etal. 2007. Further details on the UDS can be found in O. Almaini et al. in preparation. UKIDSS uses the UKIRT Wide Field Camera (WFCAM; Casali et al. 2007). The photometric system is described in Hewett et al. 2006, and the calibration is described in Hodgkin et al. 2009. The pipeline processing and science archive are described in Irwin et al. (in prep) and Hambly et al. 2008. In principle, wide area coverage should increase the detection rate for quiescent galaxies at high redshift because their spatial distribution is expected to be highly biased, which is the confirmed trend at 1 < z < 3 in the previous studies of EROs or DRGs (Daddi et al. 2000, 2003; Roche et al. 2002; Hamana et al. 2006; Tinker et al. 2010; Uchimoto et al. 2012). The UDS reaches depths of J = 24.9, H = 24.3, and K = 24.5 (5σ). Compared with the CANDELS/J- and H-band imaging, the UDS K-band imaging yields the advantage in the wide area and small wavelength interval between K and 3.6 μm. We estimated the above 5σ limiting magnitudes of the UDS by measuring the standard deviation of the background fluctuations with a 2″ diameter aperture. Our estimated limiting magnitudes are ∼0.3–0.4 mag shallower than the official values reported by the WFCAM Science Archive.¹ For the present work, we adopt our own measurements as the conservative limiting magnitudes.

While no survey achieved the same depth as SEDS in the mid-infrared range, shallower IRAC/5.8, 8.0, and MIPS/24 μm imaging from the Spitzer UKIDSS UDS (SpUDS; P.I.: J. Dunlop) may be useful for constraining the SEDs of galaxies. SpUDS is shallow, but covers a wide area (∼1 deg²), which contains the entire SEDS UDS coverage. We obtained the SpUDS imaging data through NASA/IPAC Infrared Science Archive (IRSA), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. We found that the limiting magnitudes are 21.7, 21.5, and 19.2 mag (7.6, 9.1, and 76 μJy; 5σ) in the 5.8, 8.0, and 24 μm bands by measuring the standard deviations of the background fluctuations, where we adopted 2 × Full Width at Half Maximum (FWHM) of the Point Spread Function (PSF) diameter apertures (4|${^{\prime\prime}_{.}}$|2, 4|${^{\prime\prime}_{.}}$|6, and 12″).

We also used SPIRE 250, 350, and 500 μm imaging from the Herschel Multi-tiered Extragalactic Survey (HerMES: Oliver et al. 2012). In the second HerMES data release (DR2), SPIRE imaging in the L4-UDS field (Oliver et al. 2012) is available, which contains the entire SEDS UDS coverage. The HerMES data was accessed through the Herschel Database in Marseille (HeDaM)² operated by CeSAM and hosted by the Laboratoire d'Astrophysique de Marseille. The depth of the SPIRE imaging is estimated by using the corresponding error maps (Levenson et al. 2010), finding that the limiting magnitudes are 14.7, 14.2, and 17.1 mJy (5σ) in the 250, 350, and 500 μm bands. We used the HerMES band-merged source catalog (xID catalog: Wang et al. 2014) for the 250, 350, and 500 μm photometry.

3 Candidate selection

3.1 Multi-band photometry for SEDS-detected objects

We used the SEDS catalog (Ashby et al. 2013) for source extraction and photometry in the 3.6 and 4.5 μm. The SEDS catalog yields PSF-fitted total magnitudes, from which we estimated 2|${^{\prime\prime}_{.}}$|4 aperture magnitudes using the aperture correction factors of Ashby et al. (2013). This small aperture size compared with the PSF (FWHM = 1|${^{\prime\prime}_{.}}$|8) is recommended by the SEDS survey team. We also investigated signal to noise ratio (S/N) in the 3.6 μm photometry as a function of aperture size, the results of which were that S/N decreases with increasing aperture size and the 2|${^{\prime\prime}_{.}}$|4 aperture is better. In this study, we focus on relatively isolated sources among the cataloged objects, because source confusion makes it hard to compare or combine multi-band photometry. In order to select isolated sources, we performed source extraction again on the 3.6 μm image by using SExtractor (Bertin & Arnouts 1996) version 2.5.0. It is suggested that SExtractor does not resolve heavy source blending at faint flux level in IRAC images (Sanders et al. 2007; Ashby et al. 2013). The software was configured carefully so that we could extract as faint isolated sources as possible. We show the SExtractor parameter settings in table 2. The objects detected in both the SEDS catalog and our own SExtractor catalog are selected for our analysis. The objects used for our analysis are shown in figure 2, where all objects in the SEDS catalog are also shown for comparison. We failed to detect not only confused but also very faint sources, which does not change the result of this paper because such faint objects in 3.6 μm cannot be identified as red K − [3.6] objects, which we focus on in the following sections. The completeness of our source extraction in the 3.6 μm image relative to the SEDS catalog is shown in figure 3. By multiplying this relative completeness and the completeness of the SEDS catalog (Ashby et al. 2013), we estimated the completeness of this work (figure 3).

Fig. 2.

Postage-stamp image of the SEDS 3.6 μm, where the panel size is 50″ × 40″. Sources used in our analysis are shown by red squares having sides of 2″. Objects in the SEDS catalog (Ashby et al. 2013) are also shown by green circles. (Color online.)

Open in new tabDownload slide

Fig. 3.

Completeness of our source extraction in the 3.6 μm image relative to that of the SEDS catalog (red circles). We also plotted the completeness that was estimated by multiplying that of the SEDS catalog (Ashby et al. 2013) and the relative completeness (black circles). (Color online.)

Open in new tabDownload slide

We measured the 2|${^{\prime\prime}_{.}}$|4 aperture magnitudes using the photometry task “PHOT” of Image Reduction and Analysis Facility (IRAF) for the B, V, R, i′, z′, J, H, K, 5.8 μm, and 8.0 μm images, where each aperture was centered on the position of the object detected in the 3.6 μm image. The B-, V-, R-, i′-, z′-, J-, H-, and K-band images were convolved to the same resolution as that of the 3.6 μm and 4.5 μm images (FWHM ∼ 1|${^{\prime\prime}_{.}}$|8). We corrected the B, V, R, i′, z′, J, H, and K magnitudes for the Galactic extinction with A_B = 0.091, A_V = 0.07, A_R = 0.056, |$A_{i^{\prime }} = 0.044$|⁠, |$A_{z^{\prime }} = 0.031$|⁠, A_J = 0.019, A_H = 0.012, and A_K = 0.0075. These values were estimated for the position (RA, Dec: 02^h18^m00^s, −05°00′00′′ in J2000) in the UDS field based on the work by Schlegel, Finkbeiner, and Davis (1998), assuming R_V = A_V/E(B − V) = 3.1. When we measured the 2|${^{\prime\prime}_{.}}$|4 diameter aperture magnitudes in the 5.8 and 8.0 μm images, we corrected the difference of PSFs to be matched with those of the 3.6 and 4.5 μm images assuming the same surface brightness profile in all bands (“negative aperture correction”). For the purpose, we smoothed the 4.5 μm image to the same resolution as that of the 5.8 and 8.0 μm images (FWHM ∼ 2|${^{\prime\prime}_{.}}$|1 and 2|${^{\prime\prime}_{.}}$|3), respectively. The correction factors were then obtained from the flux ratio of the original and smoothed images for each individual object. The average correction factor of the 5.8 and 8.0 μm images are −0.09 and −0.15 mag, respectively.

We checked whether our objects are sufficiently bright to be detected in the SpUDS 24 μm and HerMES imaging. For the SpUDS 24 μm image, object detection and photometry with 12″ diameter apertures was performed using SExtractor. For the HerMES images, the 250, 350, and 500 μm fluxes were obtained for the objects in the HerMES DR1 catalog (Oliver et al. 2012; Wang et al. 2014).

3.2 Color selection of BBGs

Fig. 4.

Behaviors of various spectral templates in the K − [3.6] versus [3.6] − [4.5] two-color diagram. The templates come from the BC03 stellar population synthesis models with the exponentially declining star-formation history. The star-formation timescale is set to τ = 0.01, 0.1, and 1 Gyr, and age is set to 0.1, 0.5, and 1 Gyr. In each panel, we show the models with different A_V (blue: A_V = 0; green: A_V = 2; magenta: A_V = 4) and different metallicity (solid: Z = 0.02; dashed: Z = 0.0001). The open squares superposed on each model track correspond to z = 1, 3, 5, and 7. The horizontal and diagonal black lines correspond to our BBG selection boundaries, where the diagonal boundary line is parallel to the reddening vector at z = 5.3 (black arrow). In the extreme right-hand side panels, we show the model template spectra with different star-formation timescale τ and age at z = 6, where the A_V and metallicity are fixed (A_V = 0 and Z = 0.02). (Color online.)

Open in new tabDownload slide

Fig. 5.

(Left-hand panel) All objects detected in the SEDS image (gray points) in the K − [3.6] vs. [3.6] − [4.5] two-color diagram. Black open circles are the BBG candidates selected by our color criteria (black horizontal and diagonal lines). Red filled circles are the most likely BBGs at z ≳ 5. Solid red and blue curves are the BC03 model tracks with different star-formation timescale, age, and A_V, where metallicity is fixed at Z = 0.02. (Right-hand panel) All objects (gray points), color-selected BBG candidates (black open circles), and most likely BBGs at z ≳ 5 (red filled circles) in the [3.6] versus K − [3.6] color–magnitude diagram. The green circles show the peaks of K − [3.6] color distribution as a function of 3.6 μm magnitude, and error bars correspond to the 4σ color errors estimated from the photometric uncertainties in the K and 3.6 μm images. (Color online.)

Open in new tabDownload slide

The BBG candidates have large uncertainties in the K − [3.6] color measurements because the K-band imaging is shallower than the SEDS 3.6 μm and 4.5 μm imaging. In the color–magnitude diagram of figure 5 (right), we show the peaks of K − [3.6] color distribution for all objects and the associated 4σ uncertainties as a function of 3.6 μm magnitudes. The 4σ color upper limits for the objects with the peak colors are below 1.3 and the number of objects is less than 10000 mag⁻¹ down to [3.6] = 24.5, which implies that typical foreground or background objects are less likely to satisfy our color criteria.

We also checked stellar contamination in our sample. First, red giant stars in the Milky Way (MW) were ruled out because of their bright absolute magnitudes ranging from −5.3 to −1.5 in the 3.6 μm (Wainscoat et al. 1992; Johnson 1966), which correspond to the apparent magnitudes ([3.6] < 18) brighter than our BBGs 100 kpc distance away. Second, we investigated the infrared colors of the brown dwarfs, which might have red K − [3.6] colors due to the cool temperature and methane absorption (Leggett et al. 2010), using the catalog of Dupuy and Liu (2012). We found that the T7.5 or later-type dwarf stars have K − [3.6] colors redder than 1.3 at the absolute magnitude range of [3.6] = 17.5–19.5, but their [3.6] − [4.5] colors are extremely red ([3.6] − [4.5] > 1). We conclude that the BBG candidates are not MW stars but galaxies.

3.3 Multi-wavelength constraint

For the eight BBG candidates, we combined multi-band photometry from B to 24 μm to make SEDs. We also checked the Herschel band-merged catalog (Wang et al. 2014), but none of the eight BBGs has a counterpart down to 8.8 mJy at 250 μm (3σ).

The eight BBG candidates may still include some contaminants. Several authors suggested that the strong emission lines from young star-forming galaxies potentially affect the broad-band photometry for high-redshift galaxies (Schaerer & de Barros 2009, 2010; Ono et al. 2010; Inoue 2011; Stark et al. 2013). Especially, strong [O iii] lines from z ∼ 4.5 galaxies or Hα lines from z ∼ 6 galaxies can make their K − [3.6] colors red. Also, some unforeseen effects caused by completely blending foreground objects might change the colors so that the objects satisfy the BBG color criteria. It is important to identify the most likely candidates of evolved galaxies using the multi-photometry data. For the eight color-selected BBGs, we required them to be undetected in the B, V, R, i′, z′, and 24 μm bands (>2σ). The observed B, V, and R bands correspond to shorter wavelengths than the Lyman break in the rest frame for z ≳ 5 galaxies. The i′ and z′ non-detection requirement was adopted to reject the strongest nebular line emitters that have strong UV continuum or Lyα line in those bands. The 24 μm non-detection requirement was adopted to reject local ultra-dusty galaxies (A_V > 10) and Type-2 AGNs. After excluding objects detected in the B, V, R, i′, z′, and 24 μm bands, we obtained only three BBGs as the most likely evolved galaxy candidates. The observed properties of these three most likely BBGs are shown in table 3.

4 Stacking analysis

4.1 Image stacking

We made the median stacked images in the B, V, R, i′, z′, J, H, K, 3.6 μm, 4.5 μm, 5.8 μm, 8.0 μm, and 24 μm for the three most likely BBGs to reduce the background noise fluctuations and constrain their average SEDs. The stacking of only three objects may cause accidental superposition of very faint foreground objects in the photometric aperture. In order to check such superposition effect, we made the 800 different stacked images for each band by rotating the individual images randomly. For these 800 stacked images we performed photometry to evaluate the average magnitude and standard deviation, and finally selected one which yields the average magnitude. The standard deviations in magnitude, which are caused by image rotations, are combined with those of background noise fluctuations to estimate the photometric errors. The stacked images typically become ∼0.5–1.0 mag deeper than the individual frames.

The stacked images are shown in figure 6. We have no detection in the B, V, R, i′, J, H, 5.8 μm, and 8.0 μm, but these non-detections yield a strong constraint in the SED fitting below. The marginal detection in the stacked z′-band image is consistent with the picture that they are galaxies with a relatively faint UV continuum and weak Lyman breaks at z ≳ 5, rather than actively star-forming galaxies with a strong UV continuum at the same redshift.

Fig. 6.

Stacked mosaics of the three most likely BBGs. Each panel is a 6″ × 6″ square on a side for the B–8.0 μm images, and 18″ for the 24 μm image. (Color online.)

Open in new tabDownload slide

4.2 SED fitting

For the stacked images, we re-performed the photometry in the same manner as that for individual objects (subsection 3.1). After obtaining the multi-band aperture photometry, we performed SED fitting in the wavelength range between 0.4 μm and 8.0 μm to study the averaged properties of the BBGs at z ≳ 5 (redshift, SFH, age, stellar mass, and so on). We used the public SED fitting code Hyperz (Bolzonella et al. 2000). By setting flux = 0 and Δ(flux) = 2σ for the flux fainter than 2σ, we kept the number of the measurements larger than that of the fitting parameters and avoided the obvious overfitting.

Spectral templates consist of the BC03 models with exponentially declining star formation rate (SFR; τ = 0.01, 0.1, 1, and 10 Gyr) and the three different metallicity values (Z = 0.0001, 0.004, and 0.02), which are the same as the templates used in figure 4. For them, redshift is changed from 0 to 7 and dust extinction A_V is changed from 0 to 4. Model ages are forced to be less than the universe age at a given redshift. We also used spectral templates with strong nebular emission lines (NL templates). In the NL templates, the BC03 model templates with constant SFR as the stellar continuum component. Our calculation of the nebular emission line flux is the same as that in Salmon et al. (2015). The number of ionizing photons was estimated from the SFH, age, and metallicity of each given model, from which we obtained Hβ line flux assuming the case B recombination and an escape fraction of ionizing photons f_esc = 0. The nebular emission line flux is calculated from their Hβ line flux and the line ratio table of Inoue (2011), who used Cloudy 08.00 (Ferland et al. 1998) to derive nebular spectra as a function of metallicity. We attenuated both nebular line flux and stellar continuum with the same dust extinction, A_V, ranging from 0 to 2. Absorption of ionizing photons (LyC absorption) by dust in the ionized nebulae was neglected in this study. Let us take the cases where nebular emission lines contribute broad-band flux at the maximum (f_esc = 0 and no LyC absorption) for the aim of this study. Our settings for the SED fitting are summarized in table 4. We adopted the Calzetti dust extinction law and the Chabrier IMF for all templates.

The measured SEDs with the acceptable-fitting templates and the fitting probability distribution which Hyperz software calculated are shown in figure 7. The best-fitting template is the old stellar population model at z = 5.7 with star-formation timescale τ = 0.1 Gyr, age ∼ 0.7 Gyr, SFR = 0.9 M_⊙ yr⁻¹, no dust, and no nebular line. However, the dusty model at z = 2.8 with A_V = 3.5 is also acceptable. We further checked this dusty solution using the 24 μm flux which is not used in the SED fitting. Dusty galaxies with A_V = 3.5 should have the significant dust emissions at longer wavelength and may be detectable in the 24 μm image. We estimated the infrared luminosity (L_IR) of the dusty template at z = 2.8 assuming that the luminosity attenuated by dust at 0.01 μm < λ < 2.2 μm is re-radiated at 5 μm < λ < 1000 μm (L_IR = L_abs), which results in L_IR = 2.3 × 10¹² L_⊙. We redshifted the model SED of a Luminous Infrared Galaxy (LIRG: Rieke et al. 2009) with the same L_IR to predict the 24 μm magnitude of the dusty template, and obtained [24] = 20.5 mag. This is inconsistent with the measured flux, [24] > 21.1 (2σ), which suggests that the dusty solution is unlikely.

Fig. 7.

Stacked SED of the three most likely BBGs (black circles) and the acceptable-fitting templates from the BC03 or NL spectral templates. Black arrows mean that the flux points are the 2σ upper limit. In the far right-hand panel, we also show the probability distribution as a function of redshift. The probability distribution for only the NL templates are shown by the dashed line. (Color online.)

Open in new tabDownload slide

We show the probability distribution for only the NL templates in figure 7 (dashed line in the right-hand panel). The strong nebular line emitter templates at z = 4.8 and 6.5 are also acceptable. While the strong nebular line emitter model cannot be rejected, we do not consider that all of the BBGs are nebular line emitters. The nebular line emitter solutions for the stacked SED have high stellar mass (M_* > 10¹⁰ M_⊙) and high dust extinction (A_V ∼ 2). Extremely strong nebular line emitters with M_* > 10¹⁰ M_⊙ and A_V ∼ 2 are very rare at z ∼ 2 (van der Wel et al. 2011; Maseda et al. 2014) and negligible especially in the photometrically selected passive galaxy sample at z ∼ 4 (Straatman et al. 2014).

From the above, we conclude that our BBG sample is expected to contain evolved galaxies at z ≳ 5. Their physical properties are expected from the best-fitting template for the stacked SED: age = 0.7 ± 0.4 Gyr, M_* = (7.5 ± 1.5) × 10¹⁰ M_⊙, SFR = 0.9 ± 0.2 M_⊙ yr⁻¹, and |$A_V = 0^{+1}_{-0}$| at z = 5.7 ± 0.6, where the uncertainties correspond to the 1σ confidence interval (Δχ² < 1) for each quantity. They are massive, but are consistent with the stellar mass of the analogs identified at intermediate redshifts (z = 2–4: Muzzin et al. 2013; Spitler et al. 2014).

5 Discussion

In the following discussion about the physical properties of the evolved galaxies in the early universe, we roughly assume that all of the three BBG candidates have a stellar masses of M_* > 5 × 10¹⁰ M_⊙ and ages older than 0.3 Gyr at the redshift 5 < z < 7. This assumption is well supported by the stacking analysis, while we cannot rule out the possibility that some of them are massive, dusty, but are extremely strong nebular line emitters (subsection 4.2).

5.1 Implication about the stellar mass density of evolved galaxies at z > 5

We estimated the number and stellar mass density of evolved galaxies at z ≳ 5. The survey volume was calculated from the SEDS area (0.34 deg²) and the redshift range (z = 5–7). The survey incompleteness was also corrected by using the completeness of our source extraction at the 3.6 μm magnitude range of the BBGs ([3.6] = 23.5–24.0), 56 ± 1% (see subsection 3.1 and figure 3). The estimated number and stellar mass density of evolved galaxies at z = 5–7 are (8 ± 5) × 10⁻⁷ Mpc⁻³ and (6 ± 4) × 10⁴ M_⊙ Mpc⁻³, respectively. These densities may be underestimated, because we might have missed some evolved galaxies as a result of photometric uncertainties. We found that about 30% of objects with an intrinsic color of K − [3.6] > 1.3 are not observed as BBGs after considering the random photometric errors which are extracted from Gaussian probability distributions. This effect, however, does not change the following discussions very much.

We compared this stellar mass density of evolved galaxies at z = 5–7 with those at z < 4 of Muzzin et al. (2013) and Straatman et al. (2014). Muzzin et al. (2013) used the COSMOS/UltraVISTA survey data to select both star-forming and quiescent galaxies up to z = 4 and to construct their stellar mass functions. Their classification of quiescent and star-forming galaxies using rest U − V and V − J colors is qualitatively similar to our BBG selection, in the sense that both selections aim to isolate galaxies with red color by Balmer breaks and blue color in longer wavelengths than Balmer breaks. We estimated the stellar mass densities at z = 0–4 by integrating their stellar mass functions down to log (M_*/M_⊙) = 10.7, which equals the stellar mass lower limit assumed for our BBG sample. Straatman et al. (2014) used the FourStar Galaxy Evolution Survey (ZFOURGE) data to construct a quiescent galaxy sample down to log (M_*/M_⊙) = 10.6 at z ∼ 4, where their selection method was the same as that in Muzzin et al. (2013) (rest-UVJ color selection). For the 15 quiescent galaxies at z ∼ 4, we adopted our stellar mass limit to calculate the stellar mass density. We show these stellar mass densities of quiescent galaxies with log (M_*/M_⊙) ≥ 10.7 in figure 8. The stellar mass density at z ∼ 6 revealed in this paper is much lower than that at z = 3–4 and is consistent with the trend that the stellar mass density of quiescent galaxies decreases continuously with increasing redshift.

Fig. 8.

(Top panel) Evolution of the stellar mass density in the universe between z = 0–7. Our measurements for evolved galaxies at z ∼ 6 are shown by red filled circles. Other symbols correspond to the measurements from previous studies: red squares for the quiescent galaxies at z = 0–4 (Muzzin et al. 2013), a red cross for the quiescent galaxies at z ∼ 4 (Straatman et al. 2014), a magenta triangle for the evolved galaxies at z ∼ 5.2 (Wiklind et al. 2008), blue squares for the star-forming galaxies at z = 0–4 (Muzzin et al. 2013), and cyan pentagons for the LBGs at z = 4–6 (González et al. 2011). All of them are corrected so that only galaxies with log (M_*/M_⊙) ≥ 10.7 contribute them. (Bottom panel) Quiescent fraction in stellar mass densities. Symbols are the same as those in the top panel. The sum of quiescent and star-forming galaxies is used as the total stellar mass density, where we used the LBGs of González et al. (2011) as star-forming galaxies for the quiescent galaxies of Straatman et al. (2014), Wiklind et al. (2008), and this work. (Color online.)

Open in new tabDownload slide

We also compared the stellar mass density of evolved galaxies at z ≳ 5 with those of Wiklind et al. (2008). The color selection of Wiklind et al. (2008) for BBGs at z ≳ 5 is similar to our BBG selection in which case they require red K − [3.6] colors, but different from our because they remove the lower redshift contaminations by using the blueward colors (J − K and H − [3.6]) instead of [3.6] − [4.5]. They identified four candidates for evolved galaxies at z = 5–7 without the 24 μm detection in the 0.04 deg² area of the GOODS-South field which have the physical properties similar to ours: log M_* = 10.7–11.7 and age = 0.2–1 Gyr. The stellar mass density after correcting the detection incompleteness is 1.4 × 10⁶ M_⊙ Mpc⁻³, which is much more than our measurement (figure 8). We cannot give a conclusive explanation for this inconsistency, although it is partially due to the fact that they did not use the nebular line emitter templates.

In figure 8, the stellar mass densities of star-forming galaxies are also plotted; they were estimated by integrating the stellar mass functions of star-forming galaxies at z = 0–4 (Muzzin et al. 2013) and Lyman Break Galaxies (LBGs) at z = 4–6 (González et al. 2011) down to log (M_*/M_⊙) = 10.7. From these, we estimated the quiescent fraction, which is defined as the fraction of the stellar mass density of quiescent galaxies compared with the total stellar mass density. Here we assume the sum of stellar mass density of quiescent galaxies and star-forming galaxies (LBG at z = 4–7) as the total stellar mass density. The quiescent fraction at z ∼ 6 is only 0.01 ± 0.0009, which is much lower than ∼0.15 at z ∼ 4. This suggests that evolved galaxies do not dominate the massive end of stellar mass functions at z > 5. Such rare but non-zero evolved galaxies are naturally considered to be embedded in prominent massive halos, because both theoretical works (e.g., De Lucia et al. 2006) and observations (Tanaka et al. 2010; Kubo et al. 2013) suggested that mass assembly and quenching of galaxies preferentially proceeded in the proto-cluster regions. It is worthwhile to investigate the surrounding environments of the evolved galaxy candidates individually in future.

5.2 Possible progenitor of evolved galaxies at z > 5

Finally, we discuss the progenitor of the evolved galaxies at z ≳ 5. Our finding of evolved galaxies at z > 5 strongly suggests intense starburst in these systems at much larger redshift. Submilimeter-selected galaxies (SMGs) are considered as possible progenitors, because their dust-obscured star formation is much more active than that of UV-selected star-forming galaxies which leads to rapid consumption of gas and subsequent quenching (Toft et al. 2014).

Although very few SMGs with spectroscopic redshifts are available above z = 5, the object HFLS3 (Riechers et al. 2013) may be a good example of a progenitor of the evolved galaxies at z ≳ 5. HFLS3 was identified by ultra-red colors in the Herschel SPIRE wavelength range (λ = 250–500 μm), and its redshift was determined to be z = 6.3369 ± 0.0009 based on several emission lines (Riechers et al. 2013). They found that HFL3 has an extremely high SFR (2900 M_⊙ yr⁻¹), a large amount of gas (1.0 × 10¹¹ M_⊙), and large stellar mass (3.7 × 10¹⁰ M_⊙). These quantities are not changed very much even if the lensing magnification caused by the two foreground galaxies is taken into account (Cooray et al. 2014). Although the gas depletion timescale of HFLS3 is 36 Myr, star formation may be suppressed on a shorter timescale (∼1–10 Myr) because of supernova (and AGN) feedback. If HFLS3 continues to produce stars with SFR = 2900 M_⊙ yr⁻¹ for 10 Myr and experiences subsequent passive evolution, it should become BBG with M_* = 6.6 × 10¹⁰ M_⊙ and age ∼0.3 Gyr at z ∼ 5, which is consistent with stellar mass and age of our evolved galaxy candidates within the 1σ uncertainties. This suggests that the most active starburst (or “maximum starburst”) SMGs can be progenitors of evolved galaxies at z ≳ 5.

We also checked their abundances. The expected number density of progenitors depends on their duty cycle. We here define duty cycle as the ratio of the time duration for which galaxies can be detected to the time interval that corresponds to the searched redshift range. Once galaxies are quenched, they can be observed as evolved galaxies at all times, which means that the duty cycle of evolved galaxies is 100%. If a progenitor has a duty cycle of 100%, the expected number density should be same as the observed number density of evolved galaxies, 8 × 10⁻⁷ Mpc⁻³. But the duty cycle of a star-forming progenitor may be small, because the time duration for which it was in active star-forming phase is expected to be short. The number density of such a progenitor should be 8 × 10⁻⁷ Mpc⁻³ per duty cycle. For “maximum starburst” SMGs such as HFLS3, we adopt the following simple assumptions. Intrinsic redshift distribution is flat at z = 5–7. Completeness of the Herschel observation of Riechers et al. (2013), which is a part of HerMES (Oliver et al. 2012) in the First Look Survey (FLS) field, Locman Hole, and GOODS-North fields, for galaxies analog to HFLS3 at z = 5–7 is 0.6. This completeness is consistent with that estimated in Wang et al. (2014) for objects with flux of 30 mJy at 500 μm (flux limit of Riechers et al. 2013). Observed time interval corresponding to z = 5–7 is 0.42 Gyr, and star-formation timescale of HFLS3 is 1–10 Myr. These assumptions lead to a spatial number density of 4.2 × 10⁻⁹ Mpc⁻³ and duty cycle of 0.24%–2.4% for “maximum starburst” SMGs at z = 5–7. The number densities of “maximum starburst” SMGs at z = 5–7 and expected progenitors as a function of duty cycle are shown in figure 9. The number density of “maximum starburst” SMGs is consistent with that of expected progenitors. Riechers et al. (2013) reported on the four analogs of HFLS3, which have similar red Herschel SPIRE colors. They may explain the small difference between the number densities of “maximum starburst” SMGs and expected progenitors, while there are large uncertainties caused by the simple assumptions. The picture that star-formation in the progenitor of evolved galaxies at z ≳ 5 rapidly proceeds surrounded with dust is a natural extension of the picture for lower-redshift, quiescent galaxies (Toft et al. 2014; Straatman et al. 2014).

Fig. 9.

The number densities of the evolved galaxy candidates (red filled circles) and “maximum starburst” SMGs (blue open squares) at z = 5–7 as a function of duty cycle. Duty cycles for “maximum starburst” SMGs are estimated from their star-forming timescale. The expected region for progenitor of evolved galaxies is shown by gray shading. (Color online.)

Open in new tabDownload slide

6 Summary

In this paper, we have investigated the nature of candidate evolved galaxies at z ≳ 5. We identified eight galaxies with K − [3.6] > 1.3 and K − [3.6] > 2.4 − [3.6]([3.6] − [4.5]) + 0.6 down to [3.6] ∼ 24 mag from the SEDS UDS 0.34 deg² area. Among them, we further selected the three most likely BBGs at z ≳ 5 by requiring non-detection in the B, V, R, i′, z′, and 24 μm bands. The best-fitting spectral template for the stacked SED is the old stellar population model with M_* = (7.5 ± 1.5) × 10¹⁰ M_⊙, SFR = 0.9 ± 0.2 M_⊙ yr⁻¹, and age = 0.7 ± 0.4 Gyr at z = 5.7 ± 0.6. We estimated the number and stellar mass densities of evolved galaxies at z ≳ 5 down to log (M_*/M_⊙) = 10.7 assuming all of the three most likely BBGs are really such objects, which results in (8 ± 5) × 10⁻⁷ Mpc⁻³ and (6 ± 4) × 10⁴ M_⊙ Mpc⁻³, respectively. This stellar mass density of evolved galaxies at z = 5–7 is consistent with the declining trend previously observed within a range of 0 < z < 4. A very simple abundance-matching test implies that the evolved galaxies identified in this paper seem to be the descendants of “maximum starburst” SMGs, such as HFLS3 reported in Riechers et al. (2013). Our finding of evolved galaxy candidates suggests that the first starburst and subsequent quenching process have been completed by z ∼ 6.

The BBGs are too red and too faint for us to further constrain their detailed properties against the photometric uncertainties. These galaxies will, however, be fully revealed by the next-generation infrared instruments on the future Wide-field Imaging Surveyor for High-redshift (WISH: Yamada et al. 2012) or James Webb Space Telescope (JWST: Gardner et al. 2006).

We thank Brett Salmon and Steven Finkelstein for the code they offered which calculates nebular emission line flux. We also appreciate Hisanori Furusawa for kindly offering the new Subaru z′-band imaging. This research was conducted as a part of K. M.'s phD thesis in Tohoku University. K. M. is grateful to Takashi Ichikawa, Masafumi Noguchi, and Masayuki Akiyama for valuable comments and suggestions. As this research started during K. M. stayed at the Harvard-Smithsonian Astrophysical Observatory, K. M. appreciates the SAO predoctoral program. This research is financially supported by Japan Society for the Promotion of Science through Brain Circulation Program (R2301) and Grant-in-aid for Scientific Research (13J02968 and 26400217).

References