Directly Observing the Galaxies Likely Responsible for Reionization

We carry out a number of tests to determine (a) which of the residual images we should use for object detection, (b) to what extent the wavelet decomposition procedure improves our detection efficiency, and (c) whether we can measure reliable photometry from the residual images. As we need to determine the effect of the wavelet subtraction process on the detection and measurement of galaxies, we do this with a simulation in which we add fake galaxies to the data and attempt to recover them.

We first generate fake point sources with a range of input magnitudes, convolve them with the H₁₆₀-band PSF, and add them to the H₁₆₀-band image in random positions. We then compute the wavelet transform as detailed above, construct models on each scale j = 1...9, and run SExtractor on both the input image and the residual images on each scale. In Figure 4, we compare the results. As the left panel of Figure 4 shows, the wavelet images are useful for detecting sources at $m\gt 26;$ brighter than this, all or part of the Galaxy is removed in the wavelet subtraction process. Although high-redshift galaxies brighter than this limit are rare, we therefore need to use both the original and wavelet-subtracted images for detection to ensure that we find any that may exist. We also need to ensure that we can measure accurate photometry when using the wavelet-subtracted images for detection. The middle panel of Figure 4 shows that photometry measured in the original image for sources detected in the wavelet-subtracted images (scales 3–5) is of comparable accuracy to that of sources detected in the original image (but note that more sources are included for the wavelet-subtracted images at the faint end). If we both detect and measure sources in the wavelet-subtracted images, we find that in the range where sources can be efficiently recovered ($m\gt 26$) the images subtracted on scales 4 and 5 have measurement errors comparable to sources measured in the original image (again, over more sources overall). However, there is a systematic shift toward fainter magnitudes: this is because objects that are partially removed in the wavelet subtraction process or lie close to residuals are still included. Therefore, we only use the wavelet photometry when it is within 1σ of the photometry on the original image (using the wavelet-subtracted image for detection in both cases). This results in accurate recovered photometry on scale j = 4 while maximizing the signal-to-noise ratio by using the wavelet-subtracted image where possible. If we use scale j = 3, the recovered flux is reduced even when measuring in the original image because oversubtraction of galaxy light causes the Kron apertures to be drawn too small. Scale j = 4 is therefore our preferred scale on which to model and subtract the cluster light, and this scale is used throughout the remainder of this paper.

Zoom In Zoom Out Reset image size

We also note that detection in the wavelet-subtracted images is hampered for bright $({H}_{160}\lt 25)$ galaxies as they get removed when we perform the subtraction to smaller scales. Although we do not expect to find any high-redshift galaxies in this regime, we nonetheless choose to combine catalogs from the wavelet-subtracted and original images in order to ensure that we include potential bright or extended sources.

To select high-redshift galaxies from our master catalog, we use a photometric redshift fitting technique. This has the advantage over color selection techniques that it uses all of the available photometry and results in more precise estimates of the error on the resulting redshifts. The redshift fitting code we use is EAZY (Brammer et al. 2008), with an updated set of templates that include contributions from emission lines.

We select high-redshift candidates from the master catalog using the full redshift probability (P(z)) distribution, as well as the best-fit redshift and the detection significance in the WFC3 filters. In our high-redshift selection criteria, we follow previous work on the unlensed CANDELS fields by Finkelstein et al. (2015) in order to ensure that our results are comparable. Therefore, to be selected, a galaxy must meet all of the following criteria:

• 1.  
  The best-fit redshift must be $\gt 4$.
• 2.  
  The integral of the main peak of the P(z) distribution must be $\gt 70 \%$.
• 3.  
  The integral of the P(z) distribution in one of the $z\geqslant 6$ redshift bins (where a bin is defined as an integer ±0.5) must be $\gt 25 \%$.
• 4.  
  The flux must be $\gt 3.5\sigma$ in at least two WFC3 bands, or $\gt 5\sigma$ in the H₁₆₀ band, where we define $1\sigma$ in a position-dependent fashion using the depth maps described in Section 2.

These criteria ensure that the candidates are firm detections with well-defined photometric redshifts and that the high-redshift solution is strongly preferred. To allocate each candidate to a redshift bin ${z}_{\mathrm{bin}}$, we integrate the P(z) distribution in from ${z}_{\mathrm{bin}}-0.5$ to ${z}_{\mathrm{bin}}+0.5$, and the Galaxy is assigned to the ${z}_{\mathrm{bin}}$ with the largest integral. We then keep all sources assigned to ${z}_{\mathrm{bin}}\geqslant 6$.

Once we have catalogs of high-redshift candidate galaxies in each cluster, we go through several steps to clean the sample of potential contaminants.

4.2.1. Visual Inspection

4.2.2. Contamination

We use the same selection method as Finkelstein et al. (2015), in which the contamination rate from faint low-redshift sources was extensively tested and found to be $\lt 10 \%$ with the adopted criteria. We show the stacked P(z) distributions in Figure 5 as an indication of the likely low-redshift contamination. The total fraction of the P(z) at $z\lt 3$ is 13% at $z\sim 6 \%$, 8.8% at $z\sim 7 \%$, and 5.1% at $z\sim 8$.

Zoom In Zoom Out Reset image size

