PHYSICAL PROPERTIES OF SPECTROSCOPICALLY CONFIRMED GALAXIES AT $z \ge 6$. III. STELLAR POPULATIONS FROM SED MODELING WITH SECURE Ly$\alpha$ EMISSION AND REDSHIFTS*

LINHUA JIANG¹, KRISTIAN FINLATOR²,<fnref n="11" />, SETH H. COHEN³, EICHI EGAMI⁴, ROGIER A. WINDHORST³, XIAOHUI FAN⁴, ROMEEL DAVÉ⁵, NOBUNARI KASHIKAWA⁶, MATTHEW MECHTLEY³, MASAMI OUCHI⁷,⁸, KAZUHIRO SHIMASAKU⁹, AND BENJAMIN CLÉMENT¹⁰

¹ Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, China; jiangKIAA@pku.edu.cn

²Dark Cosmology Centre, Niels Bohr Institute, University of Copenhagen, Denmark

³School of Earth and Space Exploration, Arizona State University, Tempe, AZ 85287, USA

⁴Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA

⁵Physics Department, University of the Western Cape, 7535 Bellville, Cape Town, South Africa

⁶Optical and Infrared Astronomy Division, National Astronomical Observatory, Mitaka, Tokyo 181-8588, Japan

⁷Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan

⁸Kavli Institute for the Physics and Mathematics of The Universe, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8583, Japan

⁹Department of Astronomy, University of Tokyo, Hongo, Tokyo 113-0033, Japan

¹⁰Centre de Recherche Astrophysique de Lyon, Universite Lyon 1, 9 Avenue Charles Andre, F-69561 Saint Genis Laval Cedex, France

Received 2015 August 2; accepted 2015 November 3; published 2015 December 23

ABSTRACT

We present a study of stellar populations in a sample of spectroscopically confirmed Lyman-break galaxies (LBGs) and Ly$\alpha$ emitters (LAEs) at $5.7 < z < 7$. These galaxies have deep images from *Subaru*, the *Hubble Space Telescope*, and *Spitzer*/IRAC. We focus on 27 galaxies with IRAC detections, and characterize their stellar populations based on the multi-band data and secure redshifts. By estimating nebular emission from the observed Ly$\alpha$ flux, we break the strong model degeneracy between young galaxies with prominent nebular emission and older galaxies with strong Balmer breaks. The results show that our galaxies cover a wide range of ages from several to a few hundred million years (Myr), and stellar masses from $\sim 10^8$ to $\sim 10^{11} M_\odot$. These galaxies can be roughly divided into two subsamples: an “old” subsample consisting of galaxies older than 100 Myr, with stellar masses higher than $10^9 M_\odot$, and a “young” subsample consisting of galaxies younger than $\sim 30$ Myr, with masses ranging between $\sim 10^8$ and $\sim 3 \times 10^9 M_\odot$. Both subsamples display a correlation between stellar mass and star formation rate (SFR), but with very different normalizations. The average specific SFR (sSFR) of the “old” subsample is $3-4 \text{ Gyr}^{-1}$, consistent with previous studies of “normal” star-forming galaxies at $z \ge 6$. The average sSFR of the “young” subsample is an order of magnitude higher, likely due to starburst activity. Our results also indicate little dust extinction in the majority of the galaxies, as already suggested by their steep rest-frame UV slopes. Finally, LAEs and LBGs with strong Ly$\alpha$ emission are indistinguishable in terms of age, stellar mass, and SFR.

Key words: galaxies: evolution – galaxies: high-redshift

1. INTRODUCTION

Star-forming galaxies at $z \ge 6$ are natural tools to study the early galaxy formation and explore the history of cosmic reionization. In recent years, with the advances of instrumentation on the *Hubble Space Telescope* (HST) and large ground-based telescopes, the number of galaxies found at $z \ge 6$ has increased dramatically. The majority of these galaxies were photometrically selected Lyman-break galaxies (LBGs) using the dropout technique (e.g., Oesch et al. 2010; Yan et al. 2012; Ellis et al. 2013; Willott et al. 2013; Bouwens et al. 2015; Laporte et al. 2015). Some of them, among the brightest in the optical and near-IR, have been spectroscopically confirmed with deep observations (e.g., Jiang et al. 2011; Curtis-Lake et al. 2012; Ono et al. 2012; Schenker et al. 2012; Finkelstein

et al. 2013; Pentericci et al. 2014; Matthee et al. 2015; Oesch et al. 2015; Sobral et al. 2015; Watson et al. 2015). A complementary way to find $z \ge 6$ galaxies is the narrow-band (or Ly$\alpha$) technique. This technique has a high success rate of spectroscopic confirmation. More than 200 Ly$\alpha$ emitters (LAEs) at $z \ge 6$ have been spectroscopically identified (e.g., Shimasaku et al. 2006; Hu et al. 2010; Ouchi et al. 2010; Kashikawa et al. 2011; Henry et al. 2012), including several at $z \ge 7$ (e.g., Iye et al. 2006; Rhoads et al. 2012; Shibuya et al. 2012; Konno et al. 2014). Many properties of high-redshift LAEs can be used to probe cosmic reionization (e.g., Silva et al. 2013; Treu et al. 2013; Cai et al. 2014; Dijkstra 2014; Jensen et al. 2014; Momose et al. 2014).

Meanwhile, the physical properties of $z \ge 6$ galaxies are also being investigated. At $z \ge 6$, the rest-frame UV/optical light moves to the IR range. Therefore, IR observations, including the *HST* near-IR and *Spitzer* Space Telescope (*Spitzer*) mid-IR observations, are essential to measure the properties of these galaxies. While properties such as the rest-frame UV slope and galaxy morphology can be directly measured from optical and near-IR images (e.g., McLure et al. 2011; Dunlop et al. 2013; Bouwens et al. 2014; Curtis-Lake et al. 2014; Kawamata et al. 2015), detailed physical

<footnote>* Based in part on observations made with the NASA/ESA *Hubble Space Telescope*, obtained from the data archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. Based in part on observations made with the *Spitzer Space Telescope*, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Based in part on data collected at *Subaru Telescope* and obtained from the SMOKA, which is operated by the Astronomy Data Center, National Astronomical Observatory of Japan.</footnote>

<footnote n="11">¹¹DARK Fellow.</footnote>
properties such as age and stellar mass have to come from spectral energy distribution (SED) modeling of stellar populations based on the combination of optical, near-IR, and mid-IR data (Conroy 2013). The optical data and *HST* near-IR data measure the slope of the rest-frame UV spectrum, and constrain the properties of young stellar populations. *Spitzer* IRAC provides mid-IR photometry. When combined with *HST* near-IR data, it measures the amplitude of the Balmer break and constrains the properties of older stellar populations. Soon after the launch of *Spitzer*, it was found that IRAC is sensitive enough to directly detect luminous $z \ge 6$ LBGs (e.g., Egami et al. 2005; Eyles et al. 2005; Yan et al. 2005). In these early results, $z \simeq 6$ LBGs showed strong IRAC detections, suggesting the existence of established massive stellar populations in their galaxies. Later studies then found that most galaxies were actually not detected in moderately deep IRAC images, meaning that they were considerably younger and less massive (e.g., Yan et al. 2006; Eyles et al. 2007; Pirzkal et al. 2007). Extensive studies are now being carried out with various galaxy samples, and a diversity of physical properties are found (e.g., González et al. 2010; Schaerer & de Barros 2010; McLure et al. 2011; Curtis-Lake et al. 2013; Labbé et al. 2013; de Barros et al. 2014). These studies are mostly based on photometrically selected samples, and their galaxies are not spectroscopically confirmed.

This paper is the third in a series presenting the physical properties of a large sample of spectroscopically confirmed galaxies at $z \ge 6$. In the first paper of the series (Jiang et al. 2013b, hereafter Paper I), we presented deep *HST* and *Spitzer* observations of 67 spectroscopically confirmed galaxies at $z \ge 6$. The sample is the largest collection of spectroscopically confirmed galaxies in this redshift range, including 51 LAEs and 16 LBGs. We measured basic properties of the rest-frame UV continuum and $Ly\alpha$ emission in these galaxies. In the second paper of the series (Jiang et al. 2013a, hereafter Paper II), we carried out a structural and morphological study of these galaxies. In this third paper, we will fit and model the SEDs of the galaxies and derive physical parameters such as age, stellar mass, and dust extinction. This is particularly important for high-redshift LAEs. Some early observations (e.g., Malhotra & Rhoads 2002; Pirzkal et al. 2007) found that $z \ge 4.5$ LAEs consist of young stellar populations with low stellar masses (compared to LBGs). However, this could be due to selection effects (e.g., Dayal & Ferrara 2012). For LAEs at $z \ge 5.7$, we currently have very little knowledge of their stellar populations. This is because almost all of the known LAEs were discovered by ground-based telescopes, and do not have sufficiently deep infrared observations. Our paper includes a large sample of LAEs. This allows us, for the first time, to systematically study stellar populations in $z > 5.7$ LAEs.

As mentioned in Paper I, the spectroscopic redshifts of this sample have great advantages for measuring physical properties of high-redshift galaxies. Secure redshifts remove one critical free parameter in the SED modeling. A model spectrum of a bright $z \ge 6$ galaxy is usually derived from four to five broadband photometric points, e.g., one optical band, two *HST* bands, and one to two *Spitzer* IRAC bands. Given the limited degrees of freedom, a spectroscopic redshift will significantly improve SED modeling, especially when gaseous or nebular emission is considered. It has been clear that strong nebular emission widely exists in high-redshift galaxies. For example, Stark et al. (2013) found that $H\alpha$ contributes more than 30% of

the IRAC 1 flux in $4 < z < 5$ galaxies. In the analysis of a large LBG sample at $3 < z < 6$, de Barros et al. (2014) estimated that about 60%-70% of the LBGs show prominent nebular emission lines. At $z \ge 6$, strong nebular lines such as [O III] $\lambda$5007, $H\beta$, and $H\alpha$ enter the IRAC 1 and 2 bands (3.6 and 4.5 $\mu$m), which significantly affects the measurements of stellar populations (e.g., Robertson et al. 2010; Schaerer & de Barros 2010; Finlator et al. 2011; Stark et al. 2013; de Barros et al. 2014). Photometric redshifts with large uncertainties may place these nebular lines in the wrong filter bands during SED fitting. With spectroscopic redshifts, the predicted observed wavelengths of nebular lines are securely known.

The structure of the paper is as follows. In Section 2, we briefly review our galaxy sample and the available data. In Section 3, we perform SED modeling of the galaxies using evolutionary synthesis models. We then present the derived stellar populations of the galaxies in Section 4 and discuss the results in Section 5. We summarize the paper in Section 6. Throughout the paper, we adopt a $\Lambda$-dominated flat cosmology with $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_m = 0.3$, and $\Omega_\Lambda = 0.7$. All magnitudes are on the AB system (Oke & Gunn 1983).

## 2. GALAXY SAMPLE AND DATA

Our galaxy sample consists of 67 spectroscopically confirmed galaxies at $z \ge 6$, including 62 galaxies in the *Subaru* Deep Field (SDF; Kashikawa et al. 2004) and 5 galaxies in the *Subaru* XMM-Newton Deep Survey field (SXDS; Furusawa et al. 2008). The SDF sample has 22 LAEs at $z \simeq 5.7$, 25 LAEs at $z \simeq 6.5$, one LAE at $z = 6.96$, and 14 LBGs at $5.9 \le z \le 6.5$. The SXDS sample contains 3 LAEs at $z \simeq 6.5$ and 2 LBGs at $z \simeq 6$. All of the LAEs at $z \simeq 5.7$ and 6.5 have a relatively uniform magnitude limit of 26 mag in the narrow bands NB816 and NB921, and thus make a well-defined sample. The LBGs were selected with different criteria, and have rather inhomogeneous depths, so they are not a statistically complete sample. The details of our galaxy sample can be found in Section 2 of Paper I. In the next two subsections, we will briefly describe the data used for this paper.

### 2.1. Optical and Near-IR Imaging Data

The optical imaging data for the two fields SDF and SXDS were obtained with the *Subaru* Suprime-Cam (Kashikawa et al. 2004; Furusawa et al. 2008). They consist of images in a series of broad and narrow bands. In Paper I, we produced a set of stacked images in six broad bands ($BVRI'z'y$) and three narrow bands (NB816, NB921, and NB973), by including all available data in the archive. Our stacked images have great depth with an excellent point-spread function (PSF) FWHM of $0.6''-0.7''$. In particular, the total integration times of the SDF $z'$- and y-band images are 29 hr and 24 hr, corresponding to depths of 27.1 mag and 26.2 mag, respectively ($5\sigma$ detection in a $2''$ diameter aperture). For the galaxies at $z < 6$, these two bands do not cover the $Ly\alpha$ emission line, so they provide two important photometric points for SED modeling.

We obtained near-IR imaging data for the SDF galaxies in three *HST* GO programs (11149: PI E. Egami; 12329 and 12616: PI L. Jiang). The *HST* observations were made with NICMOS and WFC3. The majority of the galaxies were observed with WFC3 in the F125W (hereafter $J_{125}$) and F160W (hereafter $H_{160}$) bands. The typical integration time was two
<table><caption>**Table 1**

