{'fraction_non_alphanumeric': 0.045793341068101116, 'fraction_numerical': 0.03176528495587667, 'mean_word_length': 4.187106774870065, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 89, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Infrared color selection of massive galaxies at z>3

We introduce a new color-selection technique to identify high-redshift, massive galaxies that are systematically missed by Lyman-break selection. The new selection is based on the H_{160} and IRAC 4.5um bands, specifically H - [4.5]>2.25 mag. These galaxies, dubbed"HIEROs", include two major populations that can be separated with an additional J - H color. The populations are massive and dusty star-forming galaxies at z>3 (JH-blue) and extremely dusty galaxies at z<3 (JH-red). The 350 arcmin^2 of the GOODS-N and GOODS-S fields with the deepest HST/WFC3 and IRAC data contain 285 HIEROs down to [4.5]<24 mag. We focus here primarily on JH-blue (z>3) HIEROs, which have a median photometric redshift z ~4.4 and stellar massM_{*}~10^{10.6} Msun, and are much fainter in the rest-frame UV than similarly massive Lyman-break galaxies (LBGs). Their star formation rates (SFRs) reaches ~240 Msun yr^{-1} leading to a specific SFR, sSFR ~4.2 Gyr^{-1}, suggesting that the sSFRs for massive galaxies continue to grow at z>2 but at a lower growth rate than from z=0 to z=2. With a median half-light radius of 2 kpc, including ~20% as compact as quiescent galaxies at similar redshifts, JH-blue HIEROs represent perfect star-forming progenitors of the most massive (M_{*}>10^{11.2} Msun) compact quiescent galaxies at z ~ 3 and have the right number density. HIEROs make up ~60% of all galaxies with M_{*}>10^{10.5} Msun identified at z>3 from their photometric redshifts. This is five times more than LBGs with nearly no overlap between the two populations. While HIEROs make up 15-25% of the total SFR density at z ~ 4-5, they completely dominate the SFR density taking place in M_{*}>10^{10.5} Msun galaxies, and are therefore crucial to understanding the very early phase of massive galaxy formation.


INTRODUCTION
Our current understanding of the cosmic star formation history at z 3 is mostly based on studies of UVselected samples, e.g., Lyman-break galaxies (LBGs). Yet the LBG selection is known to be biased significantly against massive galaxies (M * 10 11 M ⊙ , because of both the relative faintness and redder UV slopes for massive galaxies (van Dokkum et al. 2006;Bian et al. 2013). On the other hand, studies of the stellar mass function based on photometric-redshift-selected galaxies from CANDELS reveals a deficiency of galaxies at the massive end at z ∼ 4 (Grazian et al. 2015), suggesting that even the deepest near-infrared-selected sample misses some massive (dusty) galaxies. Understanding the selection bias of different samples and studying star formation in more complete massive-galaxy populations at z 3 are key to mapping the full cosmic star formation history as well as understanding the very early phases of massive galaxy formation.
One of the most efficient ways of selecting massive galaxies is color selection, a method that relies on the strong spectral breaks prevalent in massive galaxies. This technique allows a rather clean selection of galaxies at certain redshifts and enables fair comparisons between different studies. A number of color-selection methods have been proposed to select z > 3 galaxies, making use of either the Balmer break (Nayyeri et al. 2014) or the 1.6 µm stellar bump (Mancini et al. 2009) as redshift indicators. One concern of these selection methods is contamination from dusty galaxies and AGNs at lower redshifts (z ∼ 2 − 3): extremely dusty galaxies or AGNs at lower redshifts present similarly red colors (due to attenuation or intrinsic red SEDs ) which mimic those of the massive galaxies at high redshift (due to Balmer or the 1.6 µm break). This is similar to what we have learned from color-selection techniques at z ∼ 2, e.g., Distant Red Galaxies (Franx et al. 2003), IRAC-selected Extremely Red Objects (Yan et al. 2004;Wang et al. 2012), and K s and IRAC-selected Extremely Red Objects . How to separate true high-redshift (z 3) from these low-redshift contaminants remains a challenge for color-selection techniques. In practice, to recover the bulk populations with color-selection techniques, we need a large range of NIR to mid-infrared (MIR) colors. This requires deep NIR and MIR imaging, which has now become available thanks to the recently completed Cosmic Assembly Nearinfrared Deep Extragalactic Legacy Survey (CANDELS; Koekemoer et al. 2011;Grogin et al. 2011), the Spitzer Extended Deep Survey (SEDS; Ashby et al. 2013) and the S-CANDELS survey (Ashby et al. 2015).
A number of recent studies attempt to provide mass-selected samples of galaxies at z > 3 based on photometric redshifts (see, e.g., Fontana et al. 2006;Pérez-González et al. 2008;Spitler et al. 2014;Straatman et al. 2014;Marchesini et al. 2014;Pannella et al. 2015;Schreiber et al. 2015). However, the quality of photometric redshifts of galaxies at z 3 is quite uncertain due to the lack of large training samples with spectroscopic redshifts. Although photometric redshifts are believed to be more reliable in deep fields, e.g., in the HUDF and GOODS fields, massive galaxies are rare, and larger fields are needed to obtain statistically significant results. In this case, color selection that requires only a few bands with high-quality photometry has the advantage of efficiently identifying large samples of certain populations of galaxies. Most existing studies are based on near-infrared-selected samples, mostly H band, which probes the rest-frame UV. As a consequence, they may miss some of the most massive and/or dustiest galaxies, as illustrated by Huang et al. (2011) (see also, e.g., Caputi et al. 2012), who reported a sample of galaxies that are bright in the IRAC bands yet undetected in deep HST /WFC3 H-band, i.e., H-dropouts. Photometric redshifts indicate that these galaxies are likely z 4 massive and dusty galaxies. Similarly, Wang et al. (2012) identified a population of massive galaxies that are bright in IRAC but not detected in the K s -band. These are believed to be massive and dusty galaxies at high redshift. Determining the number density and star formation properties of these galaxies is essential to obtain a complete view of both galaxy stellar mass functions and cosmic star formation rate densities at high redshifts. This paper presents a new color selection of Extremely Red Objects with H and IRAC colors (HIEROs, H − [4.5] > 2.25). This color-selection technique is designed to select massive galaxies at that are systematically missed by the Lyman-break selection technique. We further show that combining J − H colors breaks the degeneracy between redshift and attenuation and enables a clean selection of z > 3 galaxies. We use the deep HST /WFC3 H-band imaging from CANDELS and IRAC imaging from the SEDS survey and select both Hdetected and H-dropout HIEROs. Utilizing the exquisite multi-wavelength dataset in the GOODS fields, we explore both their star formation and structural properties. This paper is organized as follows. We describe the data and the selection of HIEROs in Section 2. Photometric redshift and stellar population analysis are presented in Sections 3 and 4, respectively. Section 5 explores star formation properties and Section 6 structural properties of HIEROs. Section 7 discusses the complete- ness of the HIERO criteria in selecting z > 3 massive galaxies, and Section 8 summarizes. Throughout the paper, we assume cosmological parameters of H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7. All magnitudes are in the AB system, where an AB magnitude is defined as AB ≡ −2.5 log(flux density in µJy) + 23.9.
2. MULTI-WAVELENGTH DATA SETS AND SAMPLE SELECTION 2.1. A combined F160W and IRAC 4.5µm-selected catalog The GOODS-South and GOODS-North fields have been the target of some of the deepest surveys ever conducted over a broad wavelength range by space observatories and the foremost ground-based telescopes. In particular, the new HST/WFC3 near-infrared survey from CANDELS and the Spitzer/IRAC mid-infrared survey from SEDS have significantly improved measurements of galaxy properties at z > 2. Here we utilize the UV to mid-infrared multi-wavelength catalogs based on detections in the HST /WFC3 F160W band from CANDELS for both GOODS-South and GOODS-North as described by Guo et al. (2013) and Barro et al. (in preparation), respectively. Both fields include deep photometry in the F435W, F606W, F775W, F814W, F850LP, F105W, F125W, F140W, and F160W bands from HST and 3.6, 4.5, 5.8, and 8.0 µm bands from Spitzer/IRAC (Ashby et al. 2013). The GOODS-South catalog also includes photometry in the U -band from both CTIO/MOSAIC and VLT/VIMOS and Ks -band imaging from the Infrared Spectrometer and Array Camera (ISAAC) and the High Acuity Wide field K-band Imager (HAWK-I) on the VLT (HUGS survey, Fontana et al. 2014). Similarly, the GOODS-North catalog includes photometry in the U -band from both KPNO and LBT and Ks-band imaging from both the Multi-Object Infrared Camera and Spectrograph (MOIRCS) and CFHT. Both fields reach 5σ depth of H ∼ 27.2 mag and are 80% complete down to H ∼ 26.2.
We also performed a systematic search for objects that are bright in the IRAC bands but are missed in the Hselected catalog, i.e., H-dropouts. We crossmatched the CANDELS H-selected catalog with an IRAC 3.6 and 4.5 µm selected catalog (Ashby et al. 2013) from the SEDS survey. The SEDS survey covers the two GOODS fields to a depth of 26 AB mag (3σ) at both 3.6 and 4.5 µm and is 80% complete down to [4.5] ∼24 mag. We first matched sources with [4.5] < 24 mag in the SEDS catalog to the H-selected catalog and identified those without H-band counterparts within a 2 ′′ radius. This 4.5 µm magnitude cut was applied to enable sufficient color range to identify extremely red objects and also give a complete 4.5 µm-selected sample. We then visually inspected the IRAC images and excluded sources whose flux is likely contaminated by bright neighbors as well as those falling on the edge of the F160W image. We dubbed this catalog of IRAC sources with no Hband counterparts the "H-dropout" catalog. Figure 1 shows examples of H-dropouts identified in GOODS-South. With knowledge of their positions, some of these H-dropouts are marginally detected in the H band but exhibit extended profiles and are unidentifiable as real sources without that prior knowledge. We measured aperture magnitudes in H and K s -bands with a 1 ′′ radius aperture at the position of their IRAC 3.6 and 4.5 µm counterparts, and then applied an aperture correction to get the total flux. For IRAC 5.8 and 8.0 µm bands, we used a 1.2 ′′ radius aperture plus aperture correction, the same as that used for IRAC 3.6 and 4.5 µm bands (directly taken from Ashby et al. (2013)). Their measured flux densities across HST /WFC3 to IRAC bands are listed in Table 1.
The combined F160W and IRAC 4.5 µm-selected catalog is not only complete to [4.5] = 24 mag but also ensures that all the H-dropouts have H − [4.5] > 2.25. As shown by Guo et al. (2013), the agreement of the IRAC  (Ashby et al. 2013) b have X-ray or 24 µm counterparts within 2 ′′ after crossmatching with X-ray and 24 µm catalogs. These galaxies are most likely AGNs if they are indeed at z > 3. But also we note that their photometric redshifts estimate with FAST is less reliable.   4.5] colors. Evolutionary tracks of a set of theoretical galaxy SED templates between z = 5 and z = 2 are shown, including an instantaneous burst (SSP) model formed at z = 10 and a constant star formation model (CSF) of age 300 Myr with different levels of reddening. The solid horizontal line shows the HIERO selection criterion adopted.The diagonal dashed line separates pure z > 3 dusty galaxies from passive galaxies at z > 3 and extremely dusty galaxies at lower redshifts. Open circles denote galaxies with spectroscopic redshifts z > 3 in the two GOODS fields. These are mostly UV-bright galaxies with lower levels of attenuation, i.e., LBGs. Right panel: The distribution of HIEROs with detections (>5σ ) in both J and H selected in the GOODS fields in the H − [4.5] versus J − H color-color diagram, color-coded by their redshifts. The diagonal dashed line separates JH-blue and JH-red HIEROs as given by Equations 1 and 2 (the same dashed line as shown in the left panel). Galaxies detected at 24 µm (F 24 µm > 30 µJy) are shown by cyan squares. Note that 24 µm-detected sources are prevalently star-forming galaxies at z < 3 and classified as JH-red HIEROs, as expected. Galaxies not detected in the F125W (J) band are shown with their 3σ upper limits. 4.5 µm photometry between the CANDELS F160Wselected catalog and the SEDS 3.6 µm and 4.5 µm selected catalog is excellent for objects with [4.5] < 24.5 mag. The magnitude cut of our selection, [4.5] = 24, is much brighter than the detection limit in the SEDS survey, and therefore Eddington bias is also negligible (Guo et al. 2013).
We searched for infrared and X-ray counterparts within a 2 ′′ radius for all the sources in both the H-selected and H-dropout catalogs based on their H-band or IRAC positions. For infrared counterparts, we employed the MIPS 24 µm-selected catalog of Magnelli et al. (2013), which also includes 100 µm and 160 µm photometry from the combination of PACS Evolutionary Probe (PEP; Lutz et al. 2011) and GOODS-Herschel (Elbaz et al. 2011) key programs. For X-ray counterparts, we used the 4 Ms catalog (Xue et al. 2011) for GOODS-South and the 2 Ms catalog (Alexander et al. 2003) for GOODS-North.
2.2. Selection of z > 3 massive galaxies At z > 3, the Balmer/4000Å break shifts redward of the H band while the 4.5 µm band probes the restframe J band. Thus both quiescent galaxies with strong Balmer/4000Å breaks and dusty galaxies with significant UV attenuation appear red in H − [4.5]. Figure 2 plots the evolution of H − [4.5] as a function of redshift for different sets of templates. These templates are based on BC03 models, including a non-evolving constant star formation (CSF) model computed for an age of 300 Myr and various levels of reddening (using the Calzetti et al. (2000) extinction law and solar metallicity). This figure illustrates that an H − [4.5] > 2.25 color cut can effectively select old or dusty galaxies at z 3. On the other hand, most commonly used LBG selection techniques (e.g., Bouwens et al. 2012) are designed specifically to select young and less attenuated galaxies, with UV slope β 0 or equivalently E(B − V ) 0.4 − 0.5 for a young star-forming galaxy. Therefore, the proposed red galaxy selection is complementary to the LBG selection and is crucial for a complete census of galaxy populations at z 3. As a further illustration of this point from Figure 2, almost none of the spectroscopically confirmed z > 3 galaxies (mostly LBGs) have H − [4.5] > 2.25. In the following sections, we refer to galaxies with H − [4.5] > 2.25 as HIEROs. Figure 2 also reveals that while passive or dusty galaxies at z > 3 are expected to be identified as HIEROs, extremely dust-obscured galaxies (E(B − V ) 0.6) at 2 < z < 3 could also enter the HIERO selection. This is similar to other red galaxy selection technique at lower redshifts. For instance, Wuyts et al. (2009) found that 15% of their distant red galaxy (DRGs, Franx et al. 2003) sample, which is intended to select galaxies at z > 2, have spectroscopic redshifts z < 2 (also see, e.g., Grazian et al. 2007). These low-redshift DRGs are on average more obscured with A V higher by 1.2 mag than the high-redshift DRGs. The situation is expected to be more serious in selecting galaxies at z > 3 due to the prevalence of dusty galaxies at 2 < z < 3 (see, e.g., Yan et al. 2007;Dey et al. 2008;Huang et al. 2009).
To enable a cleaner selection of z > 3 galaxies, we use an additional J − H color to separate z > 3 galaxies from low-redshift contaminants, following the color track of theoretical templates. As shown in Figure 3, although a heavily attenuated galaxy (E(B − V ) 0.4) at z < 3 could have a similar H − [4.5] color as normal massive , and IRAC 4.5 µm filters used to define the criteria for selecting z > 3 galaxies. Solid lines show theoretical galaxy templates for a z = 3.5 star-forming galaxy, a z = 2.5 extremely dusty galaxy, and a passive/old galaxy with an age of 1 Gyr at z = 3.5. While the three galaxy templates all present similar red H − [4.5] colors, an additional J − H color can distinguish between them due to the differences in the position and strength of the 4000Å break as well as the UV slope.
star-forming galaxies at z > 3, the Balmer/4000Å break falling between the J and H bands at z ∼ 2-3 leads to a much redder J − H (and also J − K s ) color than for galaxies at z > 3. Similarly, passive galaxies at z > 3 also have redder J − H colors due to much redder rest UV slopes than star-forming galaxies. Therefore, an additional J − H color criterion separates these different populations and approaches a pure selection of z > 3 massive (star-forming) galaxies. Based on the color tracks of theoretical models and the photometric redshifts of HI-EROs (details of photometric redshift determinations are discussed in Section 3), we separate z 3.5 star-forming galaxies (JH-blue HIEROs) from z ∼ 2−3 contaminants (JH-red HIEROs) as: Because both H − [4.5] and J − H are poorly constrained for H-dropouts, we classify all the H-dropouts as JH-blue HIEROs (the redder H − [4.5] color of Hdropouts suggests that they are in general at higher redshift). In total, we identify 359 HIEROs (116 JH-blue and 243 JH-red) in the two GOODS fields, 18 of which are H-dropouts. After examining the reliability of the color measurements of independent sources with an approach described in Section 3.2, our final sample includes 285 (80% of the original sample) HIEROs. We list respective fractions of the two categories of HIEROs in Table 2. Caputi et al. (2012) studied a sample of extremely red galaxies with H − [4.5] > 4 in UDS. Only 15 galaxies in our sample present such extremely red colors. While nearly all of these 15 galaxies have z phot > 3, the ma-   jority of the massive z > 3 galaxies will be missed by this extreme criterion. (A more extreme color cut leads to fewer contaminants at lower redshifts but also to a lower completeness in selecting high-redshift galaxies.) Similarly, Wang et al. (2012) studied a sample of K-and IRAC-selected extremely red objects (KIEROs, K s − [4.5] > 1.6) in GOODS-North, aiming to identify specifically dusty galaxies at z > 2. They showed that the majority of KIEROs are at z ∼ 2-3.5. 46% of our HIERO sample satisfy the color criterion of KIEROs. Compared to both previous studies, the advantage of our color selection is a more complete sample of massive (including both passive and star-forming) galaxies at z > 3 because our selection was specifically designed to complement the LBG selection. This allows us to perform a complete census of massive galaxy evolution at z 3. Moreover, with the proposed H − [4.5] and J − H diagram, we are able to distinguish high-redshift galaxies from low-redshift contaminants, enabling a much cleaner (and also complete) selection of massive galaxies compared to previous studies.

REDSHIFTS OF HIEROS
3.1. Photometric Redshifts The HIEROs are extremely faint at observed UV and visible wavelengths. Even though the GOODS fields have been extensively covered by spectroscopic observations, only 11 HIEROs have spectroscopic redshifts (Kajisawa et al. 2010;Dahlen et al. 2013;Hsu et al. 2014, and references therein). Therefore we use photometric redshifts to gain insight into their nature. The unique features of these galaxies, e.g., many of them are consistent with being extremely dusty with A V exceeding 3 − 4, lead to concerns that they may be not represented in the templates used by most photometric redshift meth- Figure 6. Examples of SED fitting when including/excluding IRAC photometry. For each galaxy, from left to right the inner panels show the F160W and IRAC 4.5 µm images with the image size of 4 ′′ × 4 ′′ . When including the IRAC photometry, a bad fit was achieved with χ 2 4. This suggests that most likely the IRAC photometry for these galaxies is contaminated by neighboring sources, as can also been seen from the stamp images.
ods. Hence we used galaxy templates spanning a larger parameter space particularly a larger range of A V to determine their photometric redshifts. To this aim, we employed FAST (Kriek et al. 2009) to fit the full U -band to 8.0 µm photometry for all galaxies in the HIERO sample. FAST also provides self-consistent estimates of stellar masses. We constructed stellar templates from the Bruzual & Charlot (2003) (BC03, hereafter) stellar population synthesis model with a Chabrier (2003) initial mass function and solar metallicity, assuming exponentially declining star-formation histories (SFHs) with efolding times τ = 0.1-10 Gyr. We allowed the galaxies to be attenuated with 0 ≤ A V ≤ 6 mag with reddening following the Calzetti et al. (2000) law. To avoid strong influence on the fitting from one single band (in some cases due to an emission line or bad photometry), we restricted the maximum S/N in the photometry to be 20. Examples of the fitting are shown in Figure 4.
For the few galaxies with spectroscopic redshifts in our sample, Figure 5 shows the comparison between spectroscopic and photometric redshifts. Our photometric redshift estimates in general agree with spectroscopic redshifts with the normalized median absolute deviation (Brammer et al. 2008) σ NMAD ∼ 0.06. (The significant outlier at z spec = 1.2 is a type-1 AGN.) However, the spectroscopic sample is significantly biased with 9 of 11 objects detected in X-rays and a median H-band magnitude H ≈ 22.3 (compared to H ≈ 23.9 for the total sample). Hence, larger and deeper spectroscopic samples are needed to verify the accuracy of the photometric redshifts. On the other hand, based on photometric redshifts estimates from CANDELS (see e.g., Dahlen et al. 2013), we derive σ NMAD ∼ 0.11. Specifically for HIEROs at z 3, the photometric redshifts from CANDELS are systematically higher by ∼0.2-0.3, likely due to the absence of highly attenuated galaxy templates in the CAN-DELS photometric redshift codes, (favoring a high redshift solution to account for the red color actually due to attenuation). This trend is also confirmed based on the photometric comparisons between ours and from CAN-DELS for the whole HIERO sample.

Identifying sources with unreliable IRAC
photometry While the majority of HIEROs yield a good fit, a substantial number of galaxies can not be well fitted with FAST and the BC03 models. We checked in detail their SEDs and F160W and IRAC images and found that most of them have bright stars or galaxies within a few arcseconds, hence their photometry is likely unreliable. This is particularly a problem in the IRAC bands due to their larger PSFs. To explore whether this is the origin of the problem, we re-fit the SEDs with the same set of templates but excluding the IRAC photometry. A better fit (much smaller χ 2 ) was achieved for many sources with a significantly different redshift solution. Figure 6 presents two such examples. This illustrates that indeed the IRAC photometry for these sources is problematic (in most cases boosted) due to contamination by close neighbors. As a result, the true H − [4.5] color of these sources is likely much bluer (as can also be seen from their best-fitted SED in the second fit).
Although in principle we can simply reject all sources with a bad fit in the first run to clean our sample, it will introduce a bias towards the choice of SED-fitting methods and templates, i.e., some sources with a bad fit may be due to the wrong templates instead of bad photometry. To be conservative, we borrow the "clean index" concept in dealing with the source confusion at far-infrared wavelengths (Elbaz et al. 2011) to identify sources whose IRAC fluxes are reliable. We define a source as "clean" only if it satisfies at least one of the following two criteria: (a) it has no neighbor within 2 ′′ (rougly the size of the FWHM of the IRAC PSF) and a good SED(photometric redshift)-fitting result (χ 2 < 4, corresponding roughly to the one-tailed (right-tail) probability > 0.95) when including IRAC photometry; (b) similar redshift solutions are achieved during the two SED-fitting runs, i.e.,  |z phot (no IRAC) − z phot (IRAC)|/z phot (IRAC) < 30% (For sources that only satisfy the second criterion, it is unclear whether their IRAC photometry is contaminated or not: the bad fit could be because there are no perfect templates in the library or because the data quality is not good enough.) This leaves us 285 clean HIEROs including the 18 H-dropouts (Table 2). Among these 285 HIEROs, 223 of them (78%) have a good fit with χ 2 < 4 when including the IRAC photometry. Figure 7 shows the photometric redshift distribution for the final clean samples. Our classification based on J − H and H − [4.5] colors successfully separates galaxies at z 3.5 from those at relatively lower redshifts, i.e., 2 < z < 3. The median redshifts for JH-blue and JH-red HIEROs are z ∼ 4.4 and z ∼ 2.5, respectively. Among the 79 JH-blue HIEROs, only 14 (18%) have photometric redshifts z phot < 3. Most of these 14 are detected at MIPS 24 µm, in contrast to the z > 3 JH-blue HIEROs. On the other hand, 52 out of 206 (28%) of JH-red HIEROs are at z > 3 with the majority (30 out of 52) being at 3 < z < 3.5. Figure   HIEROs in our sample. Solid lines show the 2 < z < 5 color tracks for a passive galaxy with an age of 1.0 Gyr and a starforming galaxy with an age of 0.3 Gyr and constant star formation history. All galaxies are color-coded with their best-fit specific star formation rates (sSFR). The dashed line denotes the dividing line that separates quiescent from star-forming galaxies. The arrow shows the effect of 1 mag of dust extinction at z = 3 assuming a Calzetti law.

STELLAR POPULATIONS OF HIEROS
4.1. JH-blue HIEROs: Normal massive star-forming galaxies at z > 3 Our analysis on photometric redshifts (both from our own estimates based on FAST and those from CAN-DELS) confirms that the selection criterion of JH-blue HIEROs yields a clean selection of z > 3 galaxies. Remaining low-redshift contaminants can be removed based on SED modeling or on their strong 24 µm detections (with a median F 24 µm ∼ 98 µJy). The F160W images of these z < 2 contaminants show that most of them are either extremely faint or likely mergers, leading to uncertain J − H colors or unusual extinction properties, which explains why they enter our selection.
The advantage of using FAST to derive photometric redshifts is that it also gives self-consistent estimates of stellar masses. Figure 8 shows the stellar mass estimates of JH-blue HIEROs. The modeling assumptions are described in Section 3, i.e., the BC03 stellar library, Calzetti et al. (2000) extinction law with 0 ≤ A V ≤ 6, and exponentially declining star formation histories (single τ model) with e-folding times ranging from 0.1 to 10 Gyr. Delayed τ models give nearly identical estimates of stellar masses. Our mass estimates imply that the JH-blue HIEROs are massive star-forming galaxies at z 3.5 with a median M * ≈ 10 10.6 M ⊙ .
The best-fit ages for JH-blue HIEROs from single τ models range from t = 0.1 to t = 1.6 Gyr with a median t ≈ 1 Gyr for galaxies at z > 3 and t ≈ 0.1 Gyr for galaxies a z ≤ 3. The median attenuation is E(B − V ) = 0.33, with galaxies at z < 3 much dustier than those at z > 3 ( E(B − V ) = 0.83 vs. E(B − V ) = 0.25). These best-fit stellar properties do not change significantly by using delayed τ models. These results suggest that most JH-blue HIEROs are massive, dusty star-forming galaxies which have already assembled relatively old stellar populations. Their red H − [4.5] colors thus appear to be caused by a combination of moderately old stellar populations (strong 4000Å break) and dust attenuation. The few z < 3 contaminants tend to be less massive, younger, and also much dustier, suggesting that the red H − [4.5] colors are mostly caused by severe dust attenuation.
4.2. JH-red HIEROs: Massive dusty star-forming galaxies at 2 < z < 3 and passive galaxies at 3 < z < 4 As discussed in Section 2.2 and Section 3, JH-red HI-EROs include in general two populations, z ∼ 2-3 dusty star-forming galaxies and passive galaxies at z ∼ 3-4. As illustrated in Figure 3, it is possible to distinguish between the two populations based on their different behaviors in the observed near-infrared to mid-infrared colors. Here we propose to separate the two populations with the observed J − K versus K − [4.5] color-color diagram, which for z ∼ 3 resembles the rest-frame U − V versus V − J diagram (Williams et al. 2009). Figure 9 shows the distribution of JH-red HIEROs in this diagram. Galaxies with lower sSFR (from SED fitting) populate the upper left region while galaxies with higher sSFR locate in the lower right region. As independent evidence, most 24 µm detections fall in the star-forming region, consistent with star-forming galaxies at z < 3. Based on the distribution of galaxies with different sSFRs and the direction of reddening from the Calzetti et al. (2000) law, we define the criteria for quiescent galaxies as : J − K > (K − [4.5]) + 1.2, J − K > 2.2, and K − [4.5] < 1.6. The remaining galaxies are classified as star-forming. This diagram is not only limited to classifications of HIEROs but also provides an efficient way to identify quiescent galaxies at z 3 in general.
As shown in Figure 9, roughly 10% (21 out of 206) of the JH-red HIEROs are classified as quiescent galaxies. The median redshift of these is z ≈ 3.4, suggesting that they are among the earliest quenched systems in the Universe. Among these 21 quiescent JH-red HIEROs, 15 are at z phot > 3 out of 52 z > 3 JH-red HIEROs altogether.
The JH-red HIEROs are also primarily massive galaxies (the right panel of Figure 8) with a M * ≈ 10 10.7 M ⊙ for both quiescent and star-forming subpopulations. The median best-fit age and attenuation for quiescent galaxies are t = 1.0 Gyr and E(B − V ) = 0.17 while for starforming galaxies t = 0.7 Gyr and E(B − V ) = 0.6, respectively. These results are consistent with the classifications based on the observed J − K versus K − [4.5] diagram. The high attenuation value for star-forming galaxies in this JH-red HIERO population is also consistent with their high 24 µm flux densities with F 24 µm ≈ 90 µJy, which corresponds roughly to a total infrared luminosity L TIR ∼ 10 12 L ⊙ at z ∼ 2.5.
Based on X-ray luminosity estimates, ∼20% of JH-red HIEROs are classified as X-ray AGNs (L 0.5−8 keV > 10 42 erg s −1 ), consistent with the X-ray AGN fraction in massive galaxies at z ∼ 2 (Wang et al. 2015, submitted). The JH-blue HIEROs have similar AGN fraction, ∼18%. These results indicate that our selection criteria are not particularly biased towards AGNs compared to non-AGN massive galaxies. Figure 10. Stellar masses as a function of UV luminosity for LBGs (B-dropouts, triangles) and JH-blue HIEROs (pentagons) in GOODS fields. Only H-detected HIEROs are shown to ensure that the rest-frame UV magnitudes are reliable. The best-fit M * -L UV,1600 relation for B-dropouts (González et al. 2011) and massselected samples (Grazian et al. 2015) at z ∼ 4 are shown with the dashed and solid line, respectively.

UV properties
By selection, the HIEROs are faint in their restframe UV despite their high stellar masses. Specifically, the JH-blue HIEROs are generally massive star-forming galaxies at z > 3. Previous studies of LBGs at these redshifts reveal a tight correlation between stellar mass and UV luminosity (see, e.g., González et al. 2011). This allows the UV luminosity function to translate directly to the stellar mass function, as is commonly adopted in studying high-redshift stellar mass functions. It is therefore interesting to compare HIERO and LBG stellar masses at similar redshifts and UV luminosities to explore how the HIEROs affect the stellar mass estimate.
We selected a sample of z ∼ 4 LBGs (B-dropouts) that are on average at similar redshifts as the JH-blue HIERO population in the two GOODS fields using the same criteria as Bouwens et al. (2012). The restframe UV luminosities (M 1600 ) were derived using EaZY (Brammer et al. 2008) at fixed redshifts from CANDELS. Figure 10 compares the M * -L UV relations. (The HIERO rest-frame UV luminosities came from the best-fit SED templates at fixed redshifts estimated by FAST.) Compared to LBGs at similar masses, the HI-EROs are generally 2-3 magnitudes fainter in the restframe UV. The existence of these galaxies suggests that using a simple M * -L UV relation (or using only UVselected samples) to determine the stellar mass function underestimates the massive end. We will illustrate this point further in Section 7.

Determining total infrared luminosities of HIEROs through stacking
The general faintness of HIEROs in the UV suggests that most of their star formation is hidden by dust. The amount of this hidden star formation can be inferred through infrared emission, which originates from the thermal reradiation by dust of their absorbed ultraviolet light. Therefore, understanding the Spectral Energy Distributions (SEDs) in the infrared is essential to get a complete view of star formation in HIEROs. Alternatively, we could use dust-unbiased tracers of star formation, e.g., radio continuum, to measure the total star formation rate. However, although both the far-infrared and radio surveys in the GOODS fields are among the deepest ever conducted, only a few HIEROs are individually detected. Moreover, most of the detected ones are JH-red HIEROs. For instance, crossmatching with the VLA 1.4 GHz catalog in GOODS-North (Morrison et al. 2010), 33 out of 142 HIEROs are detected at 1.4 GHz with F 1.4 GHz > 20 uJy, but only 6 out of 33 are classified as JH-blue HIEROs. This is likely caused by the high redshifts of the JH-blue HIERO population and that most of them are not extreme starbursts (contrary to bright submillimeter galaxies). To probe lower star formation rates that are typical of HIEROs, a stacking approach is required.
For the JH-blue HIERO population, we excluded contaminants at low redshifts by applying a redshift cut of z > 3, which leaves us 66 galaxies. For the JH-red HIERO population, we separately stacked star-forming (SF) and quiescent (QS) galaxies. The exquisite multiwavelength data in GOODS allow us to perform stacking across the whole infrared wavelength range, including the 16 and 24 µm bands from Spitzer, 100, 160, 250, 350 and 500 µm from Herschel, 850 µm from SCUBA, 870 µm from LABOCA, and 1.1 mm from AzTEC. This permits a comprehensive understanding of their infrared SEDs. Moreover, the combination of stacked far-infrared and submillimeter colors provides an independent and complementary estimate of their redshift from the position of their peak far-infrared emission (Hughes et al. 2002;Daddi et al. 2009).
To avoid contamination from a few relatively bright members, we conducted a median stacking in the infrared and submillimeter bands (16 µm-1.1 mm) using the IAS stacking code library (Béthermin et al. 2010). (Using mean stacking would leave the main results unchanged, likely because the fraction of bright members is small.) We retrieved flux densities by PSF-fitting the stacked images. A correction factor ranging from 0.92   . Stacked SEDs for JH-red SF (left) and JH-blue (right) HIEROs. Median flux densities between F435W and F160W bands were derived based on stacking of the corresponding HST images, while the median K-band and IRAC flux densities were derived from the measured fluxes of individual sources. The best-fitting IR SED templates from Chary & Elbaz (2001) are shown as dashed curves, while the best-fitting full-band SEDs from CIGALE are shown with solid lines. The derived total infrared luminosities and median redshifts for the two populations are also shown. at 250 µm to 0.75 at 500 µm was applied to account for clustering, which does not change much with redshift and stellar mass of the galaxies (Béthermin et al. 2012;Schreiber et al. 2015). We determine uncertainties on the flux densities with bootstrapping. The stacked flux densities for all populations are listed in Table 3. For JH-blue HIEROs, significant detections are revealed except at 16 and 70 µm due to shallower depths. For JH-red SF HIEROs, significant detections are found at all wavelengths while there are no detections with S/N > 3 in any infrared bands for JH-red QS HIEROs. This provides independent evidence that our approach successfully separates the QS and SF populations. Moreover, the peak of the infrared SED for JH-blue and JH-red SF HIEROs falls respectively at ∼500 µm and ∼350 µm, lending evidence that they are most likely at z ∼ 4 and z ∼ 2.5, consistent with our photometric redshifts. Figure 11 shows the median SEDs of the two populations of HIEROs. With the well-constrained SED shape, their total infrared luminosities (TIR) constrain their star formation rates. We fit the stacked 160 µm−1.1 mm (to avoid rest-frame <40 µm wavelengths, which may suffer from AGN contamination) SEDs using a suite of infrared templates, including the 105 template SEDs from Chary & Elbaz (2001) (CE01, hereafter). During the fitting, we fixed the redshift at the median value of the sample and left the template normalizations as free   3) and JH-red SF HIEROs are shown as red and blue filled pentagrams, respectively. The horizontal error bars represent the 1σ scatter of the redshift distribution. We also derived sSFR for each population of HIEROs in narrower redshift bins, which are shown with red and blue open pentagrams for JH-red and JH-blue HIEROs, respectively. The dashed line indicates the best-fit of the sSFR evolution at z < 2 (Sargent et al. 2014), while the solid line indicates the predicted sSFR evolution from theoretical models (Dekel & Mandelker 2014), normalized at z ∼ 1. Given that the HIEROs represent the majority of galaxies with M * > 10 10.5 M ⊙ , our results suggest that the sSFR for massive galaxies continues to increase up to z ∼ 4 but then becomes relatively flat at higher redshifts. Right panel: The evolution of star formation rate density (SFRD) as a function of redshift. The dashed and dotted lines denote the evolution of SFRD presented by Madau & Dickinson (2014) and Behroozi et al. (2013), respectively. The blue open pentagrams show the SFRD of JH-blue HIEROs at two redshift bins: 3 < z < 4.4 and 4.4 < z < 6, the same redshift bins as shown by the blue open pentagrams in the left panel. Compared to previous SFRD measurements based on LBGs at these redshifts, the HIEROs make up 15 − 25% of that by LBGs.
parameters. The best-fit model is the template that minimizes χ 2 . The total infrared luminosity for JH-red and JH-blue HIEROs are 1.2 × 10 12 L ⊙ and 2.4 × 10 12 L ⊙ , respectively. We also fit the full-band SED from UV to far-IR for HIEROs using the code CIGALE (Noll et al. 2009;Serra et al. 2011) using BC03 models in the UVto-NIR and Draine & Li (2007) models in the infrared. This gives a consistent measurement of L T IR as shown in Figure 11.
To derive realistic uncertainties on star formation rates and specific star formation rates, we bootstraped galaxies simultaneously in all the infrared bands, i.e., bootstrapping SEDs. Each time we randomly selected a subsample of the galaxies, performed median stacking in all the bands, and determined the SFR and sSFRs from the TIR and median stellar mass of the subsample. We repeated this process 50 times and determined the dispersion of SFR and sSFR. This method avoids the drawback that the uncertainties in different bands as derived from bootstrapping galaxies in a single band are likely correlated, e.g., galaxies that are fainter in the shorter wavelength may tend to be brighter in the larger wavelength due to the variations in dust temperature and/or redshifts.