Another potential source of contamination noted by Finkelstein et al. (2015) is that SExtractor cannot efficiently separate stars from galaxies at faint (${J}_{125}\gt 25$) magnitudes. A particular concern at $z\geqslant 6$ is contamination from L and T brown dwarfs, which can have similar colors to high-redshift galaxies. However, the colors of brown dwarfs fall on well-defined sequences and can therefore be easily separated from true galaxies. In Figure 6, we compare the colors of all of our high-redshift galaxy candidates with ${J}_{125}\lt 27$ to the colors of brown dwarfs from the SpeX Prism Spectral Libraries.⁴ We find one source consistent with the brown dwarf colors, indicated in green in Figure 6. This source is the brightest in our sample (${H}_{160}=24.0$) and is unresolved ($\mathrm{FWHM}=0\buildrel{\prime\prime}\over{.} 185$), so we reclassify it as a star and remove it from the sample.

Zoom In Zoom Out Reset image size

An additional potential source of contamination in this work is from globular clusters within the foreground galaxy cluster. At the redshifts of the clusters Abell 2744 and MACS 0416, we find that only very massive ($M\gt {10}^{6}$) and very young ($\lt 1\,$ Gyr) globular clusters would be detectable in these data. We consider the colors of such globular clusters and find that they are very narrowly distributed and inconsistent with the colors of galaxies in our high-redshift sample. Given the highly unlikely presence of young, high-mass globular clusters and their color distributions, we conclude that they are unlikely to be a source of contamination for our sample.

The results of both visual inspections are combined to produce the final catalogs, which contain 167 galaxies: 94 in Abell 2744 and 73 in MACS 0416. The positions are shown in Figure 7, and the catalogs for the $z\sim 6$, 7, and 8 samples are listed in Tables 7–9 in the Appendix. In addition, we find six galaxies at $z\sim 9$, which will be the subject of future work.

Zoom In Zoom Out Reset image size

Crucially, we are able to detect galaxies close to the critical line, the line of theoretically infinite magnification close to which we benefit most from the gravitational lensing effect. Due to the critical line's location near the cluster center, there are a large number of bright foreground galaxies in this region that have hampered previous efforts to detect faint background sources. Subtraction of these foreground galaxies substantially improves the detection significance of the background source. We note that 46 of the 167 galaxies in the sample suffer from unreliable photometry in the wavelet-subtracted images, either due to partial subtraction of the source itself or due to proximity to residuals from the subtraction process. The remaining 121 (or 72% of the sample) have their photometry measured in the wavelet-subtracted images. In order to test whether the sources with bad photometry in the wavelet-subtracted images are disproportionately faint, we compare these two samples using the K-S test and find that they are consistent with being drawn from the same observed magnitude distribution with a probability of 0.992.

Magnification estimates are highly sensitive to substructure within the cluster, so for this work we make use of the full range of possible lens models produced for the Frontier Fields by seven independent teams who use different assumptions and methodologies. The details of the models and the products are available on the MAST archive website⁵ , and we summarize the characteristics of the different models in Table 2. The primary difference is that some models assume that the mass in the cluster is traced by the luminous galaxies, while others operate without that assumption. The light-traces-mass assumption has been well tested and proven to improve the accuracy of lens models, but at the expense of flexibility (Zitrin et al. 2010, 2011); models that do not assume that light traces mass allow a wider range of possible mass distributions. Both of the clusters considered here are undergoing mergers, so we might expect more deviations from light tracing mass than usual (Wang et al. 2015).

Table 2.  Details of the Lens Models Used

Team	Version(s)	Version(s)	Area	Resolution	Light Traces Mass?	${n}_{\mathrm{models}}$	References
	Pre-HFF	Post-HFF
Bradac	1	2a, 3b	35 × 3.5	005	No	100	Bradac et al. (2004)
							Bradac et al. (2009)
CATS	1	2, 3, 3.1	53 × 61	02	Yes, plus dark matter	200	Jullo & Kneib (2009)
							Jauzac et al. (2012)
							Jauzac et al. (2014)
							Richard et al. (2014)
							Jauzac et al. (2015a)
glafic	1a	3	27 × 27	03	Yes, plus dark matter	100	Oguri (2010)
							Ishigaki et al. (2015)
							Kawamata et al. (2016)
Merten	1		25' × 25'	83	No	250	Merten et al. (2009)
							Merten et al. (2011)
Williams	1	3, 3.1b	23 × 23	03	No	30	Liesenborgs et al. (2006)
							Jauzac et al. (2014)
							Mohammed et al. (2014)
							Grillo et al. (2015)
Sharon	2	3	34 × 34	003	Yes, plus dark matter	100	Jullo et al. (2007)
							Johnson et al. (2014)
Zitrin-NFW	1	3	30 × 30	006	Yes, plus dark matter	100	Zitrin et al. (2009)
							Zitrin et al. (2013)
Zitrin-LTM	1	3	30 × 30	006	Yes	100	Zitrin et al. (2009)
							Zitrin et al. (2013)

Notes. ${n}_{\mathrm{models}}$ is the number of alternative models provided by each team in addition to their best-fit model.

^aA2744 only. ^bM0416 only.

Download table as:  ASCIITypeset image