Mid-IR Photometry of the Galaxies in Our Sample</caption><thead><tr><th>No.</th><th>R.A. (J2000)</th><th>decl. (J2000)</th><th>Redshift</th><th>IRAC 1</th><th>IRAC 2</th></tr></thead><tbody><tr><td>2</td><td>13:23:54.601</td><td>+27:24:12.72</td><td>5.654</td><td>&gt;26.85</td><td>...</td></tr><tr><td>3</td><td>13:24:16.468</td><td>+27:19:07.65</td><td>5.665</td><td>25.09 ± 0.21</td><td>25.39 ± 0.33</td></tr><tr><td>4</td><td>13:24:32.885</td><td>+27:30:08.82</td><td>5.671</td><td>24.78 ± 0.22</td><td>...</td></tr><tr><td>5</td><td>13:24:11.887</td><td>+27:41:31.81</td><td>5.681</td><td>&gt;26.64</td><td>...</td></tr><tr><td>10</td><td>13:24:33.097</td><td>+27:29:38.58</td><td>5.696</td><td>&gt;26.35</td><td>...</td></tr><tr><td>15</td><td>13:24:23.705</td><td>+27:33:24.82</td><td>5.710</td><td>24.04 ± 0.07</td><td>24.29 ± 0.11</td></tr><tr><td>17</td><td>13:23:44.747</td><td>+27:24:26.81</td><td>5.716</td><td>&gt;26.56</td><td>...</td></tr><tr><td>20</td><td>13:24:40.527</td><td>+27:13:57.91</td><td>5.724</td><td>25.20 ± 0.30</td><td>...</td></tr><tr><td>21</td><td>13:24:30.633</td><td>+27:29:34.61</td><td>5.738</td><td>&gt;26.16</td><td>...</td></tr><tr><td>22</td><td>13:24:41.264</td><td>+27:26:49.09</td><td>5.743</td><td>&gt;26.43</td><td>...</td></tr><tr><td>23</td><td>13:24:18.450</td><td>+27:16:32.56</td><td>5.922</td><td>25.11 ± 0.16</td><td>...</td></tr><tr><td>24</td><td>13:25:19.463</td><td>+27:18:28.51</td><td>6.002</td><td>25.15 ± 0.22</td><td>...</td></tr><tr><td>25</td><td>13:24:26.559</td><td>+27:15:59.72</td><td>6.032</td><td>24.30 ± 0.07</td><td>24.47 ± 0.13</td></tr><tr><td>27</td><td>13:24:10.766</td><td>+27:19:03.95</td><td>6.040</td><td>25.70 ± 0.30</td><td>...</td></tr><tr><td>28</td><td>13:24:42.452</td><td>+27:24:23.35</td><td>6.042</td><td>25.59 ± 0.39</td><td>...</td></tr><tr><td>29</td><td>13:24:05.895</td><td>+27:18:37.72</td><td>6.049</td><td>&gt;26.56</td><td>...</td></tr><tr><td>30</td><td>13:24:00.301</td><td>+27:32:37.95</td><td>6.062</td><td>25.50 ± 0.35</td><td>...</td></tr><tr><td>31</td><td>13:23:45.632</td><td>+27:17:00.53</td><td>6.112</td><td>25.06 ± 0.20</td><td>...</td></tr><tr><td>33</td><td>13:24:20.628</td><td>+27:16:40.47</td><td>6.269</td><td>&gt;26.62</td><td>...</td></tr><tr><td>34</td><td>13:23:45.757</td><td>+27:32:51.30</td><td>6.315</td><td>23.82 ± 0.09</td><td>...</td></tr><tr><td>35</td><td>13:24:40.643</td><td>+27:36:06.94</td><td>6.332</td><td>23.23 ± 0.05</td><td>23.74 ± 0.11</td></tr><tr><td>36</td><td>13:23:45.937</td><td>+27:25:18.06</td><td>6.482</td><td>25.00 ± 0.19</td><td>25.20 ± 0.28</td></tr><tr><td>37</td><td>13:24:18.416</td><td>+27:33:44.97</td><td>6.508</td><td>&gt;26.61</td><td>...</td></tr><tr><td>39</td><td>13:23:43.190</td><td>+27:24:52.04</td><td>6.534</td><td>&gt;26.43</td><td>...</td></tr><tr><td>40</td><td>13:24:55.772</td><td>+27:40:15.31</td><td>6.534</td><td>&gt;26.59</td><td>...</td></tr><tr><td>43</td><td>13:23:53.054</td><td>+27:16:30.75</td><td>6.542</td><td>25.30 ± 0.24</td><td>...</td></tr><tr><td>44</td><td>13:24:15.678</td><td>+27:30:57.79</td><td>6.543</td><td>23.77 ± 0.05</td><td>24.02 ± 0.09</td></tr><tr><td>45</td><td>13:24:40.239</td><td>+27:25:53.11</td><td>6.544</td><td>&gt;26.39</td><td>...</td></tr><tr><td>46</td><td>13:23:52.680</td><td>+27:16:21.76</td><td>6.545</td><td>&gt;26.38</td><td>...</td></tr><tr><td>47</td><td>13:24:10.817</td><td>+27:19:28.08</td><td>6.547</td><td>23.70 ± 0.05</td><td>23.74 ± 0.09</td></tr><tr><td>48</td><td>13:23:48.922</td><td>+27:15:30.33</td><td>6.548</td><td>&gt;26.35</td><td>...</td></tr><tr><td>49</td><td>13:24:17.909</td><td>+27:17:45.94</td><td>6.548</td><td>25.33 ± 0.24</td><td>...</td></tr><tr><td>50</td><td>13:23:44.896</td><td>+27:31:44.90</td><td>6.550</td><td>24.56 ± 0.16</td><td>...</td></tr><tr><td>52</td><td>13:24:35.005</td><td>+27:39:57.43</td><td>6.554</td><td>&gt;26.48</td><td>...</td></tr><tr><td>54</td><td>13:24:08.313</td><td>+27:15:43.49</td><td>6.556</td><td>25.14 ± 0.25</td><td>25.21 ± 0.36</td></tr><tr><td>58</td><td>13:24:43.427</td><td>+27:26:32.62</td><td>6.583</td><td>25.44 ± 0.26</td><td>24.72 ± 0.23</td></tr><tr><td>61</td><td>13:25:22.291</td><td>+27:35:19.95</td><td>6.599</td><td>24.08 ± 0.07</td><td>24.50 ± 0.18</td></tr><tr><td>62</td><td>13:23:59.766</td><td>+27:24:55.75</td><td>6.964</td><td>24.84 ± 0.16</td><td>25.17 ± 0.21</td></tr><tr><td>63</td><td>02:18:00.899</td><td>-05:11:37.69</td><td>6.023</td><td>24.49 ± 0.12</td><td>24.97 ± 0.25</td></tr><tr><td>64</td><td>02:17:35.337</td><td>-05:10:32.50</td><td>6.116</td><td>24.90 ± 0.11</td><td>...</td></tr><tr><td>66</td><td>02:18:20.701</td><td>-05:11:09.89</td><td>6.575</td><td>25.90 ± 0.32</td><td>...</td></tr><tr><td>67</td><td>02:17:57.585</td><td>-05:08:44.72</td><td>6.595</td><td>23.72 ± 0.04</td><td>24.25 ± 0.08</td></tr></tbody></table>

Note. The sequence numbers of the galaxies in Column 1 correspond to the numbers in Column 1 of Table 1 in Paper I.

HST orbits (roughly 5400 s) per band, leading to a depth of ~27.4 mag (5σ detection) in $J_{125}$ and a depth of ~27.1 mag in $H_{160}$ (see also Windhorst et al. 2011). The remaining several SDF galaxies were observed with NICMOS in the F110W (hereafter $J_{110}$) and $H_{160}$ bands. The typical integration time was also two HST orbits, and the depths in the two bands are ~26.4 and ~26.1 mag, respectively. The five SXDS galaxies were covered by the UKIDSS Ultra-Deep Survey (UDS). Their HST WFC3 near-IR data were obtained from the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011). The exposure depth of the CANDELS UDS data is 1900 s in $J_{125}$ and 3300 s in $H_{160}$, slightly shallower than our WFC3 data for the SDF. All of the optical and near-IR photometry is listed in Table 1 of Paper I.

Our galaxies represent the most luminous galaxies at $z \ge 6$, in terms of Ly$\alpha$ luminosity (for LAEs) or UV continuum

luminosity (for LBGs). They cover the brightest UV luminosity range of $M_{1500} < -19.5$ mag, so the majority of them were detected ($>5\sigma$) in the near-IR images. In Papers I and II, we have shown that these galaxies have steep UV continuum slopes $\beta$ with a median value of $\beta = -2.3$. They have moderately strong rest-frame Ly$\alpha$ equivalent widths (EWs) in the range of ~10 to ~200 Å. Their star formation rates (SFRs) are moderate from a few to a few tens of solar masses per year. These galaxies also exhibit a wide range of rest-frame UV continuum morphology in the HST images, from compact features to multiple component systems. In this paper, we will measure stellar populations in these galaxies.

## 2.2. Spitzer Mid-IR Imaging Data

Our Spitzer IRAC imaging data for the SDF were obtained from two GO programs, 40026 (PI: E. Egami) and 70094 (PI:
L. Jiang). Program 40026 was carried out during the *Spitzer* cryogenic phase, and the other one was carried out in the Warm Mission phase. The two programs imaged roughly 70% of the SDF to a depth of 3–7 hr. They covered all 62 SDF galaxies in our sample with IRAC channel 1 (3.6 µm), and 51 (out of 62) galaxies with IRAC channel 2 (4.5 µm). The IRAC data were reduced with the *Spitzer* Science Center (SSC) pipeline MOPEX. The details of the *Spitzer* observations and IRAC data reduction were described in Paper I. The final co-added images have a pixel size of 0.′6, roughly half of the IRAC native pixel scale.

The IRAC images for the SXDS galaxies were obtained from the *Spitzer* Extended Deep Survey (SEDS; Ashby et al. 2013). SEDS is a very deep imaging survey in the IRAC 1 and 2 bands over five well-studied fields, including the UDS. With an integration time of 12 hr per pointing, SEDS reaches 26 mag in IRAC 1. The SEDS co-added images also have pixel scales of 0.′6. The IRAC thumbnail images of all galaxies are shown in Figure 12 of Paper I.

The IRAC mid-IR photometry is complicated by source confusion. In deep IRAC images, faint galaxies are often blended with nearby neighbors, so reliable photometry usually requires proper deblending and removal of neighbors. We performed source deblending using iGALFIT (Ryan 2011), an interactive tool to run GALFIT (Peng et al. 2002). For each galaxy in our sample, the basic procedure is as follows. We first modeled its bright neighbors with iGALFIT. The model neighbors were convolved with a PSF image, and were subtracted from the original image. The PSF image was constructed from a number of bright (but unsaturated) point sources. We then carried out aperture photometry for this galaxy on the residual image. We used a 3 pixel (1.′8) aperture radius and computed its background in an annulus from 5 to 10 pixels. Finally, we applied an aperture correction (roughly 0.4 mag), which was measured from the PSF image. If a galaxy was isolated from any bright neighboring objects, we performed aperture photometry directly, without removing neighbors.

The results of the above IRAC photometry are shown in Table 1, where we list the magnitudes (or 2σ upper limits) and errors for 42 (out of 67) galaxies in our sample. We have discarded the galaxies that have less than two broadband photometric points in the optical and near-IR (too few data points for SED modeling). These discarded galaxies are all very faint in terms of their rest-frame UV continuum emission. They were barely (or not) detected in the HST $J_{125}$ (or $J_{110}$) band, the deepest band that we have. This virtually puts a magnitude limit on our sample: galaxies fainter than $J_{125} \sim 27.2$ mag were discarded. In Table 1, we have also discarded the galaxies that are heavily blended with (or completely covered by) much brighter neighbors in the IRAC images. In these cases, the IRAC photometry of the galaxies is not reliable. In Table 1, the sequence numbers of the galaxies in Column 1 correspond to the numbers in Column 1 of Table 1 in Paper I. Column 2 lists the redshifts. Columns 3 and 4 list the aperture photometry in IRAC 1 and 2 channels. If a galaxy is not detected in IRAC 1, a 2σ upper limit is given. We do not give upper limits for IRAC 2. Our IRAC 2 data are often significantly shallower than the IRAC 1 data, thus the inclusion of the IRAC 2 upper limits put no significant constraints on SED modeling.

### 3. SED MODELING

In this section, we perform SED modeling using the GALEV evolutionary synthesis models (Kotulla et al. 2009). The GALEV models are similar to other evolutionary synthesis models such as BC03 (Bruzual & Charlot 2003) and STARBURST99 (Leitherer et al. 1999). One distinct feature of GALEV is that it provides an option to include metallicity-dependent gaseous or nebular emission (both continuum and line emission). As we will see, nebular emission is critical for the SED modeling of our galaxies.

For high-redshift galaxies, the quality of the SED fitting is usually dominated by data quality rather than the quality of synthesis models (e.g., Pirzkal et al. 2012), i.e., it is limited by the number of available photometric data points and photometric uncertainties. The galaxies in our sample, like $z \ge 6$ galaxies in many other samples, have only three to five broadband photometric points available, including one to two ground-based optical points, two HST near-IR points, and one to two IRAC mid-IR points. In addition, these galaxies are faint, with relatively large photometric uncertainties, particularly in the IRAC bands. Among 42 galaxies shown in Table 1, 27 galaxies were detected ($>3\sigma$) in our IRAC 1 images, and 13 galaxies were also detected in the IRAC 2 images. In order to account for the associated systematic uncertainties, we model our measurements under a broad range of assumptions, with the view that the true physical parameters of our galaxies lies somewhere within the range of possibilities that we consider.

#### 3.1. Stellar Populations with and without Nebular Emission

Given the limited number of available photometric data points, we use as few free parameters as possible in our SED fitting. Compared to photometric samples, the advantage of our sample are the spectroscopic redshifts that remove one critical free parameter. The other major parameters for SED fitting are metallicity, dust extinction (or reddening), stellar mass, age, and star formation history (SFH). We adopt a Salpeter initial mass function with a mass range of 0.1–100 $M_{\odot}$. For metallicity, the GALEV models provide two options: chemically consistent treatment or fixed metallicity values. We choose to use the fixed metallicity values as other synthesis models do. There is a strong age–metallicity degeneracy at young ages (many of our galaxies are young), and our data are not sufficient enough to break this degeneracy. So we fix metallicity to be 0.2 $Z_{\odot}$, which was suggested by the simulations in, e.g., Finlator et al. (2011).

We use two representative SFHs, an exponentially declining SFH (dSFH) and a smoothly rising SFH (rSFH). For the exponentially dSFH~exp($-t/\tau$), we fix the decline factor $\tau$ to be 200 million years (Myr; to reduce the number of free parameters). The SFR of this SFH declines slowly at young ages (younger than ~200 Myr) and is not sensitive to $\tau$ for $\tau > 200$ Myr. At high redshift, a more realistic SFH is probably a smoothly rSFH (e.g., Finlator et al. 2007, 2011; Papovich et al. 2011). GALEV has not included any rSFHs yet, so we incorporate the smoothly rSFH (SFR as a function of age) of Finlator et al. (2011) into the GALEV models. Hereafter, we denote the above exponentially declining SFH and smoothly rising SFH as “dSFH” and “rSFH,” respectively. Our purpose of using two SFH models is to explore the possible ranges of physical parameters, rather than distinguishing one model from the other. We do not use more complex
SFHs that usually introduce new free parameters. We do not consider the simple stellar population (SSP) or an instantaneous burst model. The SSP model is physically difficult to understand our $z \ge 6$ galaxies. These galaxies have strong Ly$\alpha$ emission, and presumably have strong nebular emission (we will discuss this later). However, the nebular emission produced by the SSP drops rapidly in the first few Myr, as the instantaneous burst goes off.

From SED modeling, we mainly constrain three physical quantities, including dust reddening $E(B-V)$, stellar mass $M_*$, and age. We use the reddening law of Calzetti et al. (2000), and allow $E(B-V)$ to range between 0.00 and 0.50 in steps of 0.02. The age provided by the GALEV models starts from 4 Myr in steps of 4 Myr. We allow age to vary between 4 and 1000 Myr, which approaches the age of the universe at $z \sim 6$. For $z \ge 6$ galaxies with the limited number of photometric points, age is usually poorly constrained due to various degeneracies among parameters. On the contrary, mass is the amplitude of a model spectrum, and is thus thought to be more easily constrained (e.g., Papovich et al. 2001; Shapley et al. 2005; Conroy 2013; Mobasher et al. 2015). In reality, the mass estimate depends on the mass-to-light ratio ($M/L$), which in turn depends strongly on the age of stellar populations. The measurement of the $M/L$ in most normal-SED galaxies can be accurate to a level of ~0.3 dex, if the uncertainty is dominated by systematics (Conroy 2013). However, the $M/L$ can be very uncertain for certain types of galaxies, or galaxies in certain age ranges (see Conroy 2013 for a review).

