SPITZER ULTRA FAINT SURVEY PROGRAM (SURFS UP). II. IRAC-DETECTED LYMAN-BREAK GALAXIES AT 6 ≲ z ≲ 10 BEHIND STRONG-LENSING CLUSTERS

We list the eight galaxy clusters analyzed in this work in Table 1. Among the eight clusters, six (MACS 0717, MACS 0744, MACS 1149, RXJ 1347, MACS 1423, and MACS 2129) are in the Cluster Lensing And Supernova Survey (CLASH; Postman et al. 2012) sample; therefore, each of them has HST imaging data in at least twelve ACS/WFC and WFC3/IR filters¹⁵ from the ACS and WFC3/IR cameras. The typical 5σ depths for the CLASH clusters reported by Postman et al. (2012) are between 27.0 and 27.5 mag (within 04 diameter apertures) for each filter, and the large number of filters provides unique constraints for high-z galaxy searches among the HST deep fields.

In addition to the CLASH data, full-depth images for MACS 0717 and MACS 1149 from the HFF program were also released in 2015 June. In these two clusters, HST spent a total ∼140 orbits that are roughly split between the ACS and WFC3/IR filters, and these images achieve ≈28.7–29 mag¹⁶ in the optical (ACS) and NIR (WFC3). We use the deepest HFF images for MACS 1149 and MACS 0717 for photometry in the filters where such images are available.¹⁷ The six CLASH clusters in our sample are also observed by the Grism Lens-Amplified Survey from Space (GLASS; PI: Treu; Schmidt et al. 2014; Treu et al. 2015) program (MACS 0717, MACS 1149, MACS 1423, RXJ 1347, MACS 0744, and MACS 2129) and have deep grism spectra available.

For the remaining two clusters being analyzed in this work, we have data for RCS 2327 as part of the SURFS UP HST observations (HST-GO13177 PI Bradač; Hoag et al. 2015) and previous archival data (HST-GO10846 PI Gladders; see also Sharon et al. 2015). For MACS 0454, we use the archival observations from HST-GO11591 (PI: Kneib), GO-9836 (PI: Ellis), and GO-9722 (PI: Ebeling). We list the limiting magnitudes of point sources within a 04 aperture for MACS 0454 and RCS 2327 in Table 2.

We will use the template-fitting software T-PHOT (Merlin et al. 2015)—the successor to TFIT (Laidler et al. 2007)—to measure the colors between HST and Spitzer IRAC images (see Section 2.2). To prepare the HST images for T-PHOT, we use the public 003 pixel⁻¹ CLASH images and block-sum the images to make 006 pixel⁻¹ images. We also edit the astrometric image header values (CRVALs and CRPIXs) to conform to T-PHOT's astrometric requirements¹⁸ and make sure that HST and Spitzer images are aligned to well within 01 (Paper I).

We also match the point-spread functions (PSFs) among all HST filters to get consistent colors. To do so, we identify isolated point sources in each cluster field, and we use the psfmatch task in IRAF to match all HST images to have the same PSF as the reddest band, F160W. In practice, because of the small field of view of each cluster and the crowded environment, we can only select ∼5 isolated point sources in each cluster for PSF-matching. However, we measure the curves of growth of each point source after PSF-matching and find that in most filters, the curves match to within ∼20% of that of F160W.

After pre-processing, we extract photometry on HST images using SExtractor (Bertin & Arnouts 1996; version 2.8.6). We use the combined IR images as the detection image, and the SExtractor detection/deblending settings similar to (but slightly more conservative than) the values adopted by CLASH for their high-z galaxy search (Postman et al. 2012). Because our focus in this work is on identifying IRAC-detected high-z sources, the slightly more conservative settings do not reject many potential IRAC-detected candidates but eliminates most spurious detections. We also follow the procedure outlined in Trenti et al. (2011) to rescale the flux errors reported by SExtractor. At the end of the process, we use the resulting photometric catalogs and segmentation maps for both IRAC photometry (Section 2.2) and color selection (Section 3).

The IRAC imaging data set for SURFS UP was presented in Paper I; the total exposure time for each IRAC channel is about 30 hr (see Ryan et al. 2014 and Paper I for details.) The coadded IRAC mosaics are deeper within the main cluster fields covered by WFC3/IR, and the typical 3σ limiting magnitude within 3''-radius apertures is 26.6 mag in the 3.6 μm channel (hereafter ch1) and 26.2 mag in the 4.5 μm channel (hereafter ch2) where source blending is not severe.

We use T-PHOT to measure consistent colors between HST and Spitzer IRAC images with a template-fitting approach (see also Laidler et al. 2007 for the template-fitting concept employed by T-PHOT.) The template-fitting approach has been demonstrated to work well for blank-field surveys such as CANDELS (Guo et al. 2013), but it does require images with zero mean background. Most galaxy cluster fields have considerable spatial variations in local sky background, so subtracting a constant background does not generally work. Therefore, instead of fitting all sources in the field at the same time—which is the strategy for blank-field surveys—we subtract the local background and perform the fit for each high-z candidate separately to get the cleanest residual possible.

As described in Merlin et al. (2015), T-PHOT is designed to measure the fluxes in the low-resolution image (in our case the IRAC images) for all the sources detected in the high-resolution image (in our case the F160W images). T-PHOT does so by constructing a template for each source; it convolves the cutout of each source in the F160W image by a PSF-transformation kernel that matches the F160W resolution to the IRAC resolution. Once the templates are available (and with their fluxes normalized to 1), T-PHOT solves the set of linear equations and finds the combination of coefficients for each template that most closely reproduce the pixel values in IRAC images; each coefficient is therefore the flux of the source in IRAC. T-PHOT also calculates the full covariance matrix and uses the diagonal terms of the covariance matrix to calculate flux errors. For each source, T-PHOT reports a "covariance index," defined as the ratio between the maximum covariance of the source with its neighbors (max(σ_ij)) over its flux variance (σ_ii), which serves as an indicator of how strongly correlated the source's flux is with its closest (or brightest) neighbor. Generally, a high covariance index (≳1) is associated with more severe blending and large flux errors, at least from simulations (Laidler et al. 2007; Merlin et al. 2015). Therefore, sources with high covariance indices should be treated with caution.