The magnification, μ, of a galaxy depends on the convergence (κ) and shear (γ) as

$$\begin{eqnarray}&&\mu =\displaystyle \frac{1}{{(1-\kappa )}^{2}-{\gamma }^{2}}.\end{eqnarray} \tag{ 3 }$$

Both κ and γ scale with the distance ratio ${D}_{\mathrm{LS}}/{D}_{S}$, where D_LS is the angular diameter distance between the source and the lensing cluster and D_S is the angular diameter distance between the observer and the source.

The differential magnification across a galaxy image can be a strong effect close to the critical line where the magnification gradient is high. In order to account for this, we use a flux-weighted magnification calculated from the H₁₆₀-band image for each galaxy, using the best-fit normalized κ and γ maps provided by each lens modeling team.

In addition to their best-fit maps, each lens modeling team has made available a range of κ and γ maps from alternate models to enable the estimation of errors. While the redshifts of the clusters are well known, the redshifts of the background sources are subject to the errors on their photometric redshifts. When calculating errors on the estimated magnification factors, we must therefore incorporate the errors associated with κ, γ, and ${D}_{\mathrm{LS}}/{D}_{S}$. We estimate errors on κ and γ by using the full range of lens models provided by each team. The error on ${D}_{\mathrm{LS}}/{D}_{S}$ is estimated from the 68% uncertainty range on the photometric redshift of the source (we assume that the error on the distance to the cluster is negligible).

To estimate the uncertainty on the overall magnification factors, we use a bootstrap method. We create 10,000 realizations of the magnification of each galaxy based on randomly drawn values of κ and γ from within the range of alternate lens models provided by each team and a randomly selected redshift from within the probability distribution function (PDF) for the Galaxy. We find that the uncertainty on the photometric redshift dominates over the lens model uncertainties except in cases where the magnification is high ($\mu \gtrsim 10$), where the precise position of the critical line can cause significant differences in the magnification.

In Tables 7–9 in the Appendix we list the median magnification across the full range of models that were updated after the Frontier Fields data were obtained (see Table 2 for the model version numbers). The uncertainties listed encompass the central 68% interval over all of the models, in order to separate the model uncertainty from the photometric redshift. We use this median value for our fiducial luminosity function, as it excludes the extreme magnifications, which are likely unphysical. However, in order to examine the systematic uncertainties from lens modeling, we also discuss the effects of using each lens model individually in Section 5.3.

In Figure 8, we show images of the two intrinsically faintest sources (${M}_{\mathrm{UV}}=-12.4$ and −14.2) in the sample, both before and after wavelet subtraction.

Zoom In Zoom Out Reset image size