The GALEV models provide options to include (or not) nebular line and continuum emission. We use both options for each individual model, and denote models with and without nebular/gaseous emission as “EM” and “NoEM” models, respectively. In the rest of the paper, we mostly discuss the 27 galaxies that were detected in the IRAC 1 band (galaxy parameters cannot be properly constrained for our galaxies without IRAC 1 detections). We fit the models to the SEDs of the galaxies and derive the above parameters by the minimum $\chi^2$ method. These results are shown in Table 2. The $1\sigma$ uncertainties quoted in the table are estimated in a standard way, i.e., allowing all other parameters to vary until $\Delta\chi^2 = 1$. Figure 1 compares the distributions of stellar mass, age, and $E(B-V)$ in the 27 galaxies for different models and SFHs. The two columns of panels from left to right in Figure 1 correspond to the dSFH and the rSFH models, respectively. The gray filled histograms represent the EM models, and the black unfilled histograms represent the NoEM models.

Table 2 and Figure 1 show that the EM and NoEM models produce very different stellar populations for the same galaxies. Stellar populations of many galaxies from the EM models are very young, with ages of several Myr. They also have relatively low stellar masses. On the contrary, stellar populations from the NoEM models are mostly older than a few hundred Myr. They are usually very massive, with masses close to or higher than $10^{10} M_\odot$. Based on the minimum $\chi^2$ values, the quality of the EM models is comparable to (in some cases marginally better than) that of the NoEM models. The results are expected (e.g., Schaerer & de Barros 2009), due to the degeneracy between young galaxies with prominent nebular emission and old galaxies with strong Balmer breaks. This degeneracy is particularly strong for the galaxies in our sample for the following reasons. First, our galaxies are at $5.7 < z < 6.6$, so both IRAC 1 and 2 bands cover some of the strongest emission

lines such as [O III], H$\beta$, and H$\alpha$, etc. The second reason is that the IRAC 1 and 2 data have the largest photometric uncertainties compared to other bands; many galaxies were even not detected in the IRAC 2 band. Finally, for this redshift range, the wavelength range that IRAC 1 and 2 cover mimics the Balmer continuum/break that the synthesis models rely on to constrain stellar populations. The combination of these reasons make the EM and NoEM models indistinguishable for our galaxies.

Figure 2 illustrates the degeneracy mentioned above by the SED modeling of three galaxies, one LAE at $z = 5.664$ (No. 3), one LBG at $z = 6.315$ (No. 34), and one LAE at $z = 6.543$ (No. 44). The two columns of the panels from left to right correspond to the dSFH and the rSFH models, respectively. The red points with error bars are the observed photometric data points. The black and blue SED profiles represent the EM and NoEM models. Note that the GALEV models do not include the Ly$\alpha$ emission line. IGM absorption has been applied to the model spectra. We calculate the IGM absorption using the method of Fan et al. (2001) and Jiang et al. (2008). Note that the Ly$\alpha$ emission or the IGM absorption does not affect our SED fitting, because we did not use the bands that cover Ly$\alpha$. The figure shows that the existence of emission lines in the EM models mimics the Balmer continua/breaks in the NoEM models, which significantly reduces the age and stellar mass of stellar populations needed in the EM models.

## 3.2. Nebular Emission Estimated from the Ly$\alpha$ Line Emission

With ongoing star formation in the models with rising and dSFHs, nebular emission is naturally expected in our galaxies. In particular, we have seen strong Ly$\alpha$ emission lines in these galaxies. However, the actual strength of nebular emission, including continuum and line emission, is unknown. In the GALEV models that we use, the nebular emission is metallicity dependent, but its relative strength to the continuum is fixed in individual model galaxies. In other words, at any given age and metallicity for the same SFH, all galaxies have the same strength of nebular emission. This is a model assumption under certain physical conditions (electron temperature, atomic density, etc.). Real galaxies could have a wide range of nebular emission strength due to different physical states of the gas, including geometry.

In this subsection, we take a more realistic approach to estimate the strength of nebular emission using the observed Ly$\alpha$ line flux in our galaxy sample. One advantage of our sample, other than the available spectroscopic redshifts, is the known Ly$\alpha$ flux. Here we take advantage of this to estimate nebular emission and incorporate the estimated nebular emission into our galaxy models. We first estimate the intrinsic Ly$\alpha$ flux for each galaxy based on the observed Ly$\alpha$ flux given in Paper I. Ly$\alpha$ emission is complicated and largely reduced by the resonant scattering of Ly$\alpha$ photons and the neutral IGM absorption. It is difficult to model Ly$\alpha$ radiative transfer and predict the intrinsic Ly$\alpha$ emission at high redshift. For a given apparent Ly$\alpha$ luminosity, its intrinsic luminosity could have a range of values. Zheng et al. (2010) showed that the distribution of the intrinsic to observed flux ratio ($f_{\text{Ly}\alpha}^{\text{int}}/f_{\text{Ly}\alpha}^{\text{obs}}$) peaks at ~4 in relatively bright galaxies at $z \sim 5.7$. We thus assume $f_{\text{Ly}\alpha}^{\text{int}}/f_{\text{Ly}\alpha}^{\text{obs}} = 4$ for our galaxies. As we already mentioned, the GALEV models do not include the Ly$\alpha$ emission line. We link the Ly$\alpha$ line flux to the model nebular line flux using H$\beta$ under Case B recombination, where the Ly$\alpha$
<table><caption>**Table 2**
Comparison between the SED-fitting Results from the NoEM Model and the GALEV-EM Model</caption><thead><tr><th rowspan="2">No.</th><th colspan="2">M<sub>*</sub> (10<sup>8</sup> M<sub>☉</sub>)</th><th colspan="2">Age (Myr)</th><th colspan="2">E(B - V)</th><th colspan="2">χ<sup>2</sup><sub>r</sub></th></tr><tr><th>NoEM</th><th>EM</th><th>NoEM</th><th>EM</th><th>NoEM</th><th>EM</th><th>NoEM</th><th>EM</th></tr></thead><tbody><tr><td>3</td><td>62.4<sup>+32.0</sup><sub>-24.8</sub></td><td>18.1<sup>+0.9</sup><sub>-0.8</sub></td><td>200<sup>+132</sup><sub>-104</sub></td><td>40<sup>+152</sup><sub>-36</sub></td><td>0.04<sup>+0.04</sup><sub>-0.04</sub></td><td>0.02<sup>+0.06</sup><sub>-0.02</sub></td><td>1.4</td><td>1.5</td></tr><tr><td>4</td><td>114.0<sup>+11.0</sup><sub>-10.0</sub></td><td>3.5<sup>+0.2</sup><sub>-0.2</sub></td><td>928<sup>+68</sup><sub>-316</sub></td><td>4<sup>+8</sup><sub>-0</sub></td><td>0.00<sup>+0.04</sup><sub>-0.00</sub></td><td>0.00<sup>+0.04</sup><sub>-0.00</sub></td><td>23.6</td><td>18.8</td></tr><tr><td>15</td><td>368.8<sup>+35.6</sup><sub>-16.6</sub></td><td>15.7<sup>+0.7</sup><sub>-0.7</sub></td><td>976<sup>+20</sup><sub>-84</sub></td><td>4<sup>+8</sup><sub>-0</sub></td><td>0.04<sup>+0.04</sup><sub>-0.04</sub></td><td>0.08<sup>+0.04</sup><sub>-0.04</sub></td><td>4.3</td><td>2.7</td></tr><tr><td>20</td><td>166.0<sup>+63.2</sup><sub>-45.8</sub></td><td>44.1<sup>+77.4</sup><sub>-25.7</sub></td><td>656<sup>+340</sup><sub>-404</sub></td><td>88<sup>+420</sup><sub>-84</sub></td><td>0.14<sup>+0.06</sup><sub>-0.06</sub></td><td>0.16<sup>+0.06</sup><sub>-0.06</sub></td><td>0.2</td><td>0.2</td></tr><tr><td>23</td><td>127.5<sup>+12.3</sup><sub>-5.7</sub></td><td>7.3<sup>+0.3</sup><sub>-0.3</sub></td><td>988<sup>+8</sup><sub>-4</sub></td><td>8<sup>+76</sup><sub>-276</sub></td><td>0.00<sup>+0.04</sup><sub>-0.00</sub></td><td>0.04<sup>+0.04</sup><sub>-0.04</sub></td><td>1.8</td><td>1.2</td></tr><tr><td>24</td><td>86.1<sup>+12.8</sup><sub>-11.1</sub></td><td>3.8<sup>+0.2</sup><sub>-0.2</sub></td><td>632<sup>+364</sup><sub>-264</sub></td><td>4<sup>+16</sup><sub>-0</sub></td><td>0.00<sup>+0.04</sup><sub>-0.00</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>2.4</td><td>0.8</td></tr><tr><td>25</td><td>356.2<sup>+34.4</sup><sub>-16.0</sub></td><td>14.4<sup>+0.7</sup><sub>-0.6</sub></td><td>992<sup>+4</sup><sub>-52</sub></td><td>4<sup>+8</sup><sub>-0</sub></td><td>0.08<sup>+0.04</sup><sub>-0.04</sub></td><td>0.12<sup>+0.04</sup><sub>-0.04</sub></td><td>17.1</td><td>13.3</td></tr><tr><td>27</td><td>93.1<sup>+35.4</sup><sub>-16.6</sub></td><td>3.2<sup>+0.3</sup><sub>-0.3</sub></td><td>988<sup>+8</sup><sub>-468</sub></td><td>4<sup>+40</sup><sub>-0</sub></td><td>0.06<sup>+0.06</sup><sub>-0.06</sub></td><td>0.08<sup>+0.08</sup><sub>-0.08</sub></td><td>3.4</td><td>2.4</td></tr><tr><td>28</td><td>56.8<sup>+41.9</sup><sub>-9.6</sub></td><td>2.8<sup>+0.3</sup><sub>-0.4</sub></td><td>536<sup>+460</sup><sub>-364</sub></td><td>4<sup>+616</sup><sub>-0</sub></td><td>0.00<sup>+0.08</sup><sub>-0.00</sub></td><td>0.00<sup>+0.08</sup><sub>-0.00</sub></td><td>0.4</td><td>&lt;0.1</td></tr><tr><td>30</td><td>86.9<sup>+71.2</sup><sub>-21.0</sub></td><td>3.0<sup>+1.0</sup><sub>-1.8</sub></td><td>992<sup>+4</sup><sub>-448</sub></td><td>4<sup>+24</sup><sub>-0</sub></td><td>0.02<sup>+0.10</sup><sub>-0.02</sub></td><td>0.04<sup>+0.08</sup><sub>-0.04</sub></td><td>3.7</td><td>2.7</td></tr><tr><td>31</td><td>275.6<sup>+41.5</sup><sub>-93.5</sub></td><td>8.3<sup>+3.7</sup><sub>-5.0</sub></td><td>984<sup>+12</sup><sub>-384</sub></td><td>4<sup>+16</sup><sub>-0</sub></td><td>0.16<sup>+0.08</sup><sub>-0.08</sub></td><td>0.16<sup>+0.08</sup><sub>-0.08</sub></td><td>1.8</td><td>1.0</td></tr><tr><td>34</td><td>852.2<sup>+82.2</sup><sub>-38.4</sub></td><td>41.0<sup>+1.9</sup><sub>-1.8</sub></td><td>992<sup>+4</sup><sub>-136</sub></td><td>8<sup>+12</sup><sub>-4</sub></td><td>0.14<sup>+0.04</sup><sub>-0.04</sub></td><td>0.16<sup>+0.06</sup><sub>-0.04</sub></td><td>3.2</td><td>1.4</td></tr><tr><td>35</td><td>3314.4<sup>+319.8</sup><sub>-291.6</sub></td><td>102.5<sup>+4.8</sup><sub>-4.6</sub></td><td>976<sup>+20</sup><sub>-112</sub></td><td>4<sup>+8</sup><sub>-0</sub></td><td>0.36<sup>+0.04</sup><sub>-0.04</sub></td><td>0.36<sup>+0.04</sup><sub>-0.04</sub></td><td>68.3</td><td>25.0</td></tr><tr><td>36</td><td>157.7<sup>+15.2</sup><sub>-13.9</sub></td><td>124.9<sup>+12.1</sup><sub>-11.0</sub></td><td>992<sup>+4</sup><sub>-200</sub></td><td>984<sup>+12</sup><sub>-392</sub></td><td>0.02<sup>+0.06</sup><sub>-0.02</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>3.5</td><td>2.9</td></tr><tr><td>43</td><td>129.4<sup>+41.2</sup><sub>-26.6</sub></td><td>4.5<sup>+0.4</sup><sub>-0.4</sub></td><td>992<sup>+4</sup><sub>-236</sub></td><td>4<sup>+12</sup><sub>-0</sub></td><td>0.06<sup>+0.08</sup><sub>-0.06</sub></td><td>0.08<sup>+0.06</sup><sub>-0.08</sub></td><td>6.9</td><td>3.8</td></tr><tr><td>44</td><td>983.9<sup>+94.9</sup><sub>-44.3</sub></td><td>68.0<sup>+3.2</sup><sub>-3.1</sub></td><td>992<sup>+4</sup><sub>-96</sub></td><td>12<sup>+12</sup><sub>-8</sub></td><td>0.16<sup>+0.04</sup><sub>-0.04</sub></td><td>0.20<sup>+0.04</sup><sub>-0.04</sub></td><td>10.0</td><td>2.3</td></tr><tr><td>47</td><td>758.3<sup>+53.4</sup><sub>-127.6</sub></td><td>388.1<sup>+123.5</sup><sub>-79.8</sub></td><td>440<sup>+60</sup><sub>-80</sub></td><td>312<sup>+96</sup><sub>-84</sub></td><td>0.18<sup>+0.04</sup><sub>-0.04</sub></td><td>0.14<sup>+0.04</sup><sub>-0.04</sub></td><td>7.1</td><td>6.6</td></tr><tr><td>49</td><td>207.8<sup>+121.6</sup><sub>-42.7</sub></td><td>6.1<sup>+0.6</sup><sub>-3.7</sub></td><td>996<sup>+0</sup><sub>-436</sub></td><td>4<sup>+992</sup><sub>-0</sub></td><td>0.12<sup>+0.08</sup><sub>-0.08</sub></td><td>0.12<sup>+0.08</sup><sub>-0.08</sub></td><td>0.4</td><td>0.1</td></tr><tr><td>50</td><td>1223.8<sup>+389.5</sup><sub>-251.7</sub></td><td>30.6<sup>+7.9</sup><sub>-6.3</sub></td><td>996<sup>+0</sup><sub>-324</sub></td><td>4<sup>+12</sup><sub>-0</sub></td><td>0.36<sup>+0.08</sup><sub>-0.06</sub></td><td>0.34<sup>+0.08</sup><sub>-0.06</sub></td><td>1.5</td><td>0.5</td></tr><tr><td>54</td><td>144.7<sup>+84.7</sup><sub>-29.8</sub></td><td>112.4<sup>+42.8</sup><sub>-14.5</sub></td><td>992<sup>+4</sup><sub>-268</sub></td><td>972<sup>+24</sup><sub>-412</sub></td><td>0.02<sup>+0.08</sup><sub>-0.02</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>1.9</td><td>1.6</td></tr><tr><td>58</td><td>183.4<sup>+81.7</sup><sub>-56.5</sub></td><td>178.8<sup>+56.9</sup><sub>-43.2</sub></td><td>480<sup>+332</sup><sub>-216</sub></td><td>884<sup>+112</sup><sub>-380</sub></td><td>0.12<sup>+0.06</sup><sub>-0.06</sub></td><td>0.06<sup>+0.06</sup><sub>-0.06</sub></td><td>1.9</td><td>3.1</td></tr><tr><td>61</td><td>875.4<sup>+129.7</sup><sub>-77.0</sub></td><td>29.2<sup>+1.4</sup><sub>-1.3</sub></td><td>984<sup>+12</sup><sub>-116</sub></td><td>4<sup>+8</sup><sub>-0</sub></td><td>0.20<sup>+0.04</sup><sub>-0.04</sub></td><td>0.22<sup>+0.04</sup><sub>-0.04</sub></td><td>16.4</td><td>9.3</td></tr><tr><td>62</td><td>194.5<sup>+61.9</sup><sub>-47.0</sub></td><td>163.2<sup>+62.1</sup><sub>-45.0</sub></td><td>756<sup>+240</sup><sub>-248</sub></td><td>656<sup>+304</sup><sub>-256</sub></td><td>0.04<sup>+0.06</sup><sub>-0.04</sub></td><td>0.04<sup>+0.06</sup><sub>-0.04</sub></td><td>1.6</td><td>1.5</td></tr><tr><td>63</td><td>201.7<sup>+19.5</sup><sub>-9.1</sub></td><td>7.0<sup>+0.3</sup><sub>-0.3</sub></td><td>988<sup>+8</sup><sub>-148</sub></td><td>4<sup>+8</sup><sub>-0</sub></td><td>0.00<sup>+0.04</sup><sub>-0.00</sub></td><td>0.02<sup>+0.04</sup><sub>-0.02</sub></td><td>6.3</td><td>2.9</td></tr><tr><td>64</td><td>167.1<sup>+16.1</sup><sub>-14.7</sub></td><td>5.0<sup>+0.2</sup><sub>-0.2</sub></td><td>992<sup>+4</sup><sub>-200</sub></td><td>4<sup>+8</sup><sub>-0</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>6.5</td><td>3.8</td></tr><tr><td>66</td><td>113.5<sup>+102.8</sup><sub>-53.9</sub></td><td>23.0<sup>+56.8</sup><sub>-15.0</sub></td><td>920<sup>+76</sup><sub>-636</sub></td><td>96<sup>+900</sup><sub>-92</sub></td><td>0.10<sup>+0.10</sup><sub>-0.10</sub></td><td>0.12<sup>+0.10</sup><sub>-0.10</sub></td><td>&lt;0.1</td><td>&lt;0.1</td></tr><tr><td>67</td><td>774.8<sup>+74.8</sup><sub>-68.2</sub></td><td>27.1<sup>+1.3</sup><sub>-1.2</sub></td><td>976<sup>+20</sup><sub>-108</sub></td><td>4<sup>+8</sup><sub>-0</sub></td><td>0.10<sup>+0.04</sup><sub>-0.04</sub></td><td>0.12<sup>+0.04</sup><sub>-0.04</sub></td><td>27.2</td><td>3.3</td></tr></tbody></table>