Obviously the PSF-transformation kernel that matches the F160W PSF to the IRAC PSF is a crucial element in this process. We generate IRAC PSFs by stacking point sources observed in the exposures from both the primary cluster field and the flanking field. We identify point sources using Sextractor with DEBLEND_MINCONT = 10⁻⁷, MINAREA = 9, DETECT_THRESH = ANALYSIS_THRESH = 2, and a Gaussian convolution kernel with σ = 3 pixels (defined over a 5 by 5-pixel grid). We require that all point sources have an axis ratio of b/a > 0.9, lie on the stellar locus within the box shown in Figure 1, and are sufficiently separated from neighboring objects to have reliable centroids (FLAGS ≤ 3). We recompute the PSF centroids by fitting a Gaussian profile to the inner profile (r < 4 pixels) using the Sextractor barycenters as initial guesses, and align the point sources using sinc interpolation. To mask neighboring objects, we grow the segmentation maps from Sextractor by 2 pixels. At each phase we subtract the local sky (assuming there are no local gradients) and normalize the flux of the point source to unity. We sigma-clip average the masked, registered, normalized point sources and do one further background correction only to ensure the convolutions with T-PHOT are flux conserving. As discussed in Paper I, our stacked PSFs are consistent with the IRAC handbook. Each of our clusters contains at least 40 point sources per bandpass in our PSF-making process.

Zoom In Zoom Out Reset image size

In practice, T-PHOT still breaks down in very crowded regions (e.g., near the cluster center or near bright cluster galaxies); in this case, we are limited mostly by our knowledge of the IRAC PSFs and our ability to subtract sky background underneath the sources. We also measure the "reduced χ²" for each source in IRAC within a 24 by 24 box by calculating the average difference per pixel between the model pixel values and the observed pixel values: ${\chi }_{\nu }^{2}={\displaystyle \sum }_{i}^{{N}_{\mathrm{pix}}}{({f}_{i,\mathrm{model}}-{f}_{i,\mathrm{obs}})}^{2}$/$({f}_{i,\mathrm{obs}}^{2}\times {N}_{\mathrm{pix}})$, where ${f}_{i,\mathrm{model}}$ and ${f}_{i,\mathrm{obs}}$ are the model (best-fit) and observed flux in pixel i, respectively. Later in this work, we only report the IRAC fluxes of the high-z galaxies with reliable IRAC flux measurements, i.e., ${\chi }_{\nu }^{2}\leqslant 3$.

For the sources with nominal signal-to-noise ratios (S/Ns) above 3 but with poor T-PHOT residuals, we do not trust the T-PHOT-measured fluxes and estimate the local 3σ flux limits via artificial source simulations. We insert artificial point sources into the F160W image (but not into the IRAC image) within 5'' of the high-z candidates and run T-PHOT to measure the local sky level. We repeat this process at least 100 times near each high-z candidate and use the resulting IRAC flux distribution to determine the 3σ flux limits. In our analysis in Section 5, we use the 3σ flux limit for only one IRAC filter for one source (ch1 for MACS 1423-1384); for all the other sources, their IRAC flux measurements in both IRAC channels pass the ${\chi }_{\nu }^{2}$ test.

We also run a separate set of simulations that independently estimate the magnitude errors in case T-PHOT underestimates the magnitude errors in crowded regions even when they pass the ${\chi }_{\nu }^{2}$ test. In this set of simulations, we insert fake point sources around each high-z target (within 5'') in IRAC images with the same magnitude as the T-PHOT-measured value, and measure the flux of the fake sources again with T-PHOT. We then calculate the median difference between the input and output magnitudes of the fake sources as an independent magnitude error estimate. We find that for most sources, the T-PHOT-reported magnitude error is within 0.1 mag from the simulated magnitude error, but sometimes the simulated magnitude error is much larger than the T-PHOT-reported value. In these cases, T-PHOT might underestimate the true magnitude errors, so we adopt the simulated magnitude errors in our SED modeling. We note that by adding fake sources in IRAC images, we increase the flux error due to crowding, so the simulated magnitude errors could be higher than the true magnitude errors.

For the six clusters that in the CLASH sample, we use the criteria below. To be selected as a z ∼ 6 LBG candidate, a source has to satisfy all of the following criteria from Gonzalez et al. (2011):

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\rm{F}}775{\rm{W}}-{\rm{F}}850\mathrm{LP}\gt 1.3\\ {\rm{F}}850\mathrm{LP}-{\rm{F}}125{\rm{W}}\lt 0.8\\ {\rm{S}}/{\rm{N}}\geqslant 5\quad \mathrm{in}\ {\rm{F}}850\mathrm{LP}\ \mathrm{and}\ {\rm{F}}125{\rm{W}}\\ {\rm{S}}/{\rm{N}}\lt 2\quad \mathrm{in}\ \mathrm{filters}\ \mathrm{bluer}\ \mathrm{than}\ {\rm{F}}606{\rm{W}},\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where we calculate all S/Ns within the isophotal (ISO) aperture.¹⁹ If the S/N in F775W is below one, we use the 1σ flux limit in F775W to calculate the ${\rm{F}}775{\rm{W}}-{\rm{F}}850\mathrm{LP}$ color. To select LBG candidates at z ∼ 7, we use the color criteria from Bouwens et al. (2011):

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\rm{F}}850\mathrm{LP}-{\rm{F}}105{\rm{W}}\gt 0.7\\ {\rm{F}}105{\rm{W}}-{\rm{F}}125{\rm{W}}\lt 0.45\\ {\rm{F}}850\mathrm{LP}-{\rm{F}}105{\rm{W}}\gt 1.4\times ({\rm{F}}105{\rm{W}}-{\rm{F}}125{\rm{W}})+0.42\\ {\rm{S}}/{\rm{N}}\geqslant 5\quad \mathrm{in}\ {\rm{F}}105{\rm{W}}\ \mathrm{and}\ {\rm{F}}125{\rm{W}}\\ {\rm{S}}/{\rm{N}}\lt 5\quad \mathrm{in}\ {\rm{F}}814{\rm{W}}\\ {\rm{S}}/{\rm{N}}\lt 2\quad \mathrm{in}\ \mathrm{filters}\ \mathrm{bluer}\ \mathrm{than}\ {\rm{F}}775{\rm{W}}.\end{array}\right.\end{eqnarray} \tag{ 2 }$$

To select the LBGs at z ∼ 8, we use the criteria from Bouwens et al. (2011):

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\rm{F}}105{\rm{W}}-{\rm{F}}125{\rm{W}}\gt 0.45\\ {\rm{F}}125{\rm{W}}-{\rm{F}}160{\rm{W}}\lt 0.5\\ {\rm{S}}/{\rm{N}}\geqslant 5\ \mathrm{in}\ {\rm{F}}125{\rm{W}}\ \mathrm{and}\ {\rm{F}}160{\rm{W}}\\ {\rm{S}}/{\rm{N}}\lt 2\ \mathrm{in}\ \mathrm{filters}\ \mathrm{bluer}\ \mathrm{than}\ {\rm{F}}850\mathrm{LP}.\end{array}\right.\end{eqnarray} \tag{ 3 }$$