Dust attenuation of HIEROs
A tight correlation between L IR /L UV and UV continuum slope β exists for UV-selected L * star-forming galaxies in the local universe and up to z ∼ 2.5 (Meurer et al. 1999;Kong et al. 2004;Reddy et al. 2012). At higher redshifts, it becomes more difficult to directly measure both L IR and β. Based on infrared stacking analysis, Lee et al. (2012) derived IR luminosities for a statistical sample of L >> L * , z ∼ 4 LBGs and showed that they are consistent with the IRX-β relation presented by Meurer et al. (1999). For HIEROs, because most of them are relatively faint in the rest-frame UV, we first performed a median stacking across the observedframe HST/ACS F435W to HST/WFC3 F160W bands for JH-blue HIEROs and JH-red HIEROs, respectively, and then measured their flux densities based on the stacked images. We then used the measurements for rest-frame 1400-2800Å to derive β. This process is illustrated in Figure 12. The strong break between observedframe B 435 -, V 606 -, and I 775 -band for JH-blue HIEROs provides additional independent evidence that most HI-EROs should be at similar redshifts as B-dropouts and V -dropouts, i.e., z ∼ 4-5.
The right panel of Figure 12 compares the relation between L IR /L UV and β for both HIEROs and LBGs. The HIEROs have significantly redder UV slopes (β ∼ 0) than the brightest/most massive LBGs (β ∼ −1.9) at similar redshifts. On the other hand, the effective dust attenuation, L IR /L UV , is extremely high for HIEROs, reaching ∼1000 for both JH-blue HIEROs and JH-red HIEROs. The JH-red HIEROs fall on the Meurer relation within uncertainties while JH-blue HIEROs tend to be above the Meurer relation, which is consistent with the results on a K s -selected massive galaxy sample at z ∼ 3.3 by Pannella et al. (2015). This is at odds with recent findings on UV-selected samples at z ∼ 5, which are systematically below the Meurer relation (Capak et al. 2015). This may be caused by UV selection tending to select less massive and less dusty galaxies, leading to a biased view of galaxy populations at high redshift.