Note. This table compares the SED-fitting results from the NoEM Model and from the GALEV-EM Model with rSFH. The minimum and maximum values of age are 4 and 1000 Myr. When an age is close to the two limits, its errors are not reliable because of the lack of dynamic range. In addition, the fitting results with very large or small $\chi_r^2$ ($\chi_r^2 \gg 1$ or $\chi_r^2 \ll 1$) are not reliable.

to Hβ ratio ($f_{Lyα}^{\text{int}}/f_{Hβ}^{\text{int}}$) is roughly 25 (Lyα/Hα = 8.7 and Hα/Hβ = 2.87). The combination of the above two steps leads to an assumption that the ratio of the observed Lyα to the intrinsic Hβ flux ratio ($f_{Lyα}^{\text{obs}}/f_{Hβ}^{\text{int}}$) is 6.25, in the case of no dust extinction.

In order to incorporate this ratio to the GALEV models, we take each pair of the EM and NoEM model spectra in the whole parameter space from the previous subsection, and compute their nebular emission by subtracting the NoEM spectrum from the EM spectrum. We then scale the nebular emission (both line and continuum) to match the $f_{Lyα}^{\text{obs}}/f_{Hβ}^{\text{int}}$ ratio (this ratio varies with dust extinction, and it is taken into account in this step). The details are as follows. We first take the Lyα EW from Paper I. Based on the continuum of the model spectrum, the Lyα EW, and the $f_{Lyα}^{\text{obs}}/f_{Hβ}^{\text{int}}$ ratio, we calculate the Hβ flux and EW for the model spectrum. Then the whole nebular emission spectrum is scaled and added to the model (continuum) spectrum. This model spectrum will have the required $f_{Lyα}^{\text{obs}}/f_{Hβ}^{\text{int}}$ ratio during the SED fitting procedure. It does not guarantee that we will obtain the observed Lyα flux, because we only reserve its EW, the relative line flux to the continuum. Fortunately, there are no strong lines in the wavelength range covered by the non-IRAC data points (we did not use the bands that cover Lyα). Therefore, the

continuum at the wavelength of Lyα is well modeled (see Figures 2 and 3). It ensures that the Lyα flux in the best-fitting model spectrum is very close to the observed flux. This has been tested for each object, and it is accurate to a few percent.

We next perform the SED modeling with the new EM models. We refer to these new EM models as “Lyα-EM” models, and the original GALEV EM models as “GALEV-EM” models. Table 3 shows the results for the “Lyα-EM” models with dSFH and rSFH. Column 1 lists the galaxy sequence numbers. Columns 2–3 are stellar masses $M_*$ derived for the dSFH and the rSFH models, respectively. Columns 4–5 and Columns 6–7 are the derived ages and dust extinction for the two SFHs. Columns 8–9 give $\chi^2$ from the modeling. Figure 3 illustrates the new SED modeling of the three galaxies shown in Figure 2. The blue profiles are the NoEM model spectra shown as a reference. The Lyα-EM model spectra are shown in black. The Lyα-EM model spectra for objects No. 3 and No. 44 still show strong nebular emission lines, as seen in Figure 2. The spectra for object No. 34, however, have much weaker lines compared to those in Figure 2. This is due to its weak Lyα emission line (its rest-frame EW is only 8.6 Å). Figure 4 compares the distributions of stellar mass, age, and $E(B-V)$ for the Lyα-EM and NoEM models. These distributions are roughly consistent with those for the GALEV-EM models
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Distributions of the stellar masses, ages, and $E(B-V)$ estimated for the 27 galaxies with the IRAC 1 detections. The two columns of the panels from the left to the right correspond to the declining SFH and the rising SFH, respectively. The gray filled histograms represent the EM models, and the black unfilled histograms represent the NoEM models. The NoEM and EM models produce different stellar populations (old and massive vs. young and less massive) for the same galaxies. This is due to the strong degeneracy between young galaxies with prominent nebular emission and old galaxies with strong Balmer breaks.</figcaption></figure>

shown in Figure 1. We will discuss this in detail in the next
section.

## 4. RESULTS

For galaxies with strong nebular emission, the results of SED
fitting strongly depend on the relative strength of nebular
emission. In the above section, we used two methods to include
nebular emission. One was to use the GALEV models with
nebular emission (the GALEV-EM models), and the other one
was to scale the nebular emission in the GALEV models to
match our observed Lyα flux (the Lyα-EM models). The latter
one is a more realistic approach, and consists of two
assumptions: the ratio of the observed to intrinsic Lyα flux
$f_{\text{Ly}\alpha}^{\text{int}}/f_{\text{Ly}\alpha}^{\text{obs}}$ is 4, and the intrinsic Lyα to Hβ flux ratio $f_{\text{Ly}\alpha}^{\text{int}}/f_{\text{H}\beta}^{\text{int}}$
is 25 (Case B), so that $f_{\text{Ly}\alpha}^{\text{obs}}/f_{\text{H}\beta}^{\text{int}} = 6.25$. Each assumption
involves non-negligible uncertainties. For the first assumption,
the probability distribution function of $f_{\text{Ly}\alpha}^{\text{int}}/f_{\text{Ly}\alpha}^{\text{obs}}$ depends on
the physical conditions of galaxies (Zheng et al. 2010). In
addition, it is apparently a function of redshift: it increases
(observed Lyα flux decreases) toward higher redshifts as the
neutral fraction of the IGM increases. For the second
assumption, real galaxies are more complex than the assumed
ideal Case B recombination.

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Three examples of SED modeling illustrating the strong degeneracy between young galaxies with prominent nebular emission and older galaxies with strong Balmer breaks. The two columns of the panels correspond to the declining SFH and the rising SFH. The red points with error bars are the observed photometric data points. The horizontal errors indicate the wavelength ranges of the filters. The black and blue profiles represent the EM and NoEM models, respectively. The large green crosses represent the photometric points predicted by the EM models. The existence of emission lines in the EM models exactly mimics the the Balmer continua/breaks in the NoEM models, which significantly reduces the required ages and stellar masses of stellar populations in the EM models.</figcaption></figure>

Despite the uncertainties mentioned above, our assumptions
are reasonable (see more details in the discussion section). As
we will see, the Lyα-EM and GALEV-EM models produce
roughly similar results on average, although there are large
object-to-object variations due to the large range of the
observed Lyα EW. In order to explore the possible ranges of
physical parameters, we vary the $f_{\text{Ly}\alpha}^{\text{obs}}/f_{\text{H}\beta}^{\text{int}}$ ratio by a factor of
two, i.e., we perform another two sets of SED modeling by
assuming $f_{\text{Ly}\alpha}^{\text{obs}}/f_{\text{H}\beta}^{\text{int}} = 3.125$ (referred to as EM-strong mod-
els), and $f_{\text{Ly}\alpha}^{\text{obs}}/f_{\text{H}\beta}^{\text{int}} = 12.5$ (referred to as EM-weak models).
Based on the simulations of Zheng et al. (2010), the majority of
galaxies are included in the range considered here. We will
discuss nebular emission and our assumptions in greater detail
in Section 5.

### 4.1. Stellar Mass

It is worth briefly discussing the bias of our sample before
we discuss the resulting stellar masses. The bias from sample
selection was discussed in Paper I, and in Section 2 of this
paper. The 42 galaxies in Table 1 are brighter than ~27.2
AB mag in J. The 27 galaxies shown in Figures 1 and 4 have
detections in the IRAC 1 band (see also Table 1), so they are
further limited by the 3.6 µm flux detection limit. Because the
3.6 µm flux is closely related to stellar mass for z ~ 6 galaxies,
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Three examples of SED fitting with the Lyα-EM models. The three galaxies are the same as shown in Figure 2. The red points with error bars are the observed photometric data points. The blue profiles are the NoEM model spectra shown as a reference (same as those in Figure 2). The Lyα-EM model spectra are shown in black. The large green crosses represent the photometric points predicted by the EM models. The Lyα-EM model spectra of No. 3 and No. 44 show strong nebular emission lines as seen in Figure 2. The spectra of No. 34 show much weaker lines compared to those in Figure 2, due to its weak Lyα emission line.</figcaption></figure>

the sample of the 27 galaxies is biased toward the more massive galaxies. This is the reason that the stellar masses derived from the NoEM models are mostly close to or higher than $10^{10} M_{\odot}$ (see the histograms in Figures 1 and 4).

Stellar mass is thought to be the parameter that is least sensitive to model assumptions, because it is directly measured from the scaling of galaxy model spectra (e.g., Papovich et al. 2001; Shapley et al. 2005). When nebular emission is included, however, the measurement of stellar mass becomes less straightforward. Figure 5 compares the stellar masses derived from models with different nebular emission. The filled gray histograms show the stellar mass distributions from the Lyα-EM models. The black unfilled histograms represent the mass distributions from the EM-weak models (top two panels) and the EM-strong models (middle two panels). The distributions of the EM-weak and EM-strong models do not significantly deviate from those for the Lyα-EM models. The bottom panel shows the ratios of the stellar masses derived from different EM models. Although the mass ratios of the GALEV-EM to Lyα-EM models span a wide range due to the wide distribution of the Lyα EW in our sample, the median ratio is close to 1, suggesting that the assumptions we made for the Lyα-EM models are reasonable. We use the EM-weak and EM-strong models to explore the possible mass ranges. As expected, they show that models with stronger (weaker) nebular emission generally produce stellar populations with lower (higher) masses and younger (older) ages. On the other

hand, the results for most galaxies from the three EM models (Lyα-EM, EM-weak, and EM-strong) are roughly consistent. The bottom panel shows that about 70%–80% of galaxies in our sample have similar stellar masses (within a factor of three) from the three different EM models.

The upper panel of Figure 6 shows the relation between the 3.6 µm flux and stellar mass derived from models with rSFH. There is a tight relation for the NoEM models (gray pluses), as explained above. For the EM models (the Lyα-EM models here; crosses), the correlation still exists for the whole sample, but with a much larger scatter. The relation also significantly deviates from the relation for the NoEM models, due to the contamination of strong nebular emission in the IRAC 1 band. We divide these galaxies into two subsamples: a “young” subsample (age <30 Myr; blue crosses) and an “old” subsample (age >30 Myr; red crosses) based on the age distributions shown in Figures 4 and 8. The old subsample consists of galaxies with ages of several hundred Myr, and the galaxies in the young subsample are usually younger than 30 Myr (see Section 4.2). It is obvious that the relation between the 3.6 µm flux and stellar mass is tight in either subsample. The linear Pearson correlation coefficients for the old and young subsamples are −0.87 and −0.96, respectively. The upper panel of Figure 7 shows the relation for dSFH, which is also obvious though less tight. The linear Pearson correlation coefficients for the two subsamples are −0.66 and −0.92, respectively. This relation is partially shaped by selection effects, i.e., these objects were selected in a relatively small parameter space. It may also be partially due to the following reasons. For the old subsample, the 3.6 µm flux is still dominated by stellar emission, and thus reflects the stellar mass. In the young subsample, nebular emission has a large contribution to the 3.6 µm flux. On the other hand, these galaxies are within a small range of (very young) age, and their relative strength of the nebular to stellar emission does not span a wide range (note that it is fixed for any GALEV-EM model SED). So the combination of the nebular and stellar emission still follows the correlation, though the presence of nebular emission has largely reduced the required stellar masses.