Finally, to select the LBG candidates at z ≳ 9, we use the criteria from Zheng et al. (2014):

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\rm{F}}125{\rm{W}}-{\rm{F}}160{\rm{W}}\gt 0.8\\ {\rm{S}}/{\rm{N}}\geqslant 5\ \mathrm{in}\ {\rm{F}}160{\rm{W}}\\ {\rm{S}}/{\rm{N}}\lt 2\ \mathrm{in}\ \mathrm{filters}\ \mathrm{bluer}\ \mathrm{than}\ {\rm{F}}850\mathrm{LP}.\end{array}\right.\end{eqnarray} \tag{ 4 }$$

In total, we find 43 z ≳ 6 LBG candidates from the six CLASH/HFF clusters; among them, 13 have ≥3σ detections in at least one IRAC channel.

For the two galaxy clusters that are not in the CLASH sample (MACS 0454 and RCS 2327), we design our own selection criteria for z ≳ 6 galaxy candidates. For MACS 0454, we have HST imaging data in F555W, F775W, F814W, F850LP, F110W, and F160W, although the images in F775W and F850LP are shallower than the other filters and we use both filters only for S/N rejection of low-z interlopers. We use the following criteria to select 6 ≲ z ≲ 7.5 galaxy candidates (F814W-dropouts):

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\rm{F}}814{\rm{W}}-{\rm{F}}110{\rm{W}}\geqslant 1.0\\ {\rm{F}}110{\rm{W}}-{\rm{F}}160{\rm{W}}\leqslant 0.3\\ {\rm{S}}/{\rm{N}}\geqslant 5\ \mathrm{in}\ {\rm{F}}110{\rm{W}}\ \mathrm{and}\ {\rm{F}}160{\rm{W}}\\ {\rm{S}}/{\rm{N}}\lt 2\ \mathrm{in}\ {\rm{F}}555{\rm{W}}.\end{array}\right.\end{eqnarray} \tag{ 5 }$$

We show the HST color–color diagram of the F814W-dropout selection in Figure 2. A total of ten sources satisfy the color and S/N cuts listed above; among them, MACS 0454-1251 and MAC 0454-1817 are detected in at least one IRAC channel.

Zoom In Zoom Out Reset image size

Because we have only one filter redward of F110W, we refrain from searching for F110W-dropouts in MACS 0454 as it would yield objects detected in only one HST filter.

For RCS 2327, we have deep HST images in F435W, F814W, F098M, F125W, and F160W, so we use the following criteria for 6 ≲ z ≲ 7.5 galaxy candidates (F814W-dropouts):

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\rm{F}}814{\rm{W}}-{\rm{F}}098{\rm{M}}\geqslant 2.2\ \\ {\rm{F}}098{\rm{M}}-{\rm{F}}160{\rm{W}}\leqslant 1.6\ \\ {\rm{F}}814{\rm{W}}-{\rm{F}}098{\rm{M}}\geqslant 2.0\times ({\rm{F}}098{\rm{M}}-{\rm{F}}160{\rm{W}})+1.0\\ {\rm{S}}/{\rm{N}}\geqslant 5\ \mathrm{in}\ \ {\rm{F}}098{\rm{M}}\ \ \mathrm{and}\ \ {\rm{F}}160{\rm{W}}\\ {\rm{S}}/{\rm{N}}\lt 3\ \mathrm{in}\ \ {\rm{F}}435{\rm{W}}.\end{array}\right.\end{eqnarray} \tag{ 6 }$$

When the S/N in F814W is less than unity, we use its 1σ flux limit to calculate colors. We demonstrate the F814W-dropout selection from RCS 2327 in Figure 3, and six sources pass the above color and signal-to-noise cuts. One of the six sources, RCS 2327-1218, is detected in both IRAC channels. In addition to the above criteria, we also use the criteria similar to the BoRG survey (Trenti et al. 2011) to search for z ≳ 7.5 galaxy candidates:

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\rm{F}}098{\rm{M}}-{\rm{F}}125{\rm{W}}\geqslant 1.75\\ {\rm{S}}/{\rm{N}}\geqslant 5\ \mathrm{in}\ \ {\rm{F}}125{\rm{W}}\ \ \mathrm{and}\ \ {\rm{F}}160{\rm{W}}\\ {\rm{S}}/{\rm{N}}\lt 3\ \mathrm{in}\ \ {\rm{F}}814{\rm{W}}\ \ \mathrm{and}\ \ {\rm{F}}435{\rm{W}}.\end{array}\right.\end{eqnarray} \tag{ 7 }$$

and when the S/N in F098M is less than unity, we use its 1σ flux limit to calculate colors. The F098M-dropout search yields no galaxy candidate, so we find a total of one galaxy candidate with IRAC detections at z ≳ 6 from RCS 2327 (see also Hoag et al. 2015 for more details on the dropout search in RCS 2327).

Zoom In Zoom Out Reset image size

To summarize, we find a total of 69 LBG candidates at z ≳ 6 from 9 clusters in SURFS UP; 17 of them have IRAC detections in at least one channel. Figures 4 and 5 show the cutouts of the 16 IRAC-detected LBG candidate in HST and Spitzer images (one candiate was reported by Ryan et al. 2014). We also report the IRAC photometry for the entire sample in Table 3. We use the simulated IRAC magnitude errors for MACS 1423-1384, MACS 1423-587, MACS 0744-2088, MACS 1423-2097, and MACS 0454-1817 because we find that T-PHOT likely underestimates their IRAC magnitude errors from the simulations. We keep MACS 1423-1384 in our sample because it has a nominal 5.9σ detection in ch2 from T-PHOT, although the simulated magnitude error suggests that the additional flux error due to crowding reduces it to a 2.2σ detection in ch2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We selected this object as a z ∼ 7 LBG candidate for spectroscopic follow-up on 2013 April 3 and 2014 May 26. This source was also selected by both Smit et al. (2014) and Bradley et al. (2014) as a z_phot ∼ 6.7 LBG with high [O iii]+Hβ equivalent widths (>1300 Å rest frame). We used the 830 l mm⁻¹ grating in the 2013 run and the 1200 l mm⁻¹ grating in the 2014 run, and the total integration times are roughly 6000 and 7200 s, respectively. We had good (but not photometric) conditions with ∼1'' seeing in the 2013 run, but the conditions were highly variable in the 2014 run. Therefore, we only present the line flux measurements from the 2013 run, though the line was detected significantly in both runs.