Star formation rates of HIEROs
Using a Kennicutt (1998) Table 4. Figure 13 shows the sSFR for JH-blue and JH-red SF HIEROs. The sSFR for JH-blue HIEROs, 4.2 Gyr −1 , is twice as high as that for similarly massive galaxies at z ∼ 2 (Schreiber et al. 2015). This suggests that the sSFR of massive galaxies continues to increase at z > 2.5 and up to z ∼ 4. However, the growth rate at z 2.5 is slower than that at z ∼ 0-2. For instance, Sargent et al. (2014) found that the evolution of sSFR at z 2.5 roughly follows (1 + z) 2.8 , which predicts sSFR ∼ 7 Gyr −1 , much higher than we observed here at z ∼ 4. Based on simple analytic arguments on the accretion rates into halos and the accretion of baryons into galaxies, Dekel & Mandelker (2014) showed the evolution of sSFR for typical galaxies should follow sSF R ∼ (1 + z) 2.5 , which also predicts a higher sSFR at z ∼ 4 than observed for HIEROs. Several recent studies report similar slow evolution of sSFR at z 2.5 but for less massive or for UV-selected galaxies (e.g., González et al. 2014;Tasca et al. 2015).
Our finding suggests that the slow evolution of sSFR at z 2.5 (or the fast evolution of sSFR at z < 2.5) is likely a universal behavior for all masses of galaxies.
Most previous studies on the cosmic star formation rate densities at z 4 contain only contributions from LBGs. The cosmic SFRD contributed by HIEROs needs to be added to have a complete view of cosmic SFR densities at high redshift. Using a total area of 340 arcmin 2 and assuming that individual HIEROs have the same sSFR as derived from stacking, the total SFRD contributed by HIEROs is ∼6 × 10 −3 M ⊙ yr −1 Mpc −3 at z ∼ 3.7 and ∼8 × 10 −3 M ⊙ yr −1 Mpc −3 at z ∼ 4.7. By comparison, recent studies showed that the total SFRD of LBGs is ∼0.049 and ∼0.036 M ⊙ yr −1 Mpc −3 at z ∼ 4 and z ∼ 5, respectively . Thus despite the small number of HIEROs, they contribute 15-25% to the SFRD at z ∼ 4-5, not taken into account in previous studies based on UV-selected samples. (There is essentially no overlap between the HIERO and LBG selections.) Figure 13 plots the SFRD from HI-EROs identified in this work and from LBGs in previous work Finkelstein et al. 2015).