The lower panel in Figure 6 shows the relation between stellar mass and rest-frame UV luminosity $M_{1500}$, the absolute AB mag at 1500 Å derived in Paper I. The $M_{1500}$ values have been corrected for dust extinction. As for the upper panel, the pluses and crosses represent the NoEM and EM models, and the red and blue crosses represent the old and young subsamples, respectively. For the NoEM models, there is a weak correlation in which stellar masses are higher in more luminous galaxies. Such correlation is more tight in the two subsamples in the Lyα-EM models. The linear Pearson correlation coefficients for the old and young subsamples are −0.92 and −0.99, respectively. The lower panel of Figure 7 shows the relation for dSFH, which is also obvious though less tight. The linear Pearson correlation coefficients for the two subsamples are −0.67 and −0.96, respectively. This relation has been reported in previous studies of high-redshift galaxies with different SFHs, including constant SFH and slowly varying (rising or declining) SFHs (e.g., Stark et al. 2009; González et al. 2011; McLure et al. 2011; Curtis-Lake et al. 2013). It is believed to be the high-redshift version of the mass–SFR relation (or the “so-called” main sequence of star-forming galaxies) found at lower redshifts, in which the SFR is higher in more massive galaxies, and the normalization of the relation at
<table><caption>**Table 3**
SED-fitting Results with the Lyα-EM Models</caption><thead><tr><th rowspan="2">No.</th><th colspan="2">M<sub>*</sub> (10<sup>8</sup> M<sub>☉</sub>)</th><th colspan="2">Age (Myr)</th><th colspan="2">E(B - V)</th><th colspan="2">χ<sup>2</sup><sub>r</sub></th></tr><tr><th>dSFH</th><th>rSFH</th><th>dSFH</th><th>rSFH</th><th>dSFH</th><th>rSFH</th><th>dSFH</th><th>rSFH</th></tr></thead><tbody><tr><td>3</td><td>7.2<sup>+0.3</sup><sub>-0.3</sub></td><td>7.2<sup>+0.3</sup><sub>-0.3</sub></td><td>4<sup>+12</sup><sub>-0</sub></td><td>4<sup>+20</sup><sub>-0</sub></td><td>0.04<sup>+0.06</sup><sub>-0.04</sub></td><td>0.04<sup>+0.06</sup><sub>-0.04</sub></td><td>1.5</td><td>1.4</td></tr><tr><td>4</td><td>3.4<sup>+0.2</sup><sub>-0.1</sub></td><td>3.4<sup>+0.2</sup><sub>-0.1</sub></td><td>4<sup>+8</sup><sub>-0</sub></td><td>4<sup>+8</sup><sub>-0</sub></td><td>0.00<sup>+0.04</sup><sub>-0.00</sub></td><td>0.00<sup>+0.04</sup><sub>-0.00</sub></td><td>18.1</td><td>18.1</td></tr><tr><td>15</td><td>127.3<sup>+12.3</sup><sub>-11.2</sub></td><td>21.1<sup>+1.0</sup><sub>-0.9</sub></td><td>168<sup>+32</sup><sub>-24</sub></td><td>8<sup>+12</sup><sub>-4</sub></td><td>0.00<sup>+0.04</sup><sub>-0.00</sub></td><td>0.08<sup>+0.04</sup><sub>-0.04</sub></td><td>2.8</td><td>2.9</td></tr><tr><td>20</td><td>29.5<sup>+13.1</sup><sub>-9.1</sub></td><td>32.8<sup>+21.6</sup><sub>-10.1</sub></td><td>148<sup>+164</sup><sub>-128</sub></td><td>492<sup>+504</sup><sub>-468</sub></td><td>0.02<sup>+0.06</sup><sub>-0.02</sub></td><td>0.02<sup>+0.08</sup><sub>-0.02</sub></td><td>1.5</td><td>1.6</td></tr><tr><td>23</td><td>46.8<sup>+9.5</sup><sub>-6.0</sub></td><td>57.9<sup>+11.7</sup><sub>-9.8</sub></td><td>124<sup>+52</sup><sub>-48</sub></td><td>444<sup>+336</sup><sub>-208</sub></td><td>0.00<sup>+0.04</sup><sub>-0.00</sub></td><td>0.00<sup>+0.04</sup><sub>-0.00</sub></td><td>1.7</td><td>1.5</td></tr><tr><td>24</td><td>3.8<sup>+0.2</sup><sub>-0.2</sub></td><td>3.8<sup>+0.2</sup><sub>-0.2</sub></td><td>4<sup>+20</sup><sub>-0</sub></td><td>4<sup>+28</sup><sub>-0</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>0.8</td><td>0.8</td></tr><tr><td>25</td><td>20.9<sup>+1.0</sup><sub>-0.9</sub></td><td>66.5<sup>+6.4</sup><sub>-5.8</sub></td><td>12<sup>+20</sup><sub>-8</sub></td><td>956<sup>+40</sup><sub>-184</sub></td><td>0.10<sup>+0.04</sup><sub>-0.04</sub></td><td>0.02<sup>+0.04</sup><sub>-0.02</sub></td><td>13.2</td><td>12.8</td></tr><tr><td>27</td><td>12.6<sup>+0.5</sup><sub>-1.6</sub></td><td>17.4<sup>+14.2</sup><sub>-4.2</sub></td><td>88<sup>+160</sup><sub>-84</sub></td><td>324<sup>+672</sup><sub>-320</sub></td><td>0.00<sup>+0.08</sup><sub>-0.00</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>2.2</td><td>2.2</td></tr><tr><td>28</td><td>2.7<sup>+0.3</sup><sub>-0.2</sub></td><td>2.7<sup>+0.3</sup><sub>-0.2</sub></td><td>4<sup>+20</sup><sub>-0</sub></td><td>4<sup>+28</sup><sub>-0</sub></td><td>0.00<sup>+0.04</sup><sub>-0.00</sub></td><td>0.00<sup>+0.04</sup><sub>-0.00</sub></td><td>1.9</td><td>1.9</td></tr><tr><td>30</td><td>43.4<sup>+28.6</sup><sub>-8.9</sub></td><td>26.8<sup>+15.7</sup><sub>-4.5</sub></td><td>192<sup>+184</sup><sub>-136</sub></td><td>912<sup>+84</sup><sub>-728</sub></td><td>0.00<sup>+0.08</sup><sub>-0.00</sub></td><td>0.00<sup>+0.08</sup><sub>-0.00</sub></td><td>3.2</td><td>3.1</td></tr><tr><td>31</td><td>2.4<sup>+0.6</sup><sub>-1.7</sub></td><td>2.4<sup>+0.6</sup><sub>-2.2</sub></td><td>4<sup>+168</sup><sub>-0</sub></td><td>4<sup>+668</sup><sub>-0</sub></td><td>0.04<sup>+0.08</sup><sub>-0.04</sub></td><td>0.04<sup>+0.08</sup><sub>-0.04</sub></td><td>0.2</td><td>0.2</td></tr><tr><td>34</td><td>640.7<sup>+129.6</sup><sub>-107.8</sub></td><td>250.0<sup>+24.1</sup><sub>-11.2</sub></td><td>356<sup>+44</sup><sub>-44</sub></td><td>992<sup>+4</sup><sub>-204</sub></td><td>0.06<sup>+0.04</sup><sub>-0.04</sub></td><td>0.12<sup>+0.04</sup><sub>-0.04</sub></td><td>1.7</td><td>2.1</td></tr><tr><td>35</td><td>54.2<sup>+2.5</sup><sub>-2.4</sub></td><td>81.1<sup>+3.8</sup><sub>-3.7</sub></td><td>592<sup>+68</sup><sub>-8</sub></td><td>4<sup>+8</sup><sub>-0</sub></td><td>0.00<sup>+0.04</sup><sub>-0.00</sub></td><td>0.34<sup>+0.04</sup><sub>-0.06</sub></td><td>1.1</td><td>2.8</td></tr><tr><td>36</td><td>109.8<sup>+28.4</sup><sub>-18.5</sub></td><td>46.4<sup>+4.5</sup><sub>-4.1</sub></td><td>240<sup>+72</sup><sub>-60</sub></td><td>988<sup>+8</sup><sub>-236</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>3.3</td><td>3.1</td></tr><tr><td>43</td><td>76.3<sup>+34.0</sup><sub>-15.7</sub></td><td>30.8<sup>+9.8</sup><sub>-7.4</sub></td><td>284<sup>+128</sup><sub>-120</sub></td><td>972<sup>+24</sup><sub>-528</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>0.02<sup>+0.06</sup><sub>-0.02</sub></td><td>4.4</td><td>4.4</td></tr><tr><td>44</td><td>331.1<sup>+31.9</sup><sub>-29.1</sub></td><td>172.9<sup>+16.7</sup><sub>-15.2</sub></td><td>260<sup>+32</sup><sub>-32</sub></td><td>872<sup>+124</sup><sub>-160</sub></td><td>0.06<sup>+0.04</sup><sub>-0.04</sub></td><td>0.08<sup>+0.04</sup><sub>-0.04</sub></td><td>1.1</td><td>1.0</td></tr><tr><td>47</td><td>368.7<sup>+17.4</sup><sub>-16.6</sub></td><td>391.1<sup>+101.3</sup><sub>-107.8</sub></td><td>76<sup>+20</sup><sub>-20</sub></td><td>184<sup>+68</sup><sub>-40</sub></td><td>0.18<sup>+0.04</sup><sub>-0.04</sub></td><td>0.18<sup>+0.04</sup><sub>-0.04</sub></td><td>5.8</td><td>5.9</td></tr><tr><td>49</td><td>59.3<sup>+30.5</sup><sub>-21.9</sub></td><td>35.6<sup>+35.4</sup><sub>-16.9</sub></td><td>196<sup>+144</sup><sub>-148</sub></td><td>212<sup>+748</sup><sub>-208</sub></td><td>0.04<sup>+0.08</sup><sub>-0.04</sub></td><td>0.08<sup>+0.08</sup><sub>-0.08</sub></td><td>0.1</td><td>0.1</td></tr><tr><td>50</td><td>173.5<sup>+77.3</sup><sub>-97.8</sub></td><td>9.0<sup>+0.9</sup><sub>-1.9</sub></td><td>528<sup>+40</sup><sub>-72</sub></td><td>4<sup>+548</sup><sub>-0</sub></td><td>0.06<sup>+0.06</sup><sub>-0.06</sub></td><td>0.22<sup>+0.06</sup><sub>-0.08</sub></td><td>&lt;0.1</td><td>&lt;0.1</td></tr><tr><td>54</td><td>36.5<sup>+13.9</sup><sub>-6.2</sub></td><td>51.5<sup>+23.0</sup><sub>-10.6</sub></td><td>112<sup>+100</sup><sub>-84</sub></td><td>460<sup>+524</sup><sub>-392</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>1.9</td><td>1.8</td></tr><tr><td>58</td><td>138.7<sup>+44.1</sup><sub>-33.5</sub></td><td>61.6<sup>+19.6</sup><sub>-14.9</sub></td><td>268<sup>+80</sup><sub>-76</sub></td><td>960<sup>+36</sup><sub>-384</sub></td><td>0.02<sup>+0.06</sup><sub>-0.02</sub></td><td>0.04<sup>+0.06</sup><sub>-0.04</sub></td><td>3.7</td><td>3.6</td></tr><tr><td>61</td><td>357.6<sup>+53.0</sup><sub>-31.5</sub></td><td>66.3<sup>+3.1</sup><sub>-3.0</sub></td><td>516<sup>+24</sup><sub>-24</sub></td><td>32<sup>+80</sup><sub>-28</sub></td><td>0.00<sup>+0.04</sup><sub>-0.00</sub></td><td>0.20<sup>+0.04</sup><sub>-0.04</sub></td><td>8.2</td><td>9.5</td></tr><tr><td>62</td><td>128.0<sup>+33.1</sup><sub>-21.5</sub></td><td>115.0<sup>+51.2</sup><sub>-31.7</sub></td><td>228<sup>+56</sup><sub>-64</sub></td><td>472<sup>+316</sup><sub>-216</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>0.04<sup>+0.06</sup><sub>-0.04</sub></td><td>1.2</td><td>1.2</td></tr><tr><td>63</td><td>9.8<sup>+0.5</sup><sub>-0.4</sub></td><td>9.4<sup>+0.5</sup><sub>-0.4</sub></td><td>8<sup>+20</sup><sub>-4</sub></td><td>8<sup>+36</sup><sub>-4</sub></td><td>0.02<sup>+0.04</sup><sub>-0.02</sub></td><td>0.02<sup>+0.04</sup><sub>-0.02</sub></td><td>3.5</td><td>3.5</td></tr><tr><td>64</td><td>35.2<sup>+9.1</sup><sub>-3.1</sub></td><td>41.3<sup>+30.5</sup><sub>-11.4</sub></td><td>72<sup>+68</sup><sub>-44</sub></td><td>200<sup>+260</sup><sub>-136</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>0.00<sup>+0.06</sup><sub>-0.00</sub></td><td>5.5</td><td>5.4</td></tr><tr><td>66</td><td>1.4<sup>+0.7</sup><sub>-0.8</sub></td><td>1.4<sup>+0.7</sup><sub>-1.2</sub></td><td>4<sup>+244</sup><sub>-0</sub></td><td>4<sup>+992</sup><sub>-0</sub></td><td>0.00<sup>+0.08</sup><sub>-0.00</sub></td><td>0.00<sup>+0.08</sup><sub>-0.00</sub></td><td>0.3</td><td>0.3</td></tr><tr><td>67</td><td>32.0<sup>+1.5</sup><sub>-1.4</sub></td><td>32.0<sup>+1.5</sup><sub>-1.4</sub></td><td>4<sup>+8</sup><sub>-0</sub></td><td>4<sup>+8</sup><sub>-0</sub></td><td>0.14<sup>+0.04</sup><sub>-0.04</sub></td><td>0.14<sup>+0.04</sup><sub>-0.04</sub></td><td>3.2</td><td>3.2</td></tr></tbody></table>

Note. This table shows the SED-fitting results using the Lyα-EM models, i.e., models with nebular emission estimated from the observed Lyα line flux (Section 3.2). The minimum and maximum values of age are 4 and 1000 Myr. When an age is close to the two limits, its errors are not reliable because of the lack of dynamic range.

higher redshift is higher (e.g., Daddi et al. 2007; Elbaz et al. 2007; Noeske et al. 2007). Objects with extreme star-forming activity may significantly deviate from the relation. In our sample, the UV luminosity $M_{1500}$ reflects SFR. We will discuss the SFRs of our sample in Section 4.4. The tight relations seen in the figure could also be partially due to the selection effects mentioned earlier.

From Figures 4–6, all the four EM models show a wide range of stellar masses in our galaxies, ranging from $M_{*} \sim 10^8$ to $10^{11} M_{\odot}$. In particular, a large fraction (~50%) of the galaxies have $M_{*}$ close to, or higher than $5 \times 10^9 M_{\odot}$, suggesting that very massive galaxies already exist when the universe was only ~900 Myr old.

## 4.2. Age

As already mentioned, the age of a stellar population is usually poorly constrained from SED modeling, especially for high-redshift galaxies with only three to five photometric data points. During our SED fitting, we allowed the age to vary between 4 and 1000 Myr. Figures 1 and 4 show that the stellar populations derived from the NoEM models are mostly older than 100 Myr. The inclusion of nebular emission largely reduces the derived age, leading to extremely young (a few Myr) stellar populations in some galaxies. Figure 8 compares the ages derived from models with different nebular emission. The filled gray histograms show the age distributions