In Figure 6 we show the two-dimensional spectra of RXJ 1347-1216 from both observation runs (in the top and middle panels) and the combined one-dimensional spectrum (in the bottom panel). We detect an extended emission feature with FWHM_obs ∼ 16.5 Å in the 2013 run, and although the blue side of the feature is severely contaminated by sky line residual, its asymetric profile with a tail to the red side of the spectrum strongly suggests that it is Lyα. Using the centroid of the sky line residual at 9439 Å as the peak of of line profile, we determine its Lyα redshift as z_Lyα = 6.76 ± 0.003 (the uncertainty corresponds to the width of the sky line residual). The Lyα redshift is in excellent agreement with its photometric redshift z_phot = 6.8 ± 0.1, lending additional support to the identification of the Lyα feature.

Zoom In Zoom Out Reset image size

The emission feature is also independently detected at 4σ in the GLASS grism data at ∼9440 Å. With the grism spectra in both G102 and G141, one can test the possibility of this feature being an [O ii] doublet (λ3727, λ3729) at z ∼ 1.5 by looking for the [O iii] λ5007 line at ∼1.3 μm. For a typical star-forming galaxy, the line ratio [O iii]/[O ii] should be at least ∼0.3 (Jones et al. 2015), and the [O iii]/[O ii] ratio is even higher for low-metallicity galaxies (e.g., Maiolino et al. 2008). Assuming the detected line in DEIMOS is [O ii] at z = 1.53, HST G141 grism data imply a 2σ upper limit on [O iii]/[O ii] ≲ 0.3, which is highly unlikely for a star-forming galaxy (Schmidt et al. 2015). Therefore, we conclude that the HST grism data also strongly support the z = 6.76 Lyα interpretation of this emission line.

We perform line flux measurements from the DEIMOS data obtained during the 2013 run (with the 830 l mm⁻¹ grating) following the procedure outlined in Section 2.4 of Lemaux et al. (2009). In short, we place two filters of width 20 Å on both sides of the emission line that are free of spectral features and sky line residuals to measure the background, and a central filter encompassing the emission line to measure the integrated flux. The background, which is fit to a linear function, is then subtracted from the integrated line flux. We perform spectrophotometric calibration of the DEIMOS data using two other compact sources (with half-light radii ∼03 as measured from the F775W image) on the same mask with continuum detection in the following manner: for each object, the combination of slit loss and loss due to clouds was determined by calculating a spectral magnitude, done by correcting each DEIMOS spectrum for spectral response and atmospheric extinction and convolving the F775W filter curve with the resulting spectra, and comparing this value with the magnitude measured in the HST image. The ratio of the flux densities for each of the two sources is calculated and averaged to estimate the total spectral loss for this mask, which is then applied to RXJ 1347-1216 assuming a similar half-light radius for this object. The reason for this assumption is, though the half-light radius of RXJ 1347-1216 is smaller (02), which would, in principle, mean less slit loss, the size of the Lyα nebula is known to far exceed the size of the UV continuum region (e.g., Wisotzki et al. 2015). For the two sources on the mask with which we performed the flux calibration, the total measured throughput of the slit was ∼40%, lower than the ∼60% expected for sources of this size (if they are symmetric), suggesting at least some departure from a photometric night.

Using the above procedure, we measure the line flux from the 830 l mm⁻¹ data to be 7.8 ± 0.7 × 10⁻¹⁸ erg cm⁻² s⁻¹, which translates into a Lyα luminosity of 4.1 ± 0.4 × 10⁴² erg s⁻¹. We do not detect continuum in the spectrum, so we estimate the rest-frame Lyα equivalent width using the object's broadband magnitude in F105W (on the red side of Lyα) to be 26 ± 4 Å. The equivalent width uncertainty include the Poisson noise in the central filter encompassing the sky line residual, the uncertainty in the continuum, and the uncertainty in DEIMOS absolute flux calibration.

We note that our measured Lyα line flux from DEIMOS data is roughly a factor of 3 lower than that measured from the HST grism data, which is 2.6 ± 0.5 × 10⁻¹⁷ erg cm⁻² s⁻¹ (Schmidt et al. 2015). The difference in the line flux measurements could be due to the following factors: (1) our DEIMOS data are shallower than the HST grism data (∼2 hr of total integration time in the 2013 mask, versus ∼5 orbits of HST integration time in G102 grism), with the difference in depths leading to differences in the spatial and extent of the emission detected above the noise in each observation, which impacts the total integrated line flux; (2) we might still underestimate the Lyα slit loss because the true extent of the Lyα-emitting region is much larger than the continuum-emitting region, while the HST slitless grism recovers more of the Lyα flux; (3) there is a sky line coincident with the peak wavelength of the Lyα emission in the DEIMOS data which appears slightly over-subtracted that serves to slightly lower the measured line flux; and (4) there could also be issues with contamination subtraction in the HST grism data, although it is unclear which direction it biases line flux measurement. We do expect the ground-based line flux measurement to be a lower limit to the space-based measurement, which is consistent with the measured values from DEIMOS and from HST grism data.

We also measure the inverse line asymmetry 1/a_λ as defined in Lemaux et al. (2009), and the line's inverse asymmetry value (1/a_λ = 0.22) is well within the range (1/a_λ < 0.75) typical for convincing Lyα emission.

Using Lyα line flux, one can also infer a SFR following Lemaux et al. (2009). The inferred SFR is strictly a lower limit, because our conversion assumes no attenuation of Lyα photons by dust or neutral hydrogen. The inferred SFR is 1.6 ± 0.1 M_⊙ yr⁻¹, consistent with being a lower limit to the value we derive from SED fitting in Section 5, 17.0 ± 0.5 M_⊙ yr⁻¹. The Lyα-inferred SFR roughly yields a Lyα escape fraction of ∼10%.

We observed MACS 0454 with DEIMOS on 2014 November 28 and 29, using the 1200 l mm⁻¹ grating on both nights. The total exposure time for this mask is 7200 s. We also reduce the DEIMOS data and extract one-dimensional spectrum following the procedure in Lemaux et al. (2009), in the same way as RXJ 1347-1216.