STRUCTURAL PROPERTIES OF HIEROS
Recent work has revealed a significant population of massive compact quiescent galaxies at z ∼ 1.5-3 with typical ages of ∼1 Gyr (Daddi et al. 2005;Buitrago et al. 2008;Szomoru et al. 2012;Cassata et al. 2013;van der Wel et al. 2014). It is suggested that the formation of these galaxies must have been through a compact star-forming phase (Barro et al. 2013;Patel et al. 2013;Stefanon et al. 2013;Williams et al. 2014;Zolotov et al. 2015). Indeed, a significant population of massive compact star-forming galaxies exists at z ∼ 2-3 (Barro et al. 2013(Barro et al. , 2014. These are likely to form the bulk population of compact quiescent galaxies at z 2. However, when and how these compact galaxies formed remains unclear. Moreover, a significant number of passive galaxies have already formed at z ∼ 3-4 (Gobat et al. 2012;Fan et al. 2013;Buitrago et al. 2013;Straatman et al. 2014), as also shown in this work, and the star-forming progenitors of these earliest quenched systems remain elusive (Straatman et al. 2015). Studying structural properties of JH-blue HIEROs, which are representative massive galaxies at z > 3, should provide us important insights into these questions.
For both the JH-blue and the JH-red QS HIERO populations, we limit sample galaxies to those with photometric redshifts z > 3. We then retrieved their F160Wband half-light radii from the CANDELS structural parameters catalogs as described by van der Wel et al. (2014). The median F160W magnitudes for the JH-blue HIERO and JH-red QS HIERO are H = 25.5 mag and H = 24.5 mag, respectively. Therefore, in the left panel of Figure 14 we only plot galaxies that have reliable size measurements (with flag = 0 in the catalog of structural parameters, as described in van der Wel et al. (2012)). This leaves us 22 (out of 66) JH-blue and 11 (out of 15) JH-red QS HIEROs with respective median F160W magnitudes of H = 25.0 mag and H = 24.5 mag. Quiescent galaxies at z ∼ 3.4 are as compact as their z ∼ 2.75 counterparts.
There are similar numbers of compact star-forming HI-EROs compared to quiescent HIEROs, though the starforming ones tend to be at higher redshifts. If we define galaxies below the 1σ upper bound of the mass-size relation for z ∼ 2.75 quiescent galaxies as compact, then there are 7/11 compact quiescent and 10/22 compact star-forming HIEROs, respectively. This corresponds to a number density of 9.3 × 10 −6 Mpc −3 and 6.7 × 10 −6 Mpc −3 for compact quiescent (3 < z < 4) and starforming HIEROs (3.5 < z < 5), respectively. The actual fraction of compact galaxies for star-forming HIEROs are likely lower because most of them are not shown in this figure due to unreliable size measurements, which tend to have larger sizes based on stacking. This seems at odds with the findings of Straatman et al. (2015), who revealed a paucity of compact star-forming galaxies in their z ∼ 4 galaxy samples. We argue that this is likely due to the difference in the sample selection methods and cosmic variance. (Straatman et al. selected their sample based on H and K s catalogs and focused on galaxies with M * > 10 10.6 M ⊙ at 3.4 < z < 4.2.) The majority of the JH-blue HIEROs are too faint to have reliable size measurements, so we measured the size of the stacked F160W image (Figure 14). We derived the half-light radius of the stacked image with single Sérsic model fitting using GALFIT (Peng et al. 2010), as shown in Figure 15. The best-fit Sérsic model yields n = 2.5 and r e = 2.0 kpc (assuming a median redshift of z ∼ 4.4) for star-forming HIEROs while it yields n = 4.3 and r e = 1.3 kpc (assuming a median redshift of z ∼ 3.4) for quiescent HIEROs. Based on the mass-size relation of quiescent galaxies at z ∼ 2.75 ( Figure. 14, van der Wel et al. 2014), the stellar mass of quiescent galaxies with similar sizes is M * ∼ 10 11.17 M ⊙ . With current SFR of 240 M ⊙ yr −1 , the HIEROs can reach this stellar mass within ∼0.6 Gyr, less than the time interval between the median redshift of HIEROs (z ∼ 4.4) and z ∼ 2.75 quiescent galaxies. This time scale, ∼ 0.6 Gyr, is similar to the gas depletion time for typical (main-sequence) star-forming galaxies at z ∼ 1 − 3 (Tacconi et al. 2013). We thus conclude that the HIEROs can evolve into the most massive quiescent galaxies at z ∼ 2.75 by in situ star formation and subsequent quenching. Moreover, the number density of massive quiescent galaxies with M * > 10 11.17 M ⊙ at z ∼ 2.75 is ∼2.0 × 10 −5 Mpc −3 (Muzzin et al. 2013;Ilbert et al. 2013), similar to the number density of JH-blue HI-EROs, 2.6 ×10 −5 Mpc −3 . The JH-blue HIEROs are in general at much higher redshift than the compact starforming galaxies identified by Barro et al. (2013) and Williams et al. (2014) and are also more massive than those identified by Williams et al. (2014). We propose that the HIEROs in our sample likely include the majority of the progenitors of the most massive, also likely to be the first quenched, quiescent galaxies at z ∼ 2.75. Further investigations of the HIEROs would be key to unveil the formation mechanism of the earliest-quenched massive (compact) quiescent galaxies.

MASSIVE GALAXIES AT Z > 3
The new color selection technique identifies massive red galaxies at z > 3. Both photometric analysis on individual optical-to-NIR SED and the stacked infrared SED suggest that we reveal a population of z > 3 massive galaxies that was largely missed in previous UV-selected samples. A critical question is then to what extent these HIEROs represent massive galaxy populations at z 3.
Unfortunately it is difficult to precisely estimate the selection efficiency of HIEROs because there are no complete spectroscopically confirmed massive galaxy samples at these redshifts. Bearing in mind the uncertainty in photometric redshifts (especially for those galaxies that are neither LBGs nor HIEROs), we seek a rough estimate of the selection efficiency of the HIERO criterion in selecting massive galaxies at z > 3 by comparing to photometric-redshift-selected samples. Figure 16 presents the stellar mass versus redshift for all galaxies with [4.5] < 24 and z > 3 in the two GOODS fields, which are selected based on the F160W-band selected catalogs from CANDELS and those H-dropouts identified in this work. For galaxies in the CANDELS catalog, we use the stellar mass and redshifts from the official CANDELS catalog while for the H-dropouts we use the stellar mass and redshifts estimated in this work. We emphasize again that no significant differences are found for JH-blue HIEROs in terms of their redshift and stellar mass estimates between this work and those from CAN-DELS. Figure 16 shows that ∼60% of the galaxies with M > 10 10.5 M ⊙ at z 3.5 can be identified by our HI-ERO selection, while only ∼10% are selected as LBGs (B 435 -, V 606 -, or I 775 -dropouts) using the same criteria as Bouwens et al. (2012). Moreover, ∼20% of the massive galaxies are H-dropouts, which are not included in the CANDELS H-band-selected catalog. To explore whether these results strongly depend on the photometric redshift estimation, we also derived the fraction of HIEROs and LBGs at the high mass end in the GOODS-South field using photometric redshift estimation from  Hsu et al. (2014), yielding similar results on the fraction of massive galaxies selected as H-dropouts and HIEROs in general. These results suggest that the HIEROs dominate the high-mass end of the stellar mass function at z > 3, and their properties are thus representative of the massive galaxies at these redshifts. On the other hand, the LBG selection misses the majority of the most massive galaxies at high redshifts and significantly underestimates the high-mass end of the stellar mass function. Similarly, using the H-band selected samples alone would also miss a significant fraction of massive galaxies. Figure 17 presents the comparison between the number density of HIEROs and the galaxy stellar mass functions (SMF) based on both UV and H-selected samples. We separated the HIEROs with M * > 10 10.5 M ⊙ into two stellar mass bins. Figure 17 reveals that the HI-EROs dominate at the high-mass end and represent the key population that reconciles the discrepancies between UV-based SMFs (Song et al. 2015) and SMFs based on stellar-mass-selected samples (Grazian et al. 2015) at z ∼ 4-5 (but we caution that the uncertainty at the massive end of the SMF for both estimates is quite large). On the other hand, we also compare particularly the number density of H-dropouts with the SMF based on H-band-selected samples (Grazian et al. 2015). In this way we can derive the fraction of massive galaxies missed by H-band selection. As shown in the bottom panel of Figure 17, at M * > 10 11 M ⊙ , the number density of H-dropouts is comparable to that of H-detected galaxies (down to the limit of the CANDELS survey). This suggests that, even based on H-band-selected samples  (Song et al. 2015) and H-band selected samples (Grazian et al. 2015) at z ∼ 4-5. The massive end of the stellar mass function at 3 < z < 4 from the COSMOS survey is also shown (Muzzin et al. 2013;Ilbert et al. 2013). The HIEROs are separated into two stellar mass bins.
(or more precisely, mass-selected samples based on Hband catalogs), we may still substantially underestimate the most massive end and get a steeper exponential tail while the true exponential tail may be much shallower. Unfortunately, spectroscopic confirmation of these most massive galaxies at the highest redshift is not possible with current optical and NIR facilities. Instead, JWST and ALMA would be the most promising tools to fully address this question.

CONCLUSION
This paper demonstrates a new H and IRAC color selection technique (HIEROs: H − [4.5] > 2.25) to identify massive, UV-faint galaxies at z > 3 that are systematically missed by the Lyman-break selection technique. The HIEROs also include a significant population of massive and dusty galaxies at z ∼ 2-3, which can be separated from the true z > 3 star-forming galaxies with an additional J − H color, enabling a clean selection of (dusty) star-forming galaxies at z > 3. The HIEROs dominate the high-mass end of the stellar mass function, making up 60% of galaxies with M * > 10 10.5 M ⊙ while LBGs contribute only 10%. The fact that only the J, H, and IRAC 4.5 µm bands are involved in this selection allows us to efficiently select large samples of massive galaxies at z > 3.
The high-redshift nature of HIEROs are independently confirmed through their stacked UV-to-NIR and FIR SEDs: the stacked rest-frame UV SED resembles those of B 435 -and V 606 -dropouts while the stacked FIR SED peaks at 500 µm. Based on the stacked SEDs, UV and infrared properties of HIEROs are representative of massive star-forming galaxies at z > 3. They are 2-3 magnitudes fainter in the rest-frame UV than LBGs with the same stellar mass and tend to be above the IRX-β relation. Thus both stellar mass and star formation rates based purely on UV are underestimated. The z > 3 HIEROs have typical SFR ∼240 M ⊙ yr −1 and sSFR ∼4.2 Gyr −1 , double the rates for similarly massive z ∼ 2 galaxies and suggesting that the sSFR for massive galaxies continue to increase at z > 2.5 yet with a decreased growth rate compared to that at z < 2.5. This is consistent with recent findings for less massive and for UVselected galaxies.
There are similar numbers of compact quiescent and star-forming galaxies among the z > 3 HIEROs, with the star-forming ones being at higher redshifts than quiescent galaxies (z ∼ 4.4 vs. z ∼ 3.4). This suggests that even the earliest quenched systems may have gone through a compact star-forming phase which started at an even earlier epoch, providing important constraints on the formation of compact quiescent galaxies. For the total z > 3 HIERO population, both their number densities and sizes match those of the most massive (M * > 10 11.2 M ⊙ ) quiescent galaxies at z ∼ 2.75, providing the most plausible star-forming progenitors.
The HIERO selection provides a reliable and representative sample of massive galaxies at z 3. Although rough constraints on their typical physical properties can be obtained through stacking, their general faintness in the rest-frame UV-to-optical inhibits accurate determination of the physical properties for individual galaxies with current optical and NIR facilities. On the other hand, their brightness in the mid-to far-infrared make JWST and ALMA promising tools to further explore their nature in great detail.