from the Lyα-EM models. The black unfilled histograms represent the age distributions from the EM-weak models (top two panels) and the EM-strong models (middle two panels). The distributions of the different models look similar. They all appear to show bimodal distributions. This bimodality could be real, but it could also be caused by selection effects (or sample bias) and modeling limitations. Our galaxies were selected to have strong Lyα emission, which is biased toward younger populations. They were further limited by the IRAC 1 detections, which is biased toward higher stellar masses (older populations) and/or those with strong nebular lines. A full exploration requires a complete, mass-limited sample in this redshift range. The other reason for the bimodality is modeling limitations. The minimum age allowed by our models is 4 Myr, so galaxies younger than 4 Myr would have a measured age of 4 Myr. In addition, these measured “young” ages usually have large uncertainties (typically 20 ~ 30 Myr), as shown in the bottom panel of Figure 8. Therefore, the actual age distribution could be more smooth.

Despite the selection effects and the measurement uncertainties of the ages, we may divide our sample into two subsamples, based on Figures 4 and 8. One consists of old galaxies with ages of several hundred Myr, and the other one includes young galaxies that are mostly only several Myr old. They are referred to as “old” and “young” subsamples in Section 4.1. The bottom panel in Figure 8 shows the relation
<figure><img src="image_4.png" /><figcaption>**Figure 4.** Distributions of stellar mass, age, and $E(B-V)$ derived from the Ly$\alpha$-EM models. Like Figure 1, the gray filled histograms show the results from the Ly$\alpha$-EM models, and the black histograms show the results from the NoEM models as references.</figcaption></figure>

between age and stellar mass derived from the Ly$\alpha$-EM models. It suggests that slightly more than half of the galaxies are a few hundred Myr old. In particular, ages in some galaxies from the rSFH models are older than 300–500 Myr. These galaxies usually have masses close to or higher than $10^{10} M_{\odot}$ (up to $10^{11} M_{\odot}$). These massive and old galaxies already appear to exist when the universe was only 0.8–1.0 Gyr old.

Meanwhile, Figure 8 also shows the existence of extremely young galaxies in our sample. Both dSFH and rSFH models suggest that some galaxies are only several Myr old. These galaxies are usually less massive, with masses between $10^8$ and $3 \times 10^9 M_{\odot}$. Pirzkal et al. (2007) found that some LAEs at $4 < z < 5.7$ in the Hubble Ultra Deep Field are very young (a few Myr old) with masses between $10^6$ and $10^8 M_{\odot}$. The young galaxies in our sample are similar to these galaxies in terms of age, but are certainly more massive. This implies that extremely young galaxies at high redshift already have a wide range of stellar masses.

### 4.3. Dust Extinction

The bottom panels of Figures 1 and 4 display the distributions of dust reddening $E(B-V)$. Unlike age and stellar mass, the $E(B-V)$ values derived from different models, including EM models and NoEM models, are consistent. They clearly suggest that the majority of the galaxies in our sample have little or no dust extinction. This is indeed expected. In Paper I, we reported that our galaxies have steep rest-frame UV slopes on average, with a median value of

<figure><img src="image_5.png" /><figcaption>**Figure 5.** Stellar masses derived from models with different nebular emission. In the upper four panels, the filled gray histograms show the stellar mass distributions from the Ly$\alpha$-EM models. The black unfilled histograms represent the distributions from the EM-weak models and the EM-strong models, which do not significantly deviate from those from the Ly$\alpha$-EM models. The bottom panel shows the distributions of the mass ratios for rSFH. The mass ratios of the GALEV-EM to Ly$\alpha$-EM models span a wide range due to the wide distribution of the Ly$\alpha$ EWs in our sample. The median value is close to one. Models with stronger (weaker) nebular emission generally produce stellar populations with lower (higher) masses. On the other hand, the results for most galaxies from three EM models (Ly$\alpha$-EM, EM-weak, and EM-strong) are roughly consistent. About 70%–80% galaxies in our sample have similar stellar masses (within a factor of three) from the three different EM models.</figcaption></figure>

$\beta \sim -2.3$. Since the UV slope is very sensitive to dust extinction, such steep slopes already indicate little dust extinction. The results are also broadly compatible with recent observations (e.g., Walter et al. 2012; Ouchi et al. 2013; Ota et al. 2014; Willott et al. 2015). Several galaxies in our sample have moderate dust extinction with $E(B-V) > 0.1$. These galaxies are relatively massive with masses higher than $10^9 M_{\odot}$.

### 4.4. Mass-SFR Relation

In low-redshift star-forming galaxies, there is a correlation between stellar mass and SFR, in which SFR is higher in more massive galaxies (e.g., Daddi et al. 2007; Elbaz et al. 2007; Noeske et al. 2007). This relation is referred to as the “main sequence” of star-forming galaxies. It evolves with redshift so
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Relations between stellar mass and brightness. The gray pluses represent the galaxies in the NoEM models with rSFH. For the EM models with rSFH, the galaxies are divided into two subsamples, a “young” subsample (age <30 Myr; blue crosses) and an “old” subsample (age >30 Myr; red crosses). The dotted lines are the best linear (in the log space) fits to the data points. Upper panel: stellar mass as a function of apparent magnitude at 3.6 µm. The masses derived from the NoEM models have a good correlation with the 3.6 µm magnitude, because stellar mass is the amplitude of a model spectrum. For the EM models, the relation is also tight in either subsample, due to the reasons explained in Section 4.1. However, many galaxies are significantly below the relation from the NoEM models, because of the presence of strong nebular emission in the IRAC 1 band. Lower panel: stellar mass as a function of absolute magnitude $M_{1500}$ at rest-frame 1500 Å. For the NoEM models, there is a weak correlation in which stellar masses are higher in more UV luminous galaxies. This correlation is more obvious in the two subsamples for the EM models. The correlation may reflect the mass–SFR relation seen at lower redshifts.</figcaption></figure>

that the normalization of the relation at higher redshift is higher (at least for the redshift range $z \le 2–3$). The mass–SFR relation usually applies to massive galaxies with $M_* \ge 10^9 M_\odot$, corresponding to the old subsample in Figures 6 and 8. In addition, starburst galaxies such as ULIRGs and SMGs may reside well above the relation.

Figure 9 shows the mass–SFR relations for our galaxies, derived from the Lyα-EM models with rSFH (upper panel) and dSFH (lower panel). The pluses and crosses represent the young and old subsamples, respectively. The gray symbols indicate the SFRs calculated from the UV continuum (without

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Same as Figure 6, but for dSFH.</figcaption></figure>

correction for dust extinction) in Paper I. We correct for dust extinction using the $E(B – V)$ values in Table 3 based on the Calzetti et al. (2000) law. The results are shown as blue and red symbols. The correction is large only for the most massive galaxies in each subsample, which show moderate dust extinction as seen in Section 4.3. The two subsamples both show a tight mass–SFR relation. In the upper panel for rSFH, the linear Pearson correlation coefficients for the old and young subsamples are 0.95 and 0.99, respectively. The coefficients in the lower panel for dSFH are 0.77 and 0.95. The slopes are close to 1, consistent with the results from simulations (e.g., Finlator et al. 2011). It also agrees with the slopes found in lower-redshift galaxies (e.g., Guo et al. 2013; Finkelstein et al. 2015; Sparre et al. 2015). The slope of ~1 indicates that the specific SFRs (sSFRs) are similar within either subsample. However, the average sSFR of the young subsample is much higher than that of the old subsample, so that the two subsamples are well separated in the mass–SFR diagram. This is because the two subsamples cover the similar range of SFRs, but the old subsample is 10 times more massive on average.

The correlation in either subsample is tight, due to the combination of a few reasons, including the intrinsic mass–SFR relation, the selection effects, and the relatively small parameter space that either subsample occupies. In particular,
<figure><img src="image_8.png" /></figure>

**Figure 8.** Ages derived from models with different nebular emission. The filled gray histograms show the age distributions from the Lyα-EM models. The black unfilled histograms represent the age distributions from the EM-weak models (top two panels) and the EM-strong models (middle two panels). The distributions of the different models look similar. The bottom panel shows the relation between age and stellar mass derived from the Lyα-EM models. The gray horizontal bars indicate the $1\sigma$ uncertainties. Note that we only show the uncertainty in the age direction and choose the best-fit mass; in reality, there is a strong degeneracy between the two. For the purpose of clarity, the symbols for the rSFH models have been shifted by 0.02 along the y-axis. Many galaxies are older than 200–300 Myr, with stellar masses close to, or higher than $10^{10} M_{\odot}$. Meanwhile, a significant fraction of the galaxies are dominated by extremely young stellar populations with ages of several Myr.

the correlation for the “young” subsample is very tight, mainly because of their extremely young ages. The GALEV model outputs start from 4 Myr in steps of 4 Myr. The ages measured in the “young” subsample are mostly 4 and 8 Myr, and thus these “young” galaxies have not evolved considerably.

Figure 10 shows the sSFR as a function of redshift, taken from Madau & Dickinson (2014; and references therein). It illustrates the efficiency of stellar mass growth in galaxies across cosmic time. The sSFR increases rapidly from the local universe to $z \sim 2$, and then flattens (or slightly climbs) toward higher redshifts. We focus on the high-redshift range. The $z \ge 4$ data points in the figure were from two LBG samples of Stark et al. (2013) and González et al. (2014). Our results are displayed as a blue circle (young subsample) and a red circle (old subsample). The horizontal error bars indicate the redshift range, and the vertical error bars indicate the $2\sigma$ range of the object number distribution (i.e., inclusion of 95% of the

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

**Figure 9.** Stellar mass–SFR relations, derived from the Lyα-EM models with rSFH (upper panel) and dSFH (lower panel). The pluses and crosses represent the young and old subsamples, respectively. The gray symbols indicate the SFRs calculated from the UV continuum (without correction for dust extinction) in Paper I. The blue and red symbols indicate the SFRs corrected for dust extinction. The dotted lines are the best linear (in the log space) fits to the data points. The stellar masses are taken from the Lyα-EM models with rSFH. The two subsamples show a tight mass–SFR relation, with slopes close to 1. The tight relation is likely caused by the combination of the intrinsic mass–SFR relation, selection effects, and relatively small parameter space that two subsamples occupy.

objects). The average sSFR of the old subsample agrees well with the two previous studies Stark et al. (2013) and González et al. (2014). Note that in these studies nebular emission lines were also incorporated during their SED modeling. In addition, the mass–SFR relation of the old subsample is consistent with the simulation of Finlator et al. (2011). All of these suggest that the galaxies in our old subsample are “normal” star-forming galaxies at $z \ge 6$. The sSFRs of the old subsample are roughly consistent with (or marginally higher than) those at $z \sim 2$, suggesting that the efficiency of stellar mass growth did not change much during most of the first 3 Gyr.

The average sSFR of the young subsample is about 10 times higher than the relation defined by previous studies. It has been clear that galaxies with strong starburst activity at $z \le 2$ are well above (≥10 times) the main sequence of star-forming galaxies (e.g., Rodighiero et al. 2011; Sargent et al. 2012). So the young subsample seems to be the high-redshift counterparts of $z \le 2$ starbursts. The fraction of such starbursts at $z \le 2$ is very low, while in our sample this fraction is much higher, partly due to selection effects. On the other hand, galaxies with very high sSFRs at high redshift are not rare. Section 5.3 provides more discussion on this topic.

The bimodal distribution of ages and sSFRs seen in Figures 8 and 9 were largely explained by selection effects and modeling
<figure><img src="image_10.png" /></figure>

**Figure 10.** sSFR as a function of redshift. The gray circles represent the data points taken from Madau & Dickinson (2014; and references therein). The open green circles represent some of high-sSFR galaxies at $z > 4$ taken from the literature (Pirzkal et al. 2007; Bowler et al. 2012; Finkelstein et al. 2013; Huang et al. 2015). Our results of the median sSFRs are displayed as a blue circle (young subsample) and a red point (old subsample). The horizontal error bars indicate the redshift range, and the vertical error bars indicate the $2\sigma$ range of the object number distribution (i.e., inclusion of 95% of the objects). The sSFR of the old subsample is consistent with previous studies compiled by Madau & Dickinson (2014), following the main sequence of star-forming galaxies. The sSFR of the young subsample is well above the main sequence, presumably due to starburst activity in these galaxies. Such galaxies are not rare at high redshift, as seen from the open circles.

limitations in the above sections, though we cannot rule out the
possibility of an intrinsic bimodality. Such bimodality has been
reported for galaxies at $2.5 < z < 3.5$ by Kajisawa et al.
(2010), who found that the sSFRs in their low-sSFR and high-
sSFR galaxies are $0.5-1.0 \text{ Gyr}^{-1}$ and $\sim 10 \text{ Gyr}^{-1}$, respectively.
This suggests that the bimodality seen in our sample may
partially reflect a real bimodal distribution of sSFRs.

# 5. DISCUSSION

## 5.1. *Testing Lyα-EM Models*

It is clear that SED modeling of high-redshift galaxies can be largely affected by the presence of strong nebular emission. For galaxies at $6 \le z < 7$, the existence of strong lines is often evidenced by their IRAC 1 flux excess (compared to the IRAC 2 flux). This is because in this redshift range, the IRAC 1 band covers some of the strongest lines such as [O III], $H\alpha$, and $H\beta$. A significant fraction of galaxies at $z \ge 6$, including photometrically selected and spectroscopically confirmed galaxies, show a strong IRAC 1 flux excess (e.g., González et al. 2010; McLure et al. 2011; Curtis-Lake et al. 2013). In fact, many galaxies in our sample also show such an excess. At $z > 7$, galaxies with strong nebular lines start to show red IRAC [3.6]–[4.5] colors (e.g., Roberts-Borsani et al. 2015; Zitrin et al. 2015).

The contribution of nebular emission to the IRAC 1 and 2
bands can be large or dominant. Stark et al. (2013) showed that
the mean rest-frame EW of $H\alpha$ at high redshift is a few
hundred Å. Smit et al. (2015) showed some extreme cases of
strong emission lines by searching for $z = 6.6–6.9$ galaxies
with very blue IRAC [3.6]–[4.5] colors. Over a small redshift
range of $z = 6.6–6.9$, the IRAC 1 band covers [O III] and $H\beta$,

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

**Figure 11.** Distribution of the scaling factors of nebular line strength derived for the Lyα-EM models. The Lyα-EM models were computed from the GALEV-EM models by scaling nebular emission lines to match the observed Lyα flux (Section 3). The derived scaling factors cover a range from 0.6 to 3.5, with a median value of 1.6.

but the IRAC 2 band does not cover any strong emission lines. The galaxies that they found have very large rest-frame EW of [O III] + $H\beta$ in the range of 900 to >2000 Å, with a median value of ~1400 Å. In these galaxies, line emission dominates the IRAC 1 photometry. Therefore, it is expected that the galaxies in our sample have strong nebular lines.

Our analysis was mostly based on the Lyα-EM models. When we computed the Lyα-EM models from the GALEV-EM models in Section 3, we scaled nebular emission lines using the observed Lyα flux. Figure 11 shows the distribution of the derived scaling factors. They cover a range from 0.6 to 3.5, with a median value of 1.6. It is expected that line emission derived from the Lyα-EM models is stronger than that from the GALEV-EM models, because our galaxies were selected to have strong line emission. The scaling factors span a relatively small range (within a factor of ~2 around the median value of 1.6), suggesting that the dynamic range of the EM-strong and EM-weak models is large enough to sample the majority of high-redshift galaxies.