We show the reduced two-dimensional spectra of MACS 0454-1251 in Figure 7, and an extended emission feature is clearly detected on both nights. However, a bright sky line residual cuts through the middle of the emission feature in the spectral direction and makes the line interpretation ambiguous. Given the width of the emission line, it could be either Lyα at z = 6.32 or [O ii] at z = 1.39, but the sky line residual makes it difficult to either confirm or rule out Lyα or [O ii] based on line shape alone. However, as we will show in Section 5.2, this line is more likely to be Lyα than [O ii] because its HST fluxes (and upper limits) are much better fit by a galaxy template at z = 6.32 than at z = 1.39 (Section 5.2), and its photometric redshift probability density function P(z) has very a low probability at z = 1.39 (see the P(z) curve in Figure 10). We measure the line fluxes from both nights to be 1.2 ± 0.2 × 10⁻¹⁷ erg cm⁻² s⁻¹ (first night) and 8.0 ± 1.5 × 10⁻¹⁸ erg cm⁻² s⁻¹ (second night), and these translate to Lyα luminosities of 5.5 ± 0.9 × 10⁴² erg s⁻¹ (first night) and 3.6 ± 0.3 × 10⁴² erg s⁻¹ (second night). From the line fluxes, we estimate its rest-frame equivalent widths (assuming Lyα) using the continuum on the red side of Lyα (estimated from its F110W flux density) from both nights to be 8.2 ± 1.4 and 5.4 ± 1.0 Å. We also infer SFR lower limits to be 2.3 ± 0.4 M_⊙ yr⁻¹ (first night) and 1.5 ± 0.3 M_⊙ yr⁻¹ (second night) based on Lyα fluxes, also fully consistent with being lower limits to the value derived from SED fitting in Section 5, ${17.0}_{-4.1}^{+18.0}\;{M}_{\odot }\;{\mathrm{yr}}^{-1}$. The Lyα-inferred SFR also roughly corresponds to an escape fraction of ∼10% for Lyα photons.

Zoom In Zoom Out Reset image size

The inverse asymmetry value measured for the one-dimensional emission feature is ∼0.5 for both masks, consistent with the values typical for Lyα. However, the line asymmetry estimate is less reliable than that of RXJ 1347-1216 because we have to mask the over-subtracted skyline near the central wavelength of the emission feature. Based on the object's photometric information and line asymmetry measurements, we identify this source as a Lyα emitter at z = 6.32, although we are less confident with this Lyα interpretation than we are with RXJ 1347-1216. MACS 0454 is not in the GLASS sample, so we are unable to cross check our measurements with HST grism data.

For photometric redshift and stellar population modeling, we use the photometric redshift code EAZY (Brammer et al. 2008) with stellar population templates from Bruzual & Charlot (2003, BC03) with Chabrier (2003) initial mass function (IMF) between 0.1 and 100 M_⊙ and a metallicity of 0.2 Z_⊙. There is very little direct observational evidence for galaxy metallicity at z ≥ 3, but limited results so far suggest that the majority of them have sub-solar metallicity (Maiolino et al. 2008); we choose 0.2 Z_⊙ for easy comparison with other works. The galaxy templates are generated assuming an exponentially declining star formation history with e-folding time τ ranging from 0.1 and 30 Gyr, and ages of the stellar population range from 10 Myr to 13 Gyr. For each combination of age and τ, we implement galaxy internal dust attenuation using the Calzetti et al. (2000) prescription, with the reddening parameter $E(B-V)$ ranging from 0 to 1 mag to include potential low-z dusty galaxy solutions. The $E(B-V)$ grid we use have a step size of ${\rm{\Delta }}E(B-V)=0.02$ mag from $E(B-V)=0$ to 0.5 mag and a step size of ${\rm{\Delta }}E(B-V)=0.1$ mag from $E(B-V)=0.5$–1 mag. We also use the stellar templates from the photometric code Le Phare (Ilbert et al. 2006)²⁰ in the fitting to check if stellar templates provide significantly better fits.

Recent studies have shown that for some galaxy candidates, strong nebular emission lines contribute significantly to broadband fluxes and therefore influence the inferred galaxy properties (e.g., Schaerer & De Barros 2010; Smit et al. 2014). Therefore, we use galaxy templates that include nebular emission lines in the modeling. In order to calculate the expected line fluxes for a given BC03 galaxy template, we calculate the integrated Lyman continuum flux (before dust attenuation) and use the relation from Leitherer & Heckman (1995) to calculate the expected fluxes from hydrogen recombination lines (mainly Hα, Hβ, Paβ, and Brγ) while assuming the Lyman continuum escape fraction to be zero. Non-zero Lyman continuum escape fraction will reduce the strength of optical nebular emission lines (see Inoue 2011 and Salmon et al. 2015). We then use the tabulated line ratios between Hβ and the metal lines from Anders & Alvensleben (2003) to calculate the metal line fluxes for a metallicity of 0.2Z_⊙. For templates with dust attenuation, we include the dust attenuation effects after adding nebular emission lines. The resultant equivalent widths as a function of galaxy age, for τ = 100 Myr and $E(B-V)=0$, are shown in Figure 8, in agreement with Leitherer & Heckman (1995).

Zoom In Zoom Out Reset image size

In addition to nebular emission lines, we also include nebular continuum emission that account for the bound-free emission of H i and He i as well as the two-photon emission of hydrogen from the 2 s level. We follow the prescription in Krueger et al. (1995, their Equations (7) and (8)) to calculate the nebular continuum flux as a function of Lyman continuum photon density, and we calculate the emission coefficients using the methods in Brown & Mathews (1970) and Nussbaumer & Schmutz (1984).²¹ Nebular continuum emission could be an important component for very young (∼10 Myr) starbursts and can contribute up to ∼1/3 of the total continuum just blueward of the rest-frame 4000 Å break. Nebular continuum emission also makes the rest-frame UV slope redder than expected from stars alone (Schaerer & De Barros 2010).

We do not include Lyα in our galaxy templates. Strong Lyα emission could affect the LBG color selection by changing the rest-frame UV broadband colors and could affect the derived physical properties from SED fitting (Schaerer et al. 2011; De Barros et al. 2014). But Lyα photons suffer from complicated radiative transfer processes and does not show tight correlations with the stellar population properties, so for simplicity we do not include Lyα emission in our modeling.

Strong gravitational lensing boosts galaxy fluxes and increases the apparent SFR and stellar mass. To calculate the unlensed SFR and stellar mass, we use the magnification factor μ_best estimated from the cluster mass models for each galaxy candidate at its redshift. We generate our own models for MACS 1149, MACS 0717, MACS 0454, RXJ 1347, and RCS 2327, following the procedures outlined in Bradač et al. (2005, 2009). In short, we constrain the gravitational potential on a mesh grid within a galaxy cluster field via χ² minimization, and we adaptively use denser pixel grids near the core(s) of the cluster and around multiple images. We find the minimum χ² values by iteratively solving a set of linearized equations that satisfy ∂χ²/∂ψ_k = 0, where ψ_k is the gravitational potential in the kth dimension. We then produce the magnification (μ) map from the best-fit gravitational potential map. For the rest of the clusters in our sample (MACS 1423, MACS 2129, and MACS 0744), we use the public PIEMD-eNFW²² models by Zitrin et al. (2013).²³ The magnification factors (and their errors) are estimated at the galaxy candidate positions from the z = 9 magnification maps except for MACS 1423-587, MACS 1423-774, MACS 1423-2248, and RXJ 1347-1800 (which have z_best ∼ 1 so we estimate their μ_best from the z = 1 magnification map).