The most comprehensive previous study of the $z=6\mbox{--}8$ luminosity function from the Frontier Fields to date by Atek et al. (2015a) contains 253 galaxies, of which 88 are selected in the two fields Abell 2744 and MACS 0416 (a further 74 are from these clusters' parallel fields; objects in these fields are included with our unlensed comparison sample from Finkelstein et al. 2015). They combine their $z\sim 6$ and $z\sim 7$ samples into one $z\sim 7$ luminosity function, so in order to compare this sample with ours, we use their photometric redshifts (supplied via private communication) to split the 88 galaxies from the first two cluster fields into the same redshift bins ($z\sim 6,7,8$) used in this work. The results are shown in Figure 9; it is likely that some galaxies scatter between bins, due to different magnification estimates, but we find $\sim 2\times$ more galaxies and probe systematically fainter than this previous work.

Zoom In Zoom Out Reset image size

To test whether the larger sample size is caused by more efficient detection or differences in the selection criteria, we reselect our high-redshift catalog using the wavelet subtraction method, but with the selection criteria of Finkelstein et al. (2015). Briefly, these comprise color–color cuts to select a strong break, a nondetection in the optical bands, and a signal-to-noise cut in the infrared. With these criteria—and after carrying out visual inspection to remove spurious sources—we find a total of 204 galaxies between the two cluster fields; of these, 119 are in Abell 2744 and 85 in MACS 0416. There are seven galaxies in our sample that would not be included using these criteria (all brighter than ${M}_{\mathrm{UV}}=-15$), leaving 44 new sources that are not in our sample. In all cases, these new sources do not meet the photometric redshift quality criteria we use in our sample selection. This implies that our selection criteria are comparatively strict, but that we find more sources overall, due to the subtraction of the ICL.

In order to compute the effective volume for the luminosity function, we carry out simulations to determine the recovery fraction of galaxies as a function of magnitude, position, size, and color using 100,00 CPU hr provided by the Texas Advanced Computing Center (TACC). The method is similar to that employed in Section 3.1 to determine the reliability of the wavelet decomposition process, but we now allow a range of intrinsic sizes and colors for the fake galaxies, incorporate the effect of lensing on their shapes, attempt to recover them using the full range of 22 detection images as described in Section 4, and finally carry out photo-z fitting on the recovered galaxies in the same manner as for the real galaxies and incorporate the effects of our photo-z selection.

The first step is to create the fake galaxies. For computational efficiency, we need to simulate multiple galaxies in each iteration, without adding artificial confusion. We therefore carry out preliminary tests to determine the optimal number of galaxies to add to each image. The results are unchanged when up to 500 galaxies are added at a time, so we conservatively opt for 250. For the sizes, we use a normal distribution of half-light radii with a peak at 0.5 kpc (1.5 pixels) and hard cutoffs at 0 and 5 kpc (∼14.5 pixels at $z\sim 6;$ Kawamata et al. 2014; Shibuya et al. 2015). We choose Sérsic indices (n) from a lognormal distribution in the range $1\lt n\lt 5$ with a peak at $n\sim 1.5$ (for disk-like morphologies). The axial ratio is also lognormal with a peak at 1.8, and the position angle on the sky is selected from a uniform random distribution between 0° and 360°. To obtain a range of colors, we define ranges of redshifts, ages, $E(B-V)$, and metallicities. The redshift is selected from a uniform distribution in the range $z=5\mbox{--}11$, extending above and below the redshift range of interest. In each realization of the simulation, we use the same redshift for all sources to simplify the transformation to the image plane. We assign a dust attenuation factor $E(B-V)$ from a normal distribution centered on 0.1 with $\sigma =0.15$, restricted to the range 0–0.5. The age is a lognormal distribution peaked at 7 Myr and limited to less than the age of the universe at the specified redshift. We also use a lognormal distribution in metallicity peaked at 0.2 ${Z}_{\odot }$. These parameters together combine to provide a distribution of UV slopes β, designed to match the distribution we expect in galaxies at these redshifts (Finkelstein et al. 2012; Bouwens et al. 2014). Given these distributions, we use the stellar population synthesis models of Bruzual & Charlot (2003) to calculate the colors of the galaxies. We then select the H-band magnitudes from a uniform distribution between 22 and 38 and use the colors derived from the stellar population models to obtain the magnitudes in the other six filters.

We use the magnitudes, sizes, and shapes selected above to generate a 101 × 101 pixel postage stamp image of each source with galfit. We assign positions to the fake sources uniformly in the image plane and then trace this position back to the source plane using the lens model. This method ensures complete coverage of the image plane, whereas uniformly sampling the source plane causes the more magnified regions to be undersampled. The fake galaxy images are each added to a blank image of the source plane of the field at the appropriate positions and then lensed back to the image plane using lenstool (Kneib 1993; Jullo et al. 2007; Jullo & Kneib 2009) with the latest lens models (Jauzac et al. 2015a, 2015b). This enables us to account for some lensing shear in our recovery fractions. Due to the computational time required for these simulations, we are limited to using a single lens model. Since the completeness is calculated in the image plane, the only uncertainty this introduces is in the difference between lensing shears in different models. This is not expected to be a dominant effect, but is a limitation of these time-intensive simulations. Once the galaxies have been lensed, we convolve the image with the H₁₆₀-band PSF and use SExtractor to recover them and measure the image-plane half-light radii and magnitudes. These, together with all of the input parameters above and the positions in the image and source planes, are saved into an input array. We then add the images consisting of fake, lensed galaxies to the real images in each filter. We then carry out the wavelet decomposition procedure described in Section 3, construct the full range of 22 detection images listed in Section 4, and construct new catalogs using the same method employed for the real data. From these catalogs, we match to the input fake source catalog using a 02 matching radius. We then measure the photometric redshifts of these recovered sources using EAZY in the same manner employed for the real catalog and apply the same cuts as described in Section 4.1 to determine which galaxies we would have selected in our high-redshift sample.

We repeat the above process 2000 times, resulting in 500,000 fake galaxies in each cluster. This number is required because the completeness is highly position dependent in strong-lensing fields. We then collate the results and calculate the recovery fraction as a function of position and observed magnitude.

The effective volume, V_m, for a given intrinsic magnitude m is given by

$$\begin{eqnarray}&&{V}_{m}=\displaystyle \sum _{i=0}^{N}\displaystyle \frac{v}{{\mu }_{i}}{f}_{i-{\mu }_{i,\mathrm{mag}}},\end{eqnarray} \tag{ 6 }$$

where N is the number of pixels in the image, v is the volume of one unlensed pixel, ${\mu }_{i}$ is the linear magnification factor of pixel i, ${f}_{i-{\mu }_{i,\mathrm{mag}}}$ is the completeness fraction in pixel i at image-plane magnitude $i-{\mu }_{i,\mathrm{mag}}$, and ${\mu }_{i,\mathrm{mag}}$ is the magnification in pixel i expressed in magnitudes. The resulting total volume as a function of magnitude for each cluster, lens model, and redshift bin is shown in Figure 10. Note that due to computational run time, we do not use the shear from different lens models, but we do account for the different magnifications when calculating the effective volumes.

Zoom In Zoom Out Reset image size

Calculating the volume in the image plane in this way means that regions that are multiply imaged are counted multiple times, with the appropriate recovery fraction for each of the multiple images. By double-counting both the multiply imaged galaxies and their respective volumes, we can eliminate errors in the luminosity function due to misidentified multiple images.

Using the Galaxy sample from Section 4 and effective volumes based on the completeness simulations described in Section 5.1, we can produce the luminosity functions at $6\lt z\lt 8$. The resulting luminosity functions are shown in Figure 11. The error bars in Figure 11 are calculated from a combination of the Poisson uncertainty and the uncertainty on the intrinsic magnitudes. To estimate the uncertainty on the number of galaxies in the bin, we carry out 10,000 realizations of the luminosity function in which each galaxy has a new intrinsic magnitude based on uncertainties in its photometry, magnification (using the bootstrap method described in Section 4.4), and photometric redshift. The final uncertainty on the number density is the 68% central confidence interval of the number of galaxies in each bin over all realizations. This is added in quadrature to the Poisson uncertainty, for which we use Equation (33.59) of Nakamura et al. (2010) appropriate for small numbers.

Zoom In Zoom Out Reset image size

In order to fit Schechter function parameters to our luminosity functions, we employ a Markov Chain Monte Carlo (MCMC) algorithm to search the parameter space. Since our data constrain only the faint end, we include unlensed data from the CANDELS fields that use the same selection method (Finkelstein et al. 2015) in order to add constraints at the bright end. In absolute magnitude bins of ${\rm{\Delta }}M=0.5$ and for each field individually, we calculate the number of galaxies expected for the given Schechter parameters M*, ϕ*, and α. Within each subfield we add a randomly generated value for cosmic variance based on a normal distribution with the width of the expected fractional cosmic variance in that field (Robertson et al. 2014). This is compared to the number of observed galaxies in each bin, which is perturbed on each iteration within the uncertainties derived above. It is also important to account for the systematic effect of errors in magnitudes, which can cause an artificially steepened observed slope in the luminosity function due to Eddington bias. Since the intrinsic slope of the luminosity function is steep, large magnitude errors mean that galaxies preferentially scatter into brighter bins. This has been shown to be especially important at the faint end, where low signal-to-noise ratio means that photometric errors are larger than at the bright end (Grazian et al. 2015). This effect becomes even more important in lensing fields where the faint end consists of galaxies with high magnifications, which also have high uncertainties.

To account for the Eddington bias in our sample, we construct a magnitude PDF for each galaxy. This is derived by selecting 10,000 realizations of the magnitude of each galaxy by drawing randomly from within the photometric and magnification errors and from the photo-z PDF. We combine all of the magnitude distribution functions within each subfield to produce a PDF $P({M}_{i},{M}_{j})$ that a galaxy with magnitude M_i has magnitude M_j consistent with the uncertainties in its photo-z, photometry, and magnification. These PDFs are narrow at the bright end, where the photometry is more certain and there is little or no lensing magnification, but become much broader at the faint end, where all of these uncertainties are higher.

To calculate the expected luminosity function in each subfield f and magnitude bin M_i, we have

$$\begin{eqnarray}&&{\phi }_{i}=\displaystyle \sum _{j=0}^{N}{\phi }_{j,\mathrm{int}}(1+{\mathrm{CV}}_{j})P(j,i),\end{eqnarray} \tag{ 7 }$$

where CV_j is the cosmic variance estimate in magnitude bin M_i, drawn from a random normal distribution with the width of the estimate of fractional cosmic variance from Robertson et al. (2014), and ${\phi }_{j,\mathrm{int}}$ is the intrinsic Schechter function at magnitude j.

For each combination of Schechter parameters M*, ϕ*, and α, we calculate the goodness-of-fit statistic

$$\begin{eqnarray}&&{{ \mathcal C }}^{2}(\phi )=-2\mathrm{ln}{ \mathcal L }(\phi ),\end{eqnarray} \tag{ 8 }$$

where ${ \mathcal L }$ is the likelihood that the number of galaxies observed in that field and magnitude bin matches the number expected according to Equation (7). The final goodness-of-fit ${{ \mathcal C }}^{2}$ is the sum over all fields and magnitude bins at a given redshift.

We use an IDL implementation of an affine-invariant ensemble MCMC sampler (Foreman-Mackey et al. 2013; Finkelstein 2015) to search the parameter space. For each redshift, we use 10³ independent chains of 10⁵ steps each, computing the likelihood of the given Schechter parameters at each step. Using independent chains prevents the fit from being trapped by local minima in the parameter space. To compute the final best fit, we join all of the chains together to give 10⁸ Schechter function parameters for each redshift. The best-fit parameters listed in Table 3 are the median of this distribution, with errors covering the central 68% of values.

Table 3.  Best-fitting Schechter Function Values

Redshift	M*	$\mathrm{log}\phi *$	α
6	$-{20.825}_{-0.040}^{+0.051}$ (stat.)${}_{-0.003}^{+0.004}$ (sys.)	$-{3.647}_{-0.033}^{+0.037}$ (stat.)${}_{-0.004}^{+0.002}$ (sys.)	$-{2.10}_{-0.03}^{+0.03}$ (stat.)${}_{-0.01}^{+0.02}$ (sys.)
7	$-{20.802}_{-0.050}^{+0.055}$ (stat.)${}_{-0.003}^{+0.007}$ (sys.)	$-{3.673}_{-0.043}^{+0.046}$ (stat.)${}_{-0.004}^{+0.001}$ (sys.)	$-{2.06}_{-0.04}^{+0.04}$ (stat.)${}_{-0.01}^{+0.01}$ (sys.)
8	$-{20.723}_{-0.142}^{+0.181}$ (stat.)${}_{-0.018}^{+0.014}$ (sys.)	$-{3.768}_{-0.131}^{+0.129}$ (stat.)${}_{-0.004}^{+0.022}$ (sys.)	$-{2.01}_{-0.08}^{+0.08}$ (stat.)${}_{-0.02}^{+0.02}$ (sys.)

Note. Parameters are fit using an MCMC method as described in the text, with a combination of the lensed HFF sample and unlensed CANDELS, HUDF, and HFF Year 1 parallel field sample from Finkelstein et al. (2015). Statistical errors are the 68% central range of the MCMC samples, and systematic errors are based on the range of available lens models.

Download table as:  ASCIITypeset image

It should be noted that the two faintest bins in our $z\sim 6$ luminosity function that contain data are based on just one galaxy each. These are highly magnified and have correspondingly large uncertainties in their magnifications; different lens models give a magnification factor for the faintest source of anywhere between $23\times$ (Zitrin_LTM_Gauss) and $150\times$ (Bradac v2, Sharon v2), with a median of $100\times$. The MCMC code used to fit the Schechter parameters should be more strongly constrained by bins containing large numbers of galaxies, but to ensure that these two bins are not skewing the results, we fit new parameters at $z\sim 6$ with these galaxies entirely excluded. We obtain exactly the same value for the faint-end slope, with error bars larger by 6%.

For our fiducial luminosity functions, we use the median magnification for each galaxy in the sample. Here, we consider the systematic uncertainty due to the lens modeling by using each available lens model individually and looking at the full range of results of the best-fitting Schechter parameters.

We use lens models produced both before and after the Frontier Fields data were acquired. The specific versions of each team's model used are listed in Tables 4–6, and the resulting luminosity functions for each redshift bin are shown in Figure 12. There is scatter between the luminosity functions produced with different lens models, as they each result in different magnification estimates for the sample galaxies. However, they all show a steeply rising faint-end slope.

Zoom In Zoom Out Reset image size

Table 4.  Schechter Function Fit Parameters from Individual Lens Models at z = 6

Model	Pre-HFF			Post-HFF
	M*	$\mathrm{log}\phi *$	α	M*	$\mathrm{log}\phi *$	α
Bradac (v1, 2a/3b)	$-{20.830}_{-0.038}^{+0.053}$	$-{3.650}_{-0.033}^{+0.041}$	$-{2.10}_{-0.04}^{+0.03}$	$-{20.825}_{-0.037}^{+0.046}$	$-{3.644}_{-0.033}^{+0.037}$	$-{2.09}_{-0.03}^{+0.03}$
CATS (v1, 3)	$-{20.822}_{-0.041}^{+0.050}$	$-{3.644}_{-0.036}^{+0.038}$	$-{2.10}_{-0.04}^{+0.04}$	$-{20.827}_{-0.033}^{+0.049}$	$-{3.648}_{-0.031}^{+0.035}$	$-{2.10}_{-0.03}^{+0.03}$
glafic (v1a, 3)	$-{20.827}_{-0.042}^{+0.051}$	$-{3.644}_{-0.034}^{+0.037}$	$-{2.09}_{-0.04}^{+0.03}$	$-{20.821}_{-0.041}^{+0.048}$	$-{3.645}_{-0.034}^{+0.035}$	$-{2.10}_{-0.03}^{+0.03}$
Merten (v1)	$-{20.816}_{-0.037}^{+0.043}$	$-{3.638}_{-0.035}^{+0.033}$	$-{2.09}_{-0.03}^{+0.03}$	⋯	⋯	⋯
Sharon (v2, 3)	$-{20.822}_{-0.040}^{+0.050}$	$-{3.643}_{-0.033}^{+0.036}$	$-{2.10}_{-0.03}^{+0.03}$	$-{20.826}_{-0.038}^{+0.049}$	$-{3.649}_{-0.033}^{+0.038}$	$-{2.10}_{-0.03}^{+0.03}$
Williams (v1, 3)	$-{20.820}_{-0.037}^{+0.046}$	$-{3.645}_{-0.034}^{+0.038}$	$-{2.10}_{-0.03}^{+0.03}$	$-{20.820}_{-0.042}^{+0.051}$	$-{3.642}_{-0.037}^{+0.037}$	$-{2.10}_{-0.03}^{+0.03}$
Zitrin-LTM (v1)	$-{20.824}_{-0.039}^{+0.050}$	$-{3.647}_{-0.033}^{+0.040}$	$-{2.11}_{-0.03}^{+0.04}$	⋯	⋯	⋯
Zitrin-NFW (v1, 3)	$-{20.861}_{-0.066}^{+0.076}$	$-{3.631}_{-0.047}^{+0.053}$	$-{2.03}_{-0.06}^{+0.07}$	$-{20.861}_{-0.079}^{+0.081}$	$-{3.631}_{-0.052}^{+0.056}$	$-{2.02}_{-0.07}^{+0.07}$
Zitrin-LTM-Gauss (v1b, 3)	$-{20.829}_{-0.040}^{+0.055}$	$-{3.642}_{-0.036}^{+0.041}$	$-{2.08}_{-0.04}^{+0.04}$	$-{20.821}_{-0.037}^{+0.047}$	$-{3.644}_{-0.036}^{+0.036}$	$-{2.10}_{-0.03}^{+0.03}$
Mean	−20.828 ± 0.013	−3.643 ± 0.005	−2.09 ± 0.02	−20.829 ± 0.014	−3.643 ± 0.006	−2.09 ± 0.03
Median	$-{20.824}_{-0.005}^{+0.004}$	$-{3.644}_{-0.003}^{+0.005}$	$-{2.10}_{-0.01}^{+0.01}$	$-{20.825}_{-0.003}^{+0.004}$	$-{3.644}_{-0.004}^{+0.002}$	$-{2.10}_{-0.01}^{+0.02}$

Notes. Best-fitting Schechter function parameters for each lens model, for the models produced by each team both before and after the Frontier Fields data. The model versions shown are those used for the (pre, post)-HFF fits. All uncertainties are the 68% central confidence interval, except for the error on the mean, which is given as the standard deviation of all of the models.

^aModel exists for A2744 only. ^bModel exists for M0416 only.

Download table as:  ASCIITypeset image

Table 5.  Schechter Function Fit Parameters from Individual Lens Models at z = 7

Model	Pre-HFF			Post-HFF
	M*	$\mathrm{log}\phi *$	α	M*	$\mathrm{log}\phi *$	α
Bradac (v1, 2a/3b)	$-{20.796}_{-0.045}^{+0.050}$	$-{3.673}_{-0.042}^{+0.047}$	$-{2.08}_{-0.04}^{+0.04}$	$-{20.799}_{-0.046}^{+0.050}$	$-{3.671}_{-0.042}^{+0.043}$	$-{2.08}_{-0.04}^{+0.04}$
CATS (v1, 3)	$-{20.798}_{-0.049}^{+0.053}$	$-{3.671}_{-0.047}^{+0.046}$	$-{2.07}_{-0.04}^{+0.04}$	$-{20.815}_{-0.056}^{+0.063}$	$-{3.676}_{-0.052}^{+0.052}$	$-{2.05}_{-0.05}^{+0.05}$
glafic (v1a, 3)	$-{20.799}_{-0.051}^{+0.053}$	$-{3.669}_{-0.043}^{+0.045}$	$-{2.08}_{-0.05}^{+0.05}$	$-{20.806}_{-0.052}^{+0.060}$	$-{3.675}_{-0.048}^{+0.049}$	$-{2.05}_{-0.05}^{+0.05}$
Merten (v1)	$-{20.805}_{-0.054}^{+0.061}$	$-{3.673}_{-0.049}^{+0.051}$	$-{2.05}_{-0.05}^{+0.05}$	⋯	⋯	⋯
Sharon (v2, 3)	$-{20.799}_{-0.047}^{+0.053}$	$-{3.669}_{-0.044}^{+0.046}$	$-{2.07}_{-0.05}^{+0.04}$	$-{20.799}_{-0.049}^{+0.057}$	$-{3.670}_{-0.049}^{+0.047}$	$-{2.06}_{-0.04}^{+0.05}$
Williams (v1, 3)	$-{20.802}_{-0.050}^{+0.057}$	$-{3.674}_{-0.046}^{+0.046}$	$-{2.07}_{-0.04}^{+0.04}$	$-{20.797}_{-0.047}^{+0.059}$	$-{3.671}_{-0.047}^{+0.050}$	$-{2.05}_{-0.04}^{+0.04}$
Zitrin-LTM (v1)	$-{20.800}_{-0.048}^{+0.055}$	$-{3.671}_{-0.046}^{+0.047}$	$-{2.07}_{-0.05}^{+0.05}$	⋯	⋯	⋯
Zitrin-NFW (v1, 3)	$-{20.814}_{-0.062}^{+0.067}$	$-{3.674}_{-0.053}^{+0.051}$	$-{2.04}_{-0.07}^{+0.07}$	$-{20.809}_{-0.058}^{+0.066}$	$-{3.674}_{-0.058}^{+0.051}$	$-{2.05}_{-0.06}^{+0.06}$
Zitrin-LTM-Gauss (v1b, 3)	$-{20.809}_{-0.052}^{+0.066}$	$-{3.676}_{-0.050}^{+0.053}$	$-{2.06}_{-0.06}^{+0.06}$	$-{20.806}_{-0.057}^{+0.061}$	$-{3.668}_{-0.050}^{+0.051}$	$-{2.05}_{-0.05}^{+0.05}$
Mean	−20.802 ± 0.006	−3.673 ± 0.002	−2.07 ± 0.01	−20.804 ± 0.006	−3.672 ± 0.003	−2.05 ± 0.01
Median	$-{20.800}_{-0.008}^{+0.002}$	$-{3.673}_{-0.002}^{+0.003}$	$-{2.07}_{-0.01}^{+0.01}$	$-{20.806}_{-0.003}^{+0.007}$	$-{3.671}_{-0.004}^{+0.001}$	$-{2.05}_{-0.01}^{+0.01}$

Note. See notes to Table 4.

Download table as:  ASCIITypeset image

Table 6.  Schechter Function Fit Parameters from Individual Lens Models at z = 8

Model	Pre-HFF			Post-HFF
	M*	$\mathrm{log}\phi *$	α	M*	$\mathrm{log}\phi *$	α
Bradac (v1, 2a/3b)	$-{20.719}_{-0.140}^{+0.178}$	$-{3.780}_{-0.141}^{+0.139}$	$-{2.05}_{-0.09}^{+0.09}$	$-{20.726}_{-0.146}^{+0.190}$	$-{3.788}_{-0.152}^{+0.146}$	$-{2.08}_{-0.09}^{+0.09}$
CATS (v1, 3)	$-{20.708}_{-0.146}^{+0.180}$	$-{3.764}_{-0.129}^{+0.123}$	$-{1.99}_{-0.09}^{+0.09}$	$-{20.712}_{-0.158}^{+0.193}$	$-{3.756}_{-0.133}^{+0.129}$	$-{1.98}_{-0.09}^{+0.08}$
glafic (v1a, 3)	$-{20.748}_{-0.135}^{+0.189}$	$-{3.808}_{-0.144}^{+0.156}$	$-{2.09}_{-0.10}^{+0.09}$	$-{20.711}_{-0.149}^{+0.186}$	$-{3.764}_{-0.131}^{+0.125}$	$-{2.00}_{-0.08}^{+0.08}$
Merten (v1)	$-{20.722}_{-0.159}^{+0.183}$	$-{3.767}_{-0.132}^{+0.128}$	$-{1.99}_{-0.08}^{+0.08}$	⋯	⋯	⋯
Sharon (v2, 3)	$-{20.685}_{-0.168}^{+0.199}$	$-{3.745}_{-0.137}^{+0.124}$	$-{1.98}_{-0.10}^{+0.09}$	$-{20.743}_{-0.155}^{+0.177}$	$-{3.785}_{-0.136}^{+0.135}$	$-{1.99}_{-0.08}^{+0.08}$
Williams (v1, 3)	$-{20.735}_{-0.143}^{+0.180}$	$-{3.790}_{-0.140}^{+0.139}$	$-{2.04}_{-0.08}^{+0.08}$	$-{20.718}_{-0.147}^{+0.186}$	$-{3.774}_{-0.141}^{+0.136}$	$-{2.02}_{-0.08}^{+0.08}$
Zitrin-LTM (v1)	$-{20.720}_{-0.150}^{+0.189}$	$-{3.780}_{-0.141}^{+0.137}$	$-{2.02}_{-0.09}^{+0.09}$	⋯	⋯	⋯
Zitrin-NFW (v1, 3)	$-{20.745}_{-0.199}^{+0.257}$	$-{3.826}_{-0.184}^{+0.190}$	$-{2.02}_{-0.14}^{+0.14}$	$-{20.754}_{-0.215}^{+0.260}$	$-{3.830}_{-0.193}^{+0.191}$	$-{2.01}_{-0.14}^{+0.15}$
Zitrin-LTM-Gauss (v1b, 3)	$-{20.752}_{-0.180}^{+0.233}$	$-{3.823}_{-0.171}^{+0.180}$	$-{2.03}_{-0.10}^{+0.11}$	$-{20.735}_{-0.150}^{+0.185}$	$-{3.786}_{-0.139}^{+0.136}$	$-{2.03}_{-0.08}^{+0.08}$
Mean	−20.726 ± 0.022	−3.787 ± 0.028	−2.02 ± 0.04	−20.728 ± 0.016	−3.783 ± 0.024	−2.02 ± 0.03
Median	$-{20.722}_{-0.026}^{+0.011}$	$-{3.780}_{-0.039}^{+0.015}$	$-{2.02}_{-0.03}^{+0.03}$	$-{20.726}_{-0.018}^{+0.014}$	$-{3.785}_{-0.004}^{+0.022}$	$-{2.01}_{-0.02}^{+0.02}$

Note. See notes to Table 4.

Download table as:  ASCIITypeset image

We carry out the full MCMC Schechter parameter fit to each lens model individually in the same manner as for the fiducial results, and the resulting best-fitting parameters are listed in Tables 4–6.

The different lens models do not have a strong effect on the fits of M* and $\phi *$ at z = 6 or 7, as the bright end is constrained by the unlensed CANDELS data. There are differences between the models at $z\sim 8$, where the overdensity of bright galaxies in Abell 2744 affects the fit.

We expect the lens model uncertainties to affect the measurement of α, which is primarily constrained by lensed galaxies, and we do find a wider range of values of α by looking at the different lens models. However, every lens model—those produced both before and after the HFF data—is consistent with the fiducial result to within $1\sigma$. We find no systematic difference between the models that use the light-traces-mass assumption and those that do not. At $z\sim 6$, the differences between the pre- and post-HFF models are statistically insignificant, and all of the models are in agreement within 1σ. At $z\sim 7$ the post-HFF versions of the models produce systematically shallower faint-end slopes than the originals, but the differences here and at $z\sim 8$ are again all within 1σ of the original models.

We conclude that although the different lens models indicate large differences in the magnifications of individual galaxies, the effect on the overall luminosity function is minimal. We caution, however, that although the different lens models cover different assumptions and methodologies, they all use the same input data, so this comparison is not indicative of the full errors in the lens modeling. However, it does provide an estimate of the uncertainties resulting from the lens modeling methodology, and we incorporate this as a systematic uncertainty in our fiducial Schechter function fits in Table 3.