We may test our Lyα-EM models using a few “special”
galaxies (Figure 12) not in our sample. These galaxies are
spectroscopically confirmed with measured Lyα line flux, so
that we can estimate nebular emission from their Lyα emission
based on the Lyα-EM models. In addition, they are at certain
redshift ranges so that one of the two IRAC bands covers
strong nebular emission, but the other one does not. In this
case, the difference of the photometry between the two bands
roughly reflects the strength of nebular lines. The first two
galaxies (a and b in Figure 12) are at $z = 6.74$ and $z = 6.76$
found by B. Clément et al. (2016, in preparation) and by Huang
et al. (2015), respectively. In these two galaxies, the IRAC 1
band covers $H\beta$ and [O III], while the IRAC 2 band does not
cover any strong lines. The other two galaxies in Figure 12 are
the one at $z = 7.51$ from Finkelstein et al. (2013) and the one at
$z = 7.73$ from Oesch et al. (2015). For these two galaxies, the
IRAC 2 band covers $H\beta$ and [O III], and the IRAC 1 band does
not cover strong lines.
<figure><img src="image_12.png" /><figcaption>**Figure 12.** SED modeling of four galaxies taken from the literature to test our Lyα-EM models. These galaxies are spectroscopically confirmed with measured Lyα line flux. They are at certain redshift ranges so that one of the two IRAC bands covers strong nebular emission, but the other one does not. We also show the difference between the observed and model photometry in the IRAC bands ($\Delta m_{IRAC1}$ and $\Delta m_{IRAC2}$). The fitting results for all of the galaxies except (c) are good, and their magnitude difference values are smaller than the corresponding photometric uncertainties ($1\sigma$) in the IRAC 1 and IRAC 2 bands. The best fit to galaxy (c) is not acceptable ($\chi_r^2 \gg 1$) because the nebular lines (scaled from Lyα) in the model spectra are not strong enough to account for the large difference between the IRAC 1 and IRAC 2 photometry. This is because its Lyα emission has been largely attenuated by the highly neutral IGM at very high redshift.</figcaption></figure>

We perform SED modeling for the four galaxies based on the Lyα-EM models with rSFH. The fitting results for galaxies (a) and (d) are very good ($\chi_r^2 \le 1$). We are not able to obtain acceptable results for galaxies (b) and (c), due to the reason that the nebular lines (scaled from Lyα) in the model spectra are not strong enough to account for the large difference between the IRAC 1 and IRAC 2 photometry. We note that the Lyα emission line in galaxy (b) is located on one of the strong sky OH lines, and its flux measurement could be significantly affected, as pointed out by the discovery paper (Huang et al. 2015). Nevertheless, we re-model the SEDs for galaxies (b) and (c) using the EM-strong models, and obtain a reasonable fit for galaxy (b). For galaxy (c), however, we fail to achieve an acceptable fit ($\chi_r^2 \gg 1$). The best fitting results for the four galaxies are plotted in Figure 12, in which we also show the difference between the observed and model photometry in the IRAC bands ($\Delta m_{IRAC1}$ and $\Delta m_{IRAC2}$). For all of the galaxies except (c), these magnitude difference values are smaller than the corresponding photometric uncertainties ($1\sigma$) in the IRAC 1 and IRAC 2 bands.

Galaxy (c) has very weak Lyα emission (compared to its other nebular lines), as already noted by Finkelstein et al. (2013). At $z = 7.51$, the IGM is much more neutral than that at $z \sim 6$, so Lyα emission can be largely attenuated by the neutral IGM or eaten by Lyα damping wings (e.g., Miralda-Escudé 1998). In other words, the ratio of intrinsic to observed Lyα flux is a function of redshift, as we mentioned in the previous section. So our assumption about this ratio is no longer valid. Recent simulations suggest that the Lyα damping

<figure><img src="image_13.png" /><figcaption>**Figure 13.** Physical properties of the “young” and “old” subsamples as a function of rest-frame UV luminosity $M_{1500}$. The blue and red circles represent the “young” and “old” galaxies, respectively. The gray circles (including filled and open circles) represent the galaxies that are not included for analysis in this paper. The top panel shows the rest-frame UV slope $\beta$ derived in Paper I. The two subsamples have similar mean slopes. The middle panel shows the rest-frame Lyα EW derived in Paper I. The “young” galaxies have relatively higher Lyα EWs than the “old” ones. The bottom panel shows the half-light radii $r_{hl,in}$ estimated in Paper II. The average radii of the two subsamples are roughly consistent.</figcaption></figure>

wing owing to patchy reionization should be fairly uniform at a given redshift (Mesinger et al. 2015). In this case, our ability to fit galaxy (d) but not (c) is not likely to reflect incomplete reionization. Instead, it may indicate scatter in the level of attenuation of self-shielded systems or in the intrinsic properties of the galaxies’ interstellar media.

In summary, three out of the four galaxies can be well fit with our Lyα-EM or EM-strong models. We were not able to obtain a reasonable fit for the galaxy at $z = 7.51$, due to the limitation (redshift coverage) of our models. The tests above show that our models can provide reasonable SED modeling for high-redshift galaxies with nebular emission taken into account.

## 5.2. “Young” and “Old” Populations

In Section 4, we identified two subsamples in our galaxies, one “old” subsample with ages of several hundred Myr and one “young” subsample with ages of several Myr. Either subsample follows a tight stellar mass–SFR relation. The mean sSFR of
the “old” subsample is consistent with those in many high-redshift galaxies reported in the literature, while the mean sSFR of the “young” subsample is an order of magnitude higher. In Figure 13, we compare the two subsamples in the context of physical properties as a function of rest-frame UV luminosity $M_{1500}$. These physical quantities (including rest-frame UV slopes, Ly$\alpha$ EWs, and half-light radii) were measured in Papers I and II. Their measurements are usually associated with considerable uncertainties (see Papers I and II), which are not plotted in the figure for the purpose of simplicity. The blue and red circles represent the “young” and “old” galaxies, respectively. The gray circles (including filled and open circles) represent the galaxies that are not included for analysis in this paper.

The top panel of Figure 13 shows the rest-frame UV slope $\beta$ as a function of $M_{1500}$. The median and standard deviation values of the slopes for the young and old subsamples are $-2.31 \pm 0.40$ and $-2.20 \pm 0.62$, respectively. They are consistent. UV slopes are most sensitive to dust extinction. As we have seen in Section 4.3, there is little dust extinction in these galaxies except for several of the most massive galaxies. So the “young” galaxies do not show bluer slopes than the “old” galaxies.

The middle panel of Figure 13 shows the rest-frame Ly$\alpha$ EW as a function of $M_{1500}$. The median and standard deviation values of the EWs for the two subsamples are $54.2 \pm 30.7$ and $39.3 \pm 29.5$, respectively. The “young” galaxies have relatively higher Ly$\alpha$ EWs than the “old” ones. A Ly$\alpha$ EW is the ratio of the Ly$\alpha$ line flux to the continuum flux. The Ly$\alpha$ line strength measures the strength of nebular lines, especially in our Ly$\alpha$-EM models. We expect to see stronger nebular lines in younger systems. So this panel reflects that a galaxy with stronger Ly$\alpha$ line emission tend to be younger.

The bottom panel shows the half-light radii $r_{hl,in}$ (at rest-frame $\sim 1800$ Å) as a function of $M_{1500}$. The radii have been corrected for PSF broadening with simulations. The median and standard deviation values of the radii for the two subsamples are $0.15 \pm 0.06$ and $0.20 \pm 0.06$, respectively. The old galaxies are marginally larger on average. However, it is not straightforward to compare the sizes of the two subsamples. At lower redshift, there is a mass-size relation, in which more massive galaxies (the “old” subsample in our case) tend to be larger. On the other hand, the galaxies in our “young” subsample may have strong starburst activity as indicated by their high sSFRs. A significant fraction of low-redshift starbursts are mergers, which tend to have large sizes.

### 5.3. Comparison with Previous Studies

There are a number of studies on the stellar populations of $z \ge 6$ galaxies in the literature. The majority of these galaxies are photometrically selected LBGs. There was little study for $z \ge 6$ LAEs. As we explained in the Introduction, almost all of the known LAEs were discovered by ground-based telescopes, and did not have deep infrared observations. Our program provides the largest sample so far for studying stellar populations in spectroscopically confirmed LAEs at $z \ge 6$.

Direct comparison with previous studies is very difficult if not impossible, because different studies use different galaxy samples (selected from different data sets with different selection criteria) and SED modeling methods. We mainly focus on two questions: whether extremely young populations have been reported, and whether previously reported sSFRs

cover a wide range that agrees with our “young” and “old” subsamples. As we emphasized earlier, age is poorly constrained from SED modeling, so here we broadly define “extremely young galaxies” as those with the best-fitting ages younger than $\sim 30$ Myr, like the galaxies in our “young” subsample.

In the previous studies of photometrically selected LBGs, extremely young galaxies have rarely been reported. The majority of these LBGs tend to have relatively old and mature populations with ages of a few hundred or several 10 of Myr, from SED fitting without nebular emission taken into account. When nebular emission was considered (especially in models with smoothly varying SFHs), however, extremely young populations were required to explain the SEDs of some LBGs (Schaerer & de Barros 2010; McLure et al. 2011; Curtis-Lake et al. 2012; Huang et al. 2015). For example, McLure et al. (2011) selected a sample of $z > 6$ LBGs from HST deep fields, and found that these galaxies were mostly a few hundred Myr old if nebular emission was not included in their SED modeling. When nebular emission was added, about 25% of their LBGs with IRAC 1 detections were found to be extremely young, with a median age of 9 Myr.

In the previous studies of high-redshift LAEs, extremely young galaxies were found to be common. For example, Ono et al. (2010) stacked a large number of photometrically selected LAEs at $z \sim 5.7$ and 6.5, and performed SED modeling on the stacked LAEs. Despite the fact that the measurements from stacked data could be unreliable (e.g., Vargas et al. 2014), they found that these LAEs have ages of only 1–3 Myr, little dust extinction, and strong nebular emission. In the sample of Pirzkal et al. (2007), there are nine spectroscopically confirmed LAEs at $z \sim 5$. Most of them were found to be younger than 10 Myr. These results are quite consistent with ours.

The sSFRs in our galaxies span quite a large range. Their values are about $3-4$ Gyr$^{-1}$ in the “old” subsample, which is consistent with previous studies such as Curtis-Lake et al. (2012), Stark et al. (2013), and González et al. (2014). The sSFRs in our “young” subsample are about an order of magnitude higher, likely due to starburst activity in these galaxies. Such high sSFRs are indeed also common in previous studies. For example, in a sample of seven secure high-redshift LBGs by Bowler et al. (2012), three of them were found to have sSFRs around $4-5$ Gyr$^{-1}$, and another three have sSFRs close to or higher than $30$ Gyr$^{-1}$. These numbers roughly agree with those found in our “old” and “young” subsamples. Some studies have reported even stronger starbursts in $z > 6$ galaxies, with sSFRs between one hundred and several hundred Gyr$^{-1}$ (e.g., Ono et al. 2010; Finkelstein et al. 2013; Huang et al. 2015). In Figure 10, the green circles represent some of the high-sSFR galaxies from Pirzkal et al. (2007), Bowler et al. (2012), Finkelstein et al. (2013), and Huang et al. (2015). Galaxies ($4 < z < 6$) from Pirzkal et al. (2007) all have very high sSFRs, and are included in Figure 10. All of the above suggests that strong starburst activity is common in very-high-redshift galaxies.

### 5.4. LAEs and LBGs

In this series of papers (including this paper and Papers I and II), LAEs are defined as galaxies found by the narrow-band (or Ly$\alpha$) technique, and LBGs are defined as galaxies found by the dropout technique. As we already pointed out in Papers I and II, this widely used classification only reflects the
<figure><img src="image_14.png" /><figcaption>**Figure 14.** Measured ages and SFRs as a function of stellar mass for the LAEs and LBGs in our sample. The minimum age allowed in our models is 4 Myr, and some ages in the upper panel are at this limit. The distributions of the LBGs and LAEs in the age–mass and SFR–mass diagrams are indistinguishable. The number ratios of the LAEs to LBGs agree with each other in the “young” and “old” subsamples.</figcaption></figure>

methodology that we apply to select galaxies. It does not mean that the two types of galaxies are intrinsically different. Another definition of LAEs is based on the Lyα EW, e.g., a galaxy is an LAE if its Lyα EW is greater than 20 Å. This definition is physically more meaningful, but observationally difficult, because one can easily obtain a flux-limited sample, not an EW-limited sample. In addition, it is meaningless for broadband selected galaxies. Thus, we use the former definition.

It is not entirely clear whether high-redshift LAEs and LBGs represent physically different populations. Direct comparison between LAEs and LBGs is difficult because of the very different target selection procedures. As we already explained in Papers I and II, the LAEs in our original sample of 67 galaxies are composed of a well-defined sample in terms of Lyα flux. On the other hand, the LBGs in the sample only represent the LBGs with strong Lyα emission since they are spectroscopically confirmed. In Papers I and II, we compared the LAEs and LBGs in our sample in great detail, and found that the two populations are indistinguishable in all aspects of physical properties that we considered, including the Lyα emission strength, UV continuum properties, sizes, and morphology, etc.

Figure 14 shows the distributions of the measured stellar masses, ages, and SFRs for the LAEs and LBGs in our sample

of 27 galaxies used in this paper. The blue pluses and the red crosses represent the LBGs and LAEs, respectively. The distributions of the LBGs and LAEs in the age–mass and SFR–mass diagrams are indistinguishable. In particular, there are 6 LAEs and 5 LBGs in the “young” subsample, and there are 9 LAEs and 7 LBGs in the “old” subsample. The number ratios of the LAEs to LBGs nicely agree with each other in the two subsamples. All of these are consistent with our previous conclusion that the LAEs and LBGs in our sample have common properties, suggesting that LAEs are a subset of LBGs with strong Lyα emission lines. The conclusion is also consistent with recent models and simulations (e.g., Dayal & Ferrara 2012; Dijkstra & Wyithe 2012; Garel et al. 2015).

5.5. Prospect for Spectroscopic Follow-up with the James Webb Space Telescope (JWST)

Currently there are two ways to assess the strength of rest-frame optical nebular emission lines and use that information to derive the physical properties of underlying stellar population in *z* > 6 galaxies: (1) to observe galaxies in specific redshift ranges where one of the IRAC band is line-free (e.g., 6.6 < *z* < 6.9 galaxies as studied by Smit et al. 2015), and (2) to model emission-line strengths based on the measured Lyα lines (this work). The first method is more direct and accurate, but because of the stringent constraints on redshift, the number of galaxies that can be studied is rather small. The second method suffers from relatively large uncertainties, but can be applied to a much larger sample of *z* > 6 galaxies as we have shown in this paper.