We list the best-fit galaxy properties in Table 4 and show the best-fit templates and the photometric redshift probability density function P(z) (while allowing redshift to float) in Figures 9 and 10. For each galaxy candidate, we estimate the statistical uncertainties of stellar population properties using Monte Carlo simulations: we perturb the photometry within the errors (assuming Gaussian flux errors), re-fit with the same set of galaxy templates, and collected the distributions of each best-fit property. We only perturb the fluxes where S/N ≥ 1. For upper limits, we do not perturb the fluxes in our simulations. The systematic errors related to assumptions in IMF, galaxy metallicity, and the functional form of star formation history are not represented by the error bars. We show the distributions from Monte Carlo simulations for stellar mass, SFR, and stellar population age in Appendix. From these distributions, we derive the confidence intervals that bracket 68% of the total probability in Table 4. The error bars do not include uncertainties in μ_best.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  Photometric Redshift and Stellar Population Modeling Results

Object ID	zbesta	${\mu }_{\mathrm{best}}$ b	Mstellar × fμc	SFR × fμc	Aged	sSFRe	$E{(B-V)}_{\mathrm{fit}}$ f	βg	${M}_{1600}-2.5\mathrm{log}({f}_{\mu })$ h
			(109 M⊙)	(M⊙ yr−1)	(Myr)	(Gyr−1)	(mag)		(mag)
F125W-dropouts (z ∼ 9; ${M}_{\mathrm{UV}}^{*}$ = −20.63 ± 0.36 at z ∼ 8 from Bouwens et al. 2015)
MACS 1149-JD	${9.3}_{-0.1}^{+0.1}$	${5.5}_{-0.3}^{+0.3}$	${0.9}_{-0.4}^{+0.5}$	${1.9}_{-0.5}^{+0.6}$	${200}_{-60}^{+250}$	${2.1}_{-0.7}^{+3.1}$	0.00	⋯	−19.9 ± 0.1
F105W-dropouts (z ∼ 8; ${M}_{\mathrm{UV}}^{*}$ = −20.63 ± 0.36 at z ∼ 8 from Bouwens et al. 2015)
RXJ 1347-1080	${7.3}_{-0.4}^{+0.3}$	${6.0}_{-0.6}^{+0.6}$	${0.6}_{-0.4}^{+0.2}$	${0.5}_{-0.5}^{+1.0}$	${290}_{-260}^{+350}$	${0.9}_{-0.9}^{+1.9}$	0.00	$-{2.5}_{-1.1}^{+1.2}$	−18.9 ± 0.2
MACS 1423-1384	${6.9}_{-0.1}^{+0.9}$	${4.0}_{-1.2}^{+1.2}$	${1.1}_{-0.6}^{+3.2}$	${115.6}_{-114.6}^{+16.1}$	≤130	${105.0}_{-98.9}^{+0.1}$	0.38	$-{0.5}_{-0.6}^{+1.6}$	−19.4 ± 0.2
F850LP-dropouts (z ∼ 7; ${M}_{\mathrm{UV}}^{*}$ = −20.87 ± 0.26 from Bouwens et al. 2015)
MACS 1423-1494	${7.1}_{-0.5}^{+0.3}$	${1.7}_{-0.1}^{+0.1}$	${0.2}_{-0.1}^{+1.2}$	${22.1}_{-19.3}^{+4.3}$	≤10	${105.1}_{-9.2}^{+0.1}$	0.14	$-{2.1}_{-1.4}^{+1.1}$	−20.1 ± 0.2
MACS 0744-2088	${7.0}_{-0.1}^{+0.2}$	${1.5}_{-0.1}^{+0.1}$	${0.7}_{-0.4}^{+1.2}$	${34.0}_{-28.1}^{+4.0}$	${20}_{-0}^{+180}$	${45.7}_{-42.3}^{+16.8}$	0.16	$-{0.8}_{-1.2}^{+0.7}$	−20.7 ± 0.1
RXJ 1347-1216i	6.76	${5.0}_{-0.5}^{+0.5}$	${0.2}_{-0.1}^{+0.1}$	${17.0}_{-2.7}^{+2.6}$	≤10	${105.0}_{-0.1}^{+0.1}$	0.20	$-{2.5}_{-1.0}^{+0.7}$	−19.3 ± 0.1
MACS 1423-2097	${6.8}_{-0.2}^{+0.1}$	${3.5}_{-0.1}^{+0.1}$	${2.9}_{-2.5}^{+0.3}$	${1.5}_{-0.1}^{+20.0}$	${720}_{-490}^{+90}$	${0.5}_{-0.1}^{+90.9}$	0.00	$-{0.6}_{-1.1}^{+0.6}$	−19.7 ± 0.1
Bullet-3j	${6.8}_{-0.1}^{+0.1}$	${12}_{-4.0}^{+4.0}$	${2.0}_{-0.8}^{+0.6}$	${1.3}_{-0.6}^{+1.4}$	${630}_{-230}^{+160}$	${0.7}_{-0.6}^{+0.5}$	0.00	⋯	−18.9 ± 0.4
MACS 1423-587	${0.1}_{-0.1}^{+6.7}$	${1.5}_{-0.1}^{+0.1}$	${0.1}_{-0.1}^{+0.2}$	≤0.1	${11500}_{-11490}^{+1250}$	≤0.1	0.70	⋯	⋯
RXJ 1347-1800	${0.8}_{-0.5}^{+0.7}$	${4.1}_{-0.1}^{+0.1}$	${0.1}_{-0.1}^{+0.1}$	≤1.2	${2600}_{-2570}^{+7830}$	≤36.0	0.00	⋯	⋯
MACS 1423-774	${1.2}_{-0.3}^{+5.5}$	${1.2}_{-0.1}^{+0.1}$	${0.3}_{-0.2}^{+0.5}$	≤20.6	${1430}_{-1420}^{+180}$	≤100.1	0.06	⋯	⋯
MACS 1423-2248	${1.2}_{-1.0}^{+0.3}$	${1.2}_{-0.1}^{+0.1}$	${1.0}_{-1.0}^{+2.7}$	≤5.2	${5000}_{-4980}^{+3640}$	≤52.9	0.00	⋯	⋯
F814W-dropouts (z ∼ 6–7; ${M}_{\mathrm{UV}}^{*}$ = −20.87 ± 0.26 from Bouwens et al. 2015)
RCS 2327-1282	${7.7}_{-0.4}^{+0.1}$	${4.1}_{-0.4}^{+0.5}$	${1.9}_{-0.5}^{+0.4}$	${6.1}_{-0.2}^{+0.5}$	${400}_{-120}^{+100}$	${3.2}_{-0.4}^{+1.6}$	0.00	$-{3.0}_{-0.3}^{+0.2}$	−20.9 ± 0.1
MACS 0454-1817	${6.5}_{-0.1}^{+0.2}$	${2.6}_{-0.3}^{+0.3}$	${0.4}_{-0.1}^{+3.8}$	${38.6}_{-36.3}^{+28.9}$	≤30	${105.0}_{-67.6}^{+0.1}$	0.28	$-{1.4}_{-0.8}^{+0.5}$	−19.3 ± 0.1
MACS 0454-1251	${6.1}_{-0.1}^{+0.1}$	${4.4}_{-0.4}^{+0.4}$	${2.1}_{-0.8}^{+2.8}$	${17.9}_{-5.8}^{+14.1}$	${90}_{-10}^{+40}$	${8.7}_{-2.5}^{+4.8}$	0.12	$-{1.6}_{-0.2}^{+0.2}$	−20.8 ± 0.1
MACS 0454-1251k	6.32	${4.4}_{-0.4}^{+0.4}$	${0.5}_{-0.2}^{+2.9}$	${19.0}_{-5.4}^{+7.9}$	${30}_{-0}^{+0}$	${39.6}_{-35.3}^{+65.5}$	0.08	$-{1.6}_{-0.2}^{+0.2}$	−20.9 ± 0.1
F775W-dropouts (z ∼ 6; ${M}_{\mathrm{UV}}^{*}$ = −20.94 ± 0.20 from Bouwens et al. 2015)
MACS 1149-274	${5.8}_{-0.1}^{+0.1}$	${1.6}_{-0.1}^{+0.1}$	${2.4}_{-0.5}^{+0.5}$	${17.7}_{-0.1}^{+5.2}$	${100}_{-20}^{+10}$	${7.3}_{-1.2}^{+3.2}$	0.08	$-{1.6}_{-0.1}^{+0.1}$	−21.2 ± 0.1
MACS 1149-1204	${5.7}_{-0.1}^{+0.1}$	${1.8}_{-0.1}^{+0.1}$	${1.3}_{-0.8}^{+0.3}$	${12.9}_{-2.0}^{+4.7}$	${80}_{-50}^{+30}$	${10.2}_{-2.8}^{+29.2}$	0.08	$-{1.5}_{-0.5}^{+0.4}$	−20.7 ± 0.2

Notes.

^a ${z}_{\mathrm{best}}$ is the photometric redshift using BC03 galaxy templates except for RXJ 1347-1216 and MACS 0454-1251, for which we identify Lyα emission at z_spec = 6.76 and z_spec = 6.32, respectively. ^bLensing magnification factor estimated from the galaxy cluster mass models mentioned in Section 5.1. For MACS 1423-587, MACS 1423-774, MACS 1423-2248, and RXJ 1347-1800, we estimate their μ_best from the magnification map at z = 1. ^cThe intrinsic stellar mass and SFR assuming $\mu ={\mu }_{\mathrm{best}}$. To use a different magnification factor μ, simply use ${f}_{\mu }\equiv \mu /{\mu }_{\mathrm{best}}$, where ${\mu }_{\mathrm{best}}$ is the best magnification factor we adopt for each object. When the best-fit SFR is zero, we report the 68% upper limit. ^dTime since the onset of star formation. For the sources with best-fit age equal to 10 Myr, the youngest template allowed in our models, we report the 68% upper limit of the age from Monte Carlo simulations. ^eSpecific star formation rate ≡ SFR / stellar mass. When the best-fit sSFR is zero, we report the 68% upper limit. ^fThe best-fit color excess $E(B-V)$ of the stellar emission from our SED modeling. The dust attenuation at rest-frame 1600 Å can be calculated using the dust attenuation curve from Calzetti et al. (2000) as ${A}_{1600}=9.97\times E(B-V)$. ^gMeasured rest-frame UV slope β from HST/WFC3 broadband fluxes, assuming the candidates are at z = z_best. We do not measureβ for MACS 1423-587, MACS 1423-774, MACS 1423-2248, and RXJ 1347-1800 because they have z_best ∼ 1. We also do not measure β for MACS 1149-JD because it has only 1 filter (F160W) that samples the rest-frame UV continuum. ^hRest-frame 1600 Å absolute magnitude assuming $z={z}_{\mathrm{best}}$ and $\mu ={\mu }_{\mathrm{best}}$. ⁱAll fits are performed at z = 6.76, the Lyα redshift. The confidence intervals reflect the maximal range returned from the simulations because the distributions for this object are highly skewed. ^jFirst reported by Ryan et al. (2014); included here for completeness. ^kAll fits are performed at z = 6.32, its Lyα redshift.

Download table as:  ASCIITypeset image

We also show the best-fit galaxy properties for RXJ 1347-1216 and MACS 0454-1251, the two galaxies for which we have line detections from DEIMOS data (see Section 4), when we fix their redshifts at their Lyα redshifts in Figure 11 (assuming both lines are Lyα). For RXJ 1347-1216, its photometric redshift is already sharply peaked at z = 6.7, so the best-fit template and physical properties do not change after fixing its redshift at the Lyα redshift. On the other hand, MACS 0454-1251 has slightly different best-fit photometric redshift and Lyα redshift, and we also list its best-fit properties at z_Lyα = 6.32 in Table 4. We also show the best-fit template at z = 1.39 in Figure 11 if the detected emission line is [O ii] instead of Lyα and see that the z = 1.39 solution has a higher ${\chi }_{\nu }^{2}$ (${\chi }_{\nu }^{2}=4.02$) than the z = 6.32 solution (${\chi }_{\nu }^{2}=2.90$). The likelihood ratio of these two fits, calculated using the total χ² as ${e}^{-{\rm{\Delta }}{\chi }^{2}}$, suggests that the z = 6.32 model is ∼8000 times more likely than the z = 1.36 model.

Zoom In Zoom Out Reset image size

The best-fit stellar masses in our sample range from 2 × 10⁸ M_⊙ to 2.9 × 10⁹ M_⊙ after correcting for magnification by the foreground clusters and excluding the z ∼ 1 interlopers. The stellar masses inferred from SED fitting have smaller statistical errors when HST fluxes are combined with IRAC fluxes because of the constraints on rest-frame optical emission from IRAC. We show the range of stellar mass from our Monte Carlo simulations in the Appendix (Figure 15), and we find that including IRAC fluxes tighten up the possible range of stellar mass for every object. We also see that the range of stellar mass spanned by our IRAC-detected sample are not necessarily at the high-mass end of the observed (Gonzalez et al. 2011) and simulated (Katsianis et al. 2015) stellar mass functions at z ≳ 6. In fact, the IRAC [3.6]–[4.5] color for several of our galaxy candidates (e.g., RXJ 1347-1216) suggest extremely young stellar population ages (∼10 Myr) and large equivalent widths from nebular emission lines. For these sources, stellar continuum emission might not dominate the observed IRAC fluxes, hence their true stellar masses depend sensitively on the equivalent widths of nebular emission lines. This demonstrates the combined power of strong gravitational lensing and deep IRAC images that allows one to measure the stellar mass of z ≳ 6 galaxies further down the stellar mass function.