Although the above is the best we can do at the moment, we note that we are actually on the verge of a breakthrough as the JWST becomes available and changes the emphasis from photometry to spectroscopy. Figure 15 shows the flux distributions of the [O III] 5007 Å and Hα lines predicted by the modeling described in this paper. The figure shows that we should be able to detect almost all of these lines with a line-flux sensitivity of ~10<sup>-17</sup> erg cm<sup>-2</sup> s<sup>-1</sup>. Based on the JWST Prototype Exposure Time Calculator (ETC), such a sensitivity will be easily reachable with NIRSpec. For example, to achieve a 5σ detection (per resolution element) of an [O III] 5007 Å line at z = 6.5 with a line flux of 10<sup>-17</sup> erg cm<sup>-2</sup> s<sup>-1</sup>, the required integration time will be only ~500 s. The corresponding integration time for Hα (i.e., at z = 6.5 and with a line flux of 10<sup>-17</sup> erg cm<sup>-2</sup> s<sup>-1</sup>) will be ~800 s. These numbers clearly indicate that with JWST/NIRSpec, it will be easy to detect [O III] and Hα lines with the brightest *z* > 6 galaxies like those we have studied here. In fact, NIRSpec has the sensitivity to detect much fainter lines, making these bright *z* > 6 SDF galaxies ideal targets for a detailed NIRSpec spectroscopic study. In several years, such a study will be able to test the validity of the results presented here, and will undoubtedly make many more interesting discoveries.

6. SUMMARY

This paper is the third in a series presenting the physical properties of a large sample of 67 spectroscopically confirmed galaxies at z ≥ 6. The sample consists of 51 LAEs at z ≃ 5.7, 6.5, and 7.0, and 16 LBGs at 5.9 ≤ z ≤ 6.5. They have deep optical imaging data from Subaru, near-IR data from HST, and mid-IR data from Spitzer. In this paper, we have reported a detailed study of stellar populations in these galaxies. We have
<figure><img src="image_15.png" /><figcaption>**Figure 15.** Flux distributions of the [O III] 5007 Å and Hα lines predicted by the modeling described in this paper. The flux was calculated based on the JWST prototype ETC (version P1.6), available at http://jwstetc.stsci.edu. The following parameters were used with the NIRSpec ETC: Filter/Grating = G395H/F290LP; R = 2700; MSA shutter=0′′2 × 0′′45; Flat continuum in Fν, FWHM = 40 Å ~ 300 km s⁻¹ at 3.755 µm for [O III], or FWHM = 50 Å ~ 300 km s⁻¹ at 4.922 µm for Hα. The continuum normalization was set to be 30 AB mag at 3.6 µm, so that the S/N calculation is essentially driven by the line flux. With JWST/NIRSpec, it will be easy to detect [O III] and Hα lines in bright z > 6 galaxies.</figcaption></figure>

focused on a subsample of 27 galaxies with Spitzer IRAC 1 detections at 3.6 µm. This subsample represents luminous and massive galaxies with strong Lyα emission at z ≥ 6. We used the wealth of the multi-band data and the secure Lyα redshifts and flux to model the SEDs of the 27 galaxies and characterize their stellar populations.

We used the GALEV evolutionary synthesis models with nebular continuum and line emission to mainly constrain three physical parameters: age, stellar mass, and dust extinction. We adopted two representative SFHs, an exponentially dSFH and a smoothly rSFH. We mostly used the latter two SFHs in our analysis. In order to incorporate nebular emission, we scaled the nebular emission from the GALEV models to match the observed Lyα flux (the Lyα-EM models) under two simple assumptions. With the Lyα-EM models, we were able to nicely break the strong degeneracy of model spectra between young galaxies with prominent nebular emission and mature galaxies with strong Balmer breaks.

Our best-fitting results show that the galaxies in our sample have a wide range of SED ages from several Myr to a few hundred Myr. They also have a wide range of stellar masses from ~10⁸ to ~10¹¹ M☉. Interestingly, the distribution of the measured ages (despite the large uncertainties) appear to be bimodal, likely due to selection effects and modeling limitations (though we cannot rule out the possibility of an intrinsic bimodality). Based on this bimodality, we divided the galaxies into two subsamples: an “old” subsample and a

“young” subsample. The “old” subsample mainly consists of galaxies older than 100 Myr, with stellar masses higher than 10⁹ M☉. Many galaxies are older than 300–500 Myr. The “young” subsample consists of galaxies younger than 30 Myr (usually several Myr old). These galaxies are less massive, with masses ranging between ~10⁸ and ~3 × 10⁹ M☉. The majority of the galaxies in both subsamples show little or no dust extinction, as already hinted by their steep rest-frame UV slopes.

Both subsamples show a correlation between stellar mass and SFR, but with very different normalizations. The mean sSFR of the “old” subsample is about 3–4 Gyr⁻¹, consistent with the mass–SFR relation defined by previous studies. The mean sSFR of the “young” subsample is an order of magnitude higher. Such higher sSFRs have also been frequently reported in previous studies. They are likely due to starburst activity in these galaxies. Finally, the LAEs and LBGs in our sample are indistinguishable in all physical properties that we have considered, suggesting that LAEs are a subset of LBGs with strong Lyα emission lines.

We acknowledge the support from a 985 project at Peking University. L.J., S.C., and E.E. also acknowledge the support from NASA through awards issued by STScI (HST PID: 11149,12329,12616) and by JPL/Caltech (Spitzer PID: 40026,70094). We would like to thank R. Ryan for his advice on IRAC photometry and using iGALFIT. We also thank J. Rhoads and Z. Zheng for their insightful comments, and thank P. Madau and P. A. Oesch for providing us with the data used for Figures 10 and 12.

*Facilities: HST (NICMOS,WFC3), Spitzer (IRAC), Subaru (Suprime-Cam).*

## REFERENCES

Ashby, M. L. N., Willner, S. P., Fazio, G. G., et al. 2013, ApJ, 769, 80

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2014, ApJ, 793, 115

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34

Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2012, MNRAS, 426, 2772

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

Cai, Z.-Y., Lapi, A., Bressan, A., et al. 2014, ApJ, 785, 65

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Conroy, C. 2013, ARA&A, 51, 393

Curtis-Lake, E., McLure, R. J., Dunlop, J. S., et al. 2013, MNRAS, 429, 302

Curtis-Lake, E., McLure, R. J., Dunlop, J. S., et al. 2014, arXiv:1409.1832

Curtis-Lake, E., McLure, R. J., Pearce, H. J., et al. 2012, MNRAS, 422, 1425

Daddi, E., Dickinson, M., Morrison, G., et al. 2007, ApJ, 670, 156

Dayal, P., & Ferrara, A. 2012, MNRAS, 421, 2568

Dijkstra, M. 2014, PASA, 31, 40

Dijkstra, M., & Wyithe, J. S. B. 2012, MNRAS, 419, 3181

de Barros, S., Schäerer, D., & Stark, D. P. 2014, A&A, 563, A81

Dunlop, J. S., Rogers, A. B., McLure, R. J., et al. 2013, MNRAS, 432, 3520

Egami, E., Kneib, J.-P., Rieke, G. H., et al. 2005, ApJL, 618, L5

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

Ellis, R. S., McLure, R. J., Dunlop, J. S., et al. 2013, ApJL, 763, L7

Eyles, L. P., Bunker, A. J., Ellis, R. S., et al. 2007, MNRAS, 374, 910

Eyles, L. P., Bunker, A. J., Stanway, E. R., et al. 2005, MNRAS, 364, 443

Fan, X., Narayanan, V. K., Lupton, R. H., et al. 2001, AJ, 122, 2833

Finkelstein, K. D., Finkelstein, S. L., Tilvi, V., et al. 2015, ApJ, 813, 78

Finkelstein, S. L., Papovich, C., Dickinson, M., et al. 2013, Natur, 502, 524

Finlator, K., Davé, R., & Oppenheimer, B. D. 2007, MNRAS, 376, 1861

Finlator, K., Oppenheimer, B. D., & Davé, R. 2011, MNRAS, 410, 1703

Furusawa, H., Kosugi, G., Akiyama, M., et al. 2008, ApJS, 176, 1

Garel, T., Blaizot, J., Guiderdoni, B., et al. 2015, MNRAS, 450, 1279

González, V., Bouwens, R., Illingworth, G., et al. 2014, ApJ, 781, 34

González, V., Labbé, I., Bouwens, R. J., et al. 2010, ApJ, 713, 115

González, V., Labbé, I., Bouwens, R. J., et al. 2011, ApJL, 735, L34

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Guo, K., Zheng, X. Z., & Fu, H. 2013, ApJ, 778, 23
Henry, A. L., Martin, C. L., Dressler, A., Sawicki, M., & McCarthy, P. 2012, ApJ, 744, 149

Hu, E. M., Cowie, L. L., Barger, A. J., et al. 2010, ApJ, 725, 394

Huang, K.-H., Bradač, M., Lemaux, B. C., et al. 2015, arXiv:1504.02099

Iye, M., Ota, K., Kashikawa, N., et al. 2006, Natur, 443, 186

Jensen, H., Hayes, M., Iliev, I. T., et al. 2014, MNRAS, 444, 2114

Jiang, L., Egami, E., Fan, X., et al. 2013a, ApJ, 773, 153

Jiang, L., Egami, E., Kashikawa, N., et al. 2011, ApJ, 743, 65

Jiang, L., Egami, E., Mechtley, M., et al. 2013b, ApJ, 772, 99

Jiang, L., Fan, X., Annis, J., et al. 2008, AJ, 135, 1057

Labbé, I., Oesch, P. A., Bouwens, R. J., et al. 2013, ApJL, 777, L19

Laporte, N., Streblyanska, A., Kim, S., et al. 2015, A&A, 575, A92

Kajisawa, M., Ichikawa, T., Yamada, T., et al. 2010, ApJ, 723, 129

Kashikawa, N., Shimasaku, K., Matsuda, Y., et al. 2011, ApJ, 734, 119

Kashikawa, N., Shimasaku, K., Yasuda, N., et al. 2004, PASJ, 56, 1011

Kawamata, R., Ishigaki, M., Shimasaku, K., Oguri, M., & Ouchi, M. 2015, ApJ, 804, 103

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Konno, A., Ouchi, M., Ono, Y., et al. 2014, ApJ, 797, 16

Kotulla, R., Fritze, U., Weilbacher, P., & Anders, P. 2009, MNRAS, 396, 462

Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Malhotra, S., & Rhoads, J. E. 2002, ApJL, 565, L71

Matthee, J., Sobral, D., Santos, S., et al. 2015, MNRAS, 451, 400

McLure, R. J., Dunlop, J. S., de Ravel, L., et al. 2011, MNRAS, 418, 2074

Mesinger, A., Aykutalp, A., Vanzella, E., et al. 2015, MNRAS, 446, 566

Miralda-Escudé, J. 1998, ApJ, 501, 15

Mobasher, B., Dahlen, T., Ferguson, H. C., et al. 2015, ApJ, 808, 101

Momose, R., Ouchi, M., Nakajima, K., et al. 2014, MNRAS, 442, 110

Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJL, 660, L43

Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2010, ApJL, 709, L16

Oesch, P. A., van Dokkum, P. G., Illingworth, G. D., et al. 2015, ApJL, 804, L30

Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713

Ono, Y., Ouchi, M., Mobasher, B., et al. 2012, ApJ, 744, 83

Ono, Y., Ouchi, M., Shimasaku, K., et al. 2010, ApJ, 724, 1524

Ota, K., Walter, F., Ohta, K., et al. 2014, ApJ, 792, 34

Ouchi, M., Ellis, R., Ono, Y., et al. 2013, ApJ, 778, 102

Ouchi, M., Shimasaku, K., Furusawa, H., et al. 2010, ApJ, 723, 869

Papovich, C., Dickinson, M., & Ferguson, H. C. 2001, ApJ, 559, 620

Papovich, C., Finkelstein, S. L., Ferguson, H. C., Lotz, J. M., & Giavalisco, M. 2011, MNRAS, 412, 1123

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266

Pentericci, L., Vanzella, E., Fontana, A., et al. 2014, ApJ, 793, 113

Pirzkal, N., Malhotra, S., Rhoads, J. E., & Xu, C. 2007, ApJ, 667, 49

Pirzkal, N., Rothberg, B., Nilsson, K. K., et al. 2012, ApJ, 748, 122

Rhoads, J. E., Hibon, P., Malhotra, S., Cooper, M., & Weiner, B. 2012, ApJL, 752, L28

Roberts-Borsani, G. W., Bouwens, R. J., Oesch, P. A., et al. 2015, arXiv:1506.00854

Robertson, B. E., Ellis, R. S., Dunlop, J. S., McLure, R. J., & Stark, D. P. 2010, Natur, 468, 49

Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJL, 739, L40

Ryan, R. E., Jr 2011, arXiv:1110.1090

Sargent, M. T., Béthermin, M., Daddi, E., & Elbaz, D. 2012, ApJL, 747, L31

Schaerer, D., & de Barros, S. 2009, A&A, 502, 423

Schaerer, D., & de Barros, S. 2010, A&A, 515, A73

Schenker, M. A., Stark, D. P., Ellis, R. S., et al. 2012, ApJ, 744, 179

Shapley, A. E., Steidel, C. C., Erb, D. K., et al. 2005, ApJ, 626, 698

Shibuya, T., Kashikawa, N., Ota, K., et al. 2012, ApJ, 752, 114

Shimasaku, K., Kashikawa, N., Doi, M., et al. 2006, PASJ, 58, 313

Silva, M. B., Santos, M. G., Gong, Y., Cooray, A., & Bock, J. 2013, ApJ, 763, 132

Smit, R., Bouwens, R. J., Franx, M., et al. 2015, ApJ, 801, 122

Sobral, D., Matthee, J., Darvish, B., et al. 2015, ApJ, 808, 139

Sparre, M., Hayward, C. C., Springel, V., et al. 2015, MNRAS, 447, 3548

Stark, D. P., Ellis, R. S., Bunker, A., et al. 2009, ApJ, 697, 1493

Stark, D. P., Schenker, M. A., Ellis, R., et al. 2013, ApJ, 763, 129

Treu, T., Schmidt, K. B., Trenti, M., Bradley, L. D., & Stiavelli, M. 2013, ApJ, 775, 29

Vargas, C. J., Bish, H., Acquaviva, V., et al. 2014, ApJ, 783, 26

Walter, F., Decarli, R., Carilli, C., et al. 2012, ApJ, 752, 93

Watson, D., Christensen, L., Kraiberg Knudsen, K., et al. 2015, Natur, 519, 327

Willott, C. J., Carilli, C. L., Wagg, J., & Wang, R. 2015, ApJ, 807, 180

Willott, C. J., McLure, R. J., Hibon, P., et al. 2013, AJ, 145, 4

Windhorst, R. A., Cohen, S. H., Hathi, N. P., et al. 2011, ApJS, 193, 27

Yan, H., Dickinson, M., Giavalisco, M., et al. 2006, ApJ, 651, 24

Yan, H., Dickinson, M., Stern, D., et al. 2005, ApJ, 634, 109

Yan, H., Finkelstein, S. L., Huang, K.-H., et al. 2012, ApJ, 761, 177

Zheng, Z., Cen, R., Trac, H., & Miralda-Escudé, J. 2010, ApJ, 716, 574

Zitrin, A., Labbé, I., Belli, S., et al. 2015, ApJL, 810, L12