On the other hand, the SFRs and stellar population ages are not necessarily well constrained by SED fitting even when IRAC fluxes are included (see Figures 16 and 17 in the Appendix). The SFRs of high-z galaxies are often calculated from their rest-frame UV fluxes (after correcting for dust attenuation), and these are often the only constraints available from observations. However, the UV-derived SFR depends critically on the amount of dust attenuation inside each galaxy, and the effect of dust on the rest-frame UV color is degenerate with the effect of stellar population age. Furthermore, UV-derived SFRs probe the star formation activity over the past ∼100 Myr, so it could underestimate the instantaneous SFR if the stellar population is younger than ∼100 Myr; for these systems, nebular emission line fluxes (e.g., Hα or [O ii]) are better proxies for SFRs (Kennicutt 1998). We consider SFRs and stellar population ages more poorly constrained compared with stellar mass, and we will discuss the degeneracies in SED fitting in Section 5.3.

IRAC fluxes also reveal four z ∼ 1 interlopers from our sample—MACS 1423-587, MACS 1423-774, MACS 1423-2248, and RXJ 1347-1800—as shown in the three bottom-right panels in Figure 10. All four sources have significant integrated probabilities at z ≥ 6 when only HST photometry is used in the fitting, but the addition of their IRAC fluxes pushes their photometric redshifts down to z ∼ 1, suggesting that the observed breaks between F850LP and F105W are more likely the rest-frame 4000 Å break instead of the Lyα break. This demonstrates the value of IRAC detections in discriminating between genuine z ≳ 6 galaxies and lower-z interlopers.

In Figures 9 and 10 we also show the best-fit stellar templates from Le Phare. Most of the sources are better fit by galaxy templates than by stellar templates, although there is only one case, RXJ 1347-1800, where both templates provide similarly good fits (${\chi }_{\nu }^{2}=0.57$ for galaxy templates and ${\chi }_{\nu }^{2}=0.62$ for stellar templates). For all of the galaxy candidates, their best-fit stellar templates are either brown dwarfs or low-mass stars from Chabrier & Baraffe (2000). We also check if our sample contains X-ray detected sources that could have significant contributions from AGNs and find no evidence that any of our galaxy candidates have AGNs. However, low-level AGN activity is still possible, and the stellar population properties we infer from SED modeling will depend on how much (if any) AGN contribution there is to their broadband fluxes.

The biases and uncertainties of SED fitting are well documented in the literature (e.g., Papovich et al. 2001; Lee et al. 2009), and most sources of systematic errors come from model assumptions. In general, stellar mass has the smallest systematic errors (∼0.3 dex), but uncertainties in galaxy star formation history can lead to large biases in galaxy age and SFR (Lee et al. 2009). In additional to the usual culprits of systematic errors (e.g., star formation history, IMF), another important source of systematic error is the nebular emission line ratios. We find that nebular emission line contributions to broadband fluxes are important for a subset of our sample, but direct observational constraints of rest-frame optical nebular emission line ratios of z > 6 galaxies will not arrive until the launch of the James Webb Space Telescope. Uncertainties in the amount of dust attenuation also complicates the interpretation of the best-fit parameters, as we demonstrate below.

We use MACS 1149-JD to show how uncertainties in dust attenuation can lead to uncertainties in SFR in Figure 12. We estimate the number density contours from our Monte Carlo simulations (with 1000 realizations) when IRAC fluxes are included in the modeling. We show the number density contour of each $E(B-V)$ value in the central panel and the marginalized histograms for galaxy age (on top) and SFR (on the right). The SFR histogram shows a single peak at 1.9 M_⊙ yr⁻¹ assuming a magnification factor of 5.5, but there is a long tail to higher SFRs that extends one order of magnitude. The long tail in SFR corresponds to higher dust attenuation templates ($E(B-V)\gt 0.1;$ green dashed contours in the central panel) compared to the peak of the histogram ($E(B-V)\lt 0.1;$ solid blue contours in the central panel). We note that because of the high photometric redshift of MACS 1149-JD (z_phot = 9.3), its UV continuum between rest-frame 1250 and 2600 Å is not well sampled by HST and IRAC photometry; the addition of K-band photometry should help constrain the amount of dust attenuation inside this galaxy.

Zoom In Zoom Out Reset image size

As a model-independent check on the inferred dust attenuation, we measure the rest-frame UV slope β for each galaxy candidate and list the values in Table 4. We can only measure β when a galaxy candidate has at least two filters sampling the UV continuum (between rest-frame 1250 and 2600 Å); therefore, we do not measure β for MACS 1149-JD (which only has F160W that samples UV continuum) or for MACS 1423-587, RXJ 1347-1800, MACS 1423-774, and MACS 1423-2248 (which have photometric redshifts ∼1). We measure β by using a power-law spectrum f_λ ∝ λ^β, convolving the spectrum with the filter curves that sample the UV continuum, and finding the β that best matches the observed fluxes. The uncertainties are quantified in bootstrap Monte Carlo simulations, and we show the $E(B-V)$ values from SED fitting versus β in Figure 13.

Zoom In Zoom Out Reset image size

In Figure 13, we also show the expected values of β given values of $E(B-V)$ using two different empirical calibrations. To calculate the expected β, we use the Calzetti et al. (2000) dust attenuation law to calculate the amount of dust attenuation at rest-frame 1600 Å from $E(B-V)$: ${A}_{1600}=9.97\times E(B-V)$. Then we use the relation between A₁₆₀₀ and β from Meurer et al. (1999, for solar metallicity; ${A}_{1600}=4.43+1.99\beta$) and Castellano et al. (2014, for sub-solar metallicity; ${A}_{1600}={5.32}_{-0.37}^{+0.41}+1.99\beta$) to calculate the expected β. We show the Meurer et al. (1999) relation as a solid line and the Castellano et al. (2014) relation as a dashed line in Figure 13. The scatter of the measured β of our sample is larger than the difference between the two empirical calibrations, but the distribution is roughly consistent with both calibrations. The agreement means that the $E(B-V)$ values derived from SED fitting is not strongly biased on average, but the large scatter also suggests that for each individual galaxy, the dust attenuation is still poorly constrained.