THE EVOLUTION OF THE GALAXY REST-FRAME ULTRAVIOLET LUMINOSITY FUNCTION OVER THE FIRST TWO BILLION YEARS

Studying galaxies in the early universe requires extremely deep imaging, necessitating space-based data. Additionally, we need to combine deep studies over small areas with larger-area surveys with shallower limiting magnitudes to probe a large dynamic range in luminosities. Our study used imaging data from a number of surveys covering both the northern and southern fields from the Great Observatories Origins Deep Survey (Giavalisco et al. 2004), with both the HST and the Spitzer Space Telescope.

The deepest imaging comes from three surveys of the HUDF: the original HUDF survey, which obtained optical imaging with the Advanced Camera for Surveys (ACS; Beckwith et al. 2006), and the more recent HUDF09 (PI Illingworth; e.g., Bouwens et al. 2010a; Oesch et al. 2010b) and UDF12 surveys (PI Ellis; Ellis et al. 2013; Koekemoer et al. 2013), which obtained near-infrared imaging with the Wide Field Camera 3 (WFC3). The full HST data set over the HUDF comprises imaging in eight bands, F435W, F606W, F775W, and F850LP with ACS; and F105W, F125W, F140W, and F160W with WFC3 (hereafter ${B}_{435},$ V₆₀₆, i₇₇₅, z₈₅₀, Y₁₀₅, J₁₂₅, JH₁₄₀, and H₁₆₀, respectively), which cover an area of ∼5 arcmin². The HUDF09 survey also obtained deep WFC3 imaging over two similarly sized flanking fields, first observed with ACS in the UDF05 survey (PI Stiavelli; Oesch et al. 2007), referred to as the HUDF09-01 and HUDF09-02 fields (Bouwens et al. 2011b). These fields each have imaging in the V₆₀₆, i₇₇₅, z₈₅₀, Y₁₀₅, J₁₂₅, and H₁₆₀ bands.

The majority of our candidate galaxy sample comes from CANDELS (PIs Faber and Ferguson; Grogin et al. 2011; Koekemoer et al. 2011). The largest HST project ever, CANDELS comprises 902 orbits over five extragalactic deep fields, including the two GOODS fields (Giavalisco et al. 2004). CANDELS, which finished in 2013 August, is composed of a deep and a wide survey. The deep survey covers the central ∼50% of each of the two GOODS fields, while the wide survey covers the remainder of the GOODS-N field, as well as the southern ∼25% of the GOODS-S field to depths ∼1 mag shallower than the deep survey (the wide survey also covers three additional fields not used in this study; see Section 6.4.1 and Figure 16). We use ACS imaging from the original GOODS survey in the B₄₃₅, V₆₀₆, i₇₇₅, and z₈₅₀ bands. The most recent ACS mosaics in these fields were used; in GOODS-S this includes all ACS imaging in that field prior to the ACS repair on Servicing Mission 4 in 2009, and the GOODS-N field includes all ACS imaging from the GOODS survey (CANDELS internal team release versions 3 and 2, respectively). The CANDELS imaging in both the deep and wide regions of the GOODS fields includes the Y₁₀₅, J₁₂₅, and H₁₆₀ bands. We add to our GOODS-S data set imaging over the northern ∼25% of the GOODS-S field from the WFC3 Science Oversight Committee's Early Release Science (ERS) program (PI O'Connell; Windhorst et al. 2011), which also includes J₁₂₅ and H₁₆₀ imaging, as well as the F098M (hereafter Y₀₉₈) band. Unless otherwise noted, hereafter we will refer to Y₀₉₈ and Y₁₀₅ together as the Y-band (both filters probe observed 1 μm light, but the Y₀₉₈ filter is narrower and thus has a higher spectral resolution).

Finally, we complete our data set with the recently obtained deep HST observations near the galaxy clusters Abell 2744 and MACS J0416.1-2403 (hereafter MACS0416) from the HFF program (PI Lotz). For this study, we use only the parallel (unlensed) fields. Both fields have been observed in the B₄₃₅, V₆₀₆, I₈₁₄, Y₁₀₅, J₁₂₅, JH₁₄₀, and H₁₆₀ bands. We use these data to complement our candidate galaxy samples at z = 5, 6, 7, and 8 (excluding z = 4 due to the reduced number of optical bands).

In parallel to the primary WFC3 observations, CANDELS obtained extremely deep imaging in the F814W band (hereafter I₈₁₄) in both of the GOODS fields. As these data were obtained recently, they suffer from poor charge transfer efficiency. Although algorithms have been devised to correct for this (Anderson & Bedin 2010), we do not include the CANDELS I₈₁₄ photometry in the initial photometric redshift fitting (although we do explore its inclusion in Section 3.6), because the CANDELS fields have imaging in both the i₇₇₅ and z₈₅₀ bands. However, we did use these very deep data during our visual inspection step, which was highly useful at z = 8, where true z = 8 galaxies should be completely undetected in the I₈₁₄-band. In the HFF parallel fields, where the I₈₁₄ band is the only imaging covering the red end of the optical, we used these data in the full analysis.

The description of the CANDELS HST imaging reduction is available from Koekemoer et al. (2011). These reduction steps were also followed for the ERS, HUDF (Koekemoer et al. 2013), and HFF data we use here. We use imaging mosaics with 006 pixels, and make use of their associated weight and rms maps. The combined imaging data set covers an area of 301.2 arcmin², with 5σ limiting magnitudes in the H₁₆₀ band ranging from 27.4 to 29.7 mag (measured in 04 diameter apertures). These data sets are summarized in Table 1.

Table 1.  Summary of Data—Limiting Magnitudes

Field	Area	B435	V606	i775	I814	z850	${Y}_{098/105}$	J125	JH140	H160
	(arcmin2)	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)
GOODS-S Deep	61.6	28.2	28.6	27.9	28.1	27.8	28.2	28.1	⋯	27.9
GOODS-S ERS	41.4	28.2	28.5	27.9	27.9	27.6	27.6	28.0	⋯	27.8
GOODS-S Wide	35.6	28.2	28.7	28.1	27.9	27.9	27.3	27.6	⋯	27.4
GOODS-N Deep	67.6	28.1	28.3	27.9	⋯	27.7	28.1	28.3	⋯	28.1
GOODS-N Wide	71.7	28.1	28.4	27.8	⋯	27.6	27.3	27.4	⋯	27.4
HUDF Main	5.1	29.5	30.0	29.7	⋯	29.1	29.9	29.6	29.6	29.7
HUDF PAR1	4.7	⋯	29.0	28.8	⋯	28.5	28.9	29.0	⋯	28.8
HUDF PAR2	4.8	⋯	29.0	28.7	⋯	28.3	28.9	29.2	⋯	28.9
MACS0416 PAR	4.4	28.8	28.9	⋯	29.2	⋯	29.2	29.0	29.0	29.0
Abell 2744 PAR	4.3	29.0	29.1	⋯	29.2	⋯	29.1	28.8	28.8	28.9
Zeropoints	⋯	25.68	26.51	25.67	25.95	24.87	26.27	26.23	26.45	25.95

Note. The magnitudes quoted are 5σ limits measured in 04-diameter apertures on non-PSF-matched images.

Download table as:  ASCIITypeset image

The HST imaging used here spans more than a factor of three in wavelength, thus the differences in PSF FWHM across that range are significant. For example, the PSF in the GOODS-S Deep field has a FWHM = 0193 in the H₁₆₀-band, but only 0119 in the B₄₃₅-band. A point source will thus have more of its flux contained within a 04 aperture in the B₄₃₅-band compared with the H₁₆₀-band. Because the selection of distant galaxies relies very heavily on accurate colors, and we are using apertures of fixed sizes (determined by the detection image, see Section 2.3) to measure photometry in all bands, the changing PSF needs to be addressed.

We corrected for this by matching the PSF of the HST imaging to the H₁₆₀-band image (which has the largest PSF FWHM) in each field by using the IDL deconv_tool Lucy–Richardson deconvolution routine in the same way as Finkelstein et al. (2010, 2012b). This routine requires the PSF for a given band, as well as a reference PSF (in this case, the H₁₆₀-band), and it generates a kernel. The PSFs were generated by stacking stars in each field in each band, where the stars were selected via identifying the stellar locus in a half-light radius versus magnitude plane. Each star was then visually inspected to ensure that there were no bright near-neighbors, and then the stars were stacked, subsampling by a factor of 10 to ensure an accurate centroiding of each star (i.e., to avoid smearing the PSF during the stacking). Using these PSFs, the deconvolution routine performed an iterative process, relying on the user to determine the number of iterations. We did this by making a guess as to the correct number of iterations, and then changing this number until the stars in the PSF-matched images in a given band had curves of growth that matched the H₁₆₀-band curves of growth to within 1% at a radius of 04. The images were then convolved with the final kernel to generate PSF-matched images.

Photometry was measured on the PSF-matched data set with a modified version of the Source Extractor software (v2.8.6, Bertin & Arnouts 1996). Our modified version adds a buffer between the source and the local background cell, and removes spurious sources associated with the distant wings of bright objects. Catalogs were generated independently in each of our 10 subfields, using Source Extractor in two-image mode, where the same detection image was used to measure photometry from all available HST filters. For most of our fields, we used a weighted sum of the F125W and F160W images as the detection image to increase our sensitivity to faint objects. In the HUDF main field and the MACS0416 and A2744 HFF parallel fields, we supplemented this catalog with catalogs using 10 additional detection images, derived by stacking all possible combinations of adjacent WFC3 filters. In these three fields, a combined catalog was made up of all unique sources in the catalogs, using a 02 matching radius. This allowed very blue sources to be included that may have been too faint in the H₁₆₀ image to be selected in the original F125W+F160W-selected catalog. This procedure was replicated in our completeness simulations (Section 4). To derive accurate flux uncertainties, Source Extractor relies on both an accurate rms map and a realistic estimate of the effective gain. The provided rms map has been shown to produce accurate uncertainties, and it has been corrected for pixel-to-pixel correlations that occur as a result of the drizzling process (see Guo et al. 2013), which are typically on the order of 10%–15% of the total rms. The effective gains were computed for each band separately, as the the instrument gain (1 for ACS, 2.5 for WFC3/IR) multiplied by the total exposure time for a given image. We have previously verified that the uncertainties measured in this manner on HST imaging are accurate (Finkelstein et al. 2012b). The zero-points to convert the observed fluxes into AB magnitudes are given in Table 1, and are appropriate for the dates when these data were taken.

Following our previous work (Finkelstein et al. 2010, 2012a, 2012b, 2013), colors were measured in small Kron apertures with the Source Extractor Kron aperture parameter PHOT_AUTOPARAMS set to values of 1.2 and 1.7. Finkelstein et al. (2012b) found that these apertures result in more reliable colors for faint galaxies when compared to isophotal or small circular apertures. An aperture correction to the total flux was derived in the H-band and was computed as the ratio between this small Kron aperture flux, and the default Source Extractor MAG_AUTO flux, which is computed with PHOT_AUTOPARAMS = 2.5, 3.5. These aperture corrections were then applied to the fluxes in all filters. To see if our aperture corrections accurately recovered the total flux, we examined our completeness simulations (discussed in Section 4) and found that after applying this aperture correction, recovered fluxes were typically 5% fainter in each band than their input fluxes (with the exception of the HUDF main field, where the measured correction was 2%). We thus increased the flux in each band by the appropriate factor to derive our best estimate of the total flux.

The Source Extractor catalogs from each band were combined into a master catalog for each field. At this step, the observed fluxes were corrected for Galactic extinction using the color excess $E(B-V)$ from Schlafly & Finkbeiner (2011) appropriate for a given field,¹⁶ and using the Cardelli et al. (1989) Milky Way reddening curve to derive the corrections based on each filter's central wavelength. We used a mask image to remove objects in regions of bad data, where the mask was generated using a threshold value from the weight map. This mask primarily trims off the noisier edges of the imaging, but it also excludes the "death star" region on the WFC3 array where the number of dithers was low (i.e., in the CANDELS Wide regions). The areas quoted in Table 1 are those of the good regions in these masks. Objects were also removed from the catalog if they had a negative aperture correction, which was the case for a very small number of sources, primarily restricted to areas near very bright sources where the flux in the larger aperture was unreliable. The remaining objects comprised our final catalog in each field.

We selected our candidate high-redshift galaxy sample via a photometric redshift fitting technique. This has the advantage in that it uses all of the available photometry simultaneously, rather than the multi-step Lyman break galaxy (LBG) method, which selects galaxies using two colors and then subsequently imposes a set of optical non-detection criteria (e.g., Bouwens et al. 2007, 2011b). Another advantage of photometric redshifts is that one obtains a redshift probability distribution function (PDF), which not only allows one to have a better estimate of the redshift uncertainty (${\sigma }_{z}$ typically ∼0.2–0.3 versus 0.5 for the LBG technique), but can also be used as a tool in the construction of the sample itself. A potential disadvantage of photometric redshift techniques is that the results are based on a set of assumed template spectra; if these templates do not encompass properties similar to the galaxies being studying, systematic offsets may occur (although we also note that similar templates are used to construct LBG color selection criteria). That being said, initial work comparing the differences between galaxy samples selected via LBG and photometric redshift techniques found that the resulting sample properties are fairly similar (e.g., McLure et al. 2013; Schenker et al. 2013).

Photometric redshifts for all sources in the catalogs for each field were measured using the EAZY software (Brammer et al. 2008). The input catalog used all available HST photometry, with the exception of the F814W imaging in the CANDELS fields, which was used solely for visual inspection (see Section 2.1). We used an updated set of templates provided with EAZY, based on the PÉGASE stellar population synthesis models (Fioc & Rocca-Volmerange 1997), which now include an increased contribution from emission lines, as recent evidence points to strong rest-frame optical emission lines being ubiquitous among star-forming galaxy populations at high redshift (e.g., Atek et al. 2011; Finkelstein et al. 2011, 2013; van der Wel et al. 2011; Stark et al. 2013; Smit et al. 2014). EAZY assumes the IGM prescription of Madau (1995), and does have the option to include magnitude priors when fitting photometric redshifts, using the luminosity functions as a prior for whether a galaxy at a given apparent magnitude resides at a given redshift. As we show later, there is still non-negligible uncertainty at the bright end of the luminosity function, therefore we did not include these magnitude priors during our photometric redshift fitting process.

We selected candidate galaxy samples in five redshift bins centered at $z\;\sim$ 4, 5, 6, 7, and 8 with ${\rm{\Delta }}\;z\;=$ 1, using criteria similar to our previous work (Finkelstein et al. 2012b, 2013). The cosmic time elapsed between our last two bins at $z\;\approx$ 7 and $z\;\approx$ 8 is ∼125 Myr. This time is much longer than the dynamical time of the systems we study, and thus leaves significant time for evolution. However, this will not always be the case as studies of galaxy evolution move toward higher redshift (e.g., ${\rm{\Delta }}\;{t}_{z=13\to 12}\;=$ 40 Myr), thus future studies with the James Webb Space Telescope (JWST) will need to pay careful attention to the choice of sample redshifts when studying galaxy evolution.

Rather than relying solely on the best-fit redshift value, we used the full redshift probability distribution curves P(z) calculated by EAZY (where $P(z)\;\propto$ exp($-{\chi }^{2}$), normalized to unity). Our selection criteria are:

• 1.  
  A ≥3.5 significance detection in both the J₁₂₅ and H₁₆₀ bands. A requirement of a significant detection in two bands removes nearly all spurious sources, since the chances of a noise peak occurring in two images at the same position are very small (Section 3.8.1). This requirement also limits our analysis to galaxies with $z\;\lt$ 8.5, because the Lyman break shifts into the J₁₂₅ band at z = 8.1.
• 2.  
  The integral of the redshift PDF under the primary redshift peak must comprise at least 70% of the total integral. This ensures that no more than 30% of the integrated redshift PDF can be in a secondary redshift solution.
• 3.  
  The integral under the redshift PDF in the redshift corresponding to a given sample (i.e., 6.5–7.5 for the z = 7 sample) must be at least 25%, which ensures that the redshift PDF is not too broad.
• 4.  
  The area under the curve in the redshift range of interest must be higher than the area in any other redshift range (i.e., for a galaxy in the z = 7 sample, the integral of $P(6.5\lt \;z\;\lt 7.5)$ must be higher than the integral in any other redshift bin). This criterion ensures that a given galaxy cannot be included in more than one redshift sample.
• 5.  
  At least 50% of the redshift PDF must be above ${z}_{\mathrm{sample}}-1$ (i.e., $\int P(z\gt 6)\gt 0.5$ for ${z}_{\mathrm{sample}}\;=$ 7), and the best-fit redshift must be above ${z}_{\mathrm{sample}}-2$.
• 6.  
  The minimum of ${\chi }^{2}$ from the fit must be less than or equal to 60. This criterion ensures that EAZY provides a reasonable fit, although in practice it does not reject many sources.
• 7.  
  Magnitude in the H₁₆₀ band must be ≥22. This effectively cleans many stars from our sample, but the limit is still more than two magnitudes brighter than our brightest $z\;\geqslant$ 6 galaxy candidate. At z = 4, we do have a few sources close to this limit, but only two sources are brighter than H = 22.4. This fact—coupled with the observation that the very few sources at $H\;\lt$ 22 that satisfy our z = 4 selection criteria are either obvious stars or diffraction spikes—implies that this criterion should not significantly affect our luminosity function results.

Of these criteria, items #1 and #2 are by far the most constraining, as most galaxies that meet these criteria with ${z}_{\mathrm{best}}\gt 3.5$ make it into our sample. Items #3 and #4 are responsible for putting a candidate galaxy in a given redshift sample. While some of these cuts are arbitrary, they will be corrected for when we apply the same criteria to our completeness simulations in Section 4. In Figure 2, we compare the photometric redshifts for 171 galaxies in our sample to available spectroscopic redshifts in the literature.¹⁷ The agreement is excellent, with ${\sigma }_{{\rm{\Delta }}z/(1+z)}\;=$ 0.031 (derived by taking an iterative 3σ-clipped standard deviation), although the number of confirmed redshifts at $z\;\gt$ 6.5 is small (only five galaxies). The number of outliers is also small, with only six out of 171 galaxies (3.5%) having a photometric redshift differing from the spectroscopic redshift by ${\rm{\Delta }}z\;\gt$ 1 at ≥3σ significance. All six galaxies have ${z}_{\mathrm{spec}}\;\gtrsim$ 4, thus no galaxies in our sample have a catastrophically lower spectroscopic redshift. In comparison, when defining outliers in the same way, we find that the published CANDELS team photometric redshift catalog has 13 outliers out of 174 total spectroscopic redshifts, for a somewhat higher outlier fraction of 7.5% (Dahlen et al. 2013). Although the fraction of galaxies with confirmed redshifts is relatively small, the available spectroscopy confirms that our selection methods yield an accurate high-redshift sample. In the remainder of this paper, we will therefore refer to our candidate galaxies solely as galaxies, with the caveat that spectroscopic follow up of a much larger sample, particularly at $z\;\gt$ 6, is warranted.

Zoom In Zoom Out Reset image size

The candidate selection process is automated, so for a truly robust galaxy sample we required a visual inspection of each of our ∼7500 candidate high-redshift galaxies. During the visual inspection, we examined the following features:

• 1.  
  Is the source a real galaxy? Objects were inspected to ensure that they were not an artifact, such as part of a diffraction spike (which frequently appear in different places in the ACS and WFC3 imaging due to different roll angles during the respective observations), oversplit regions of bright galaxies, or noise near the edge of the images.
• 2.  
  Is the aperture drawn correctly? While the small Kron apertures yield the most reliable colors, they are also susceptible to "stretching" (i.e., becoming highly elongated) in regions of high noise or near very bright objects. For each source, we compared the ratio of the flux between the Kron aperture and a 04-diameter circular aperture to that same quantity for objects of a similar magnitude from the full photometry catalog. If an object had a value ≳30% higher than similarly bright sources in the full photometry catalog and the aperture looked affected by noise/bright sources, we adjusted the photometry of the object in question accordingly, using the 04-to-total correction of similarly bright objects in the catalog. In practice, these issues affected <10% of the galaxies in our high-redshift sample.
• 3.  
  Is there significant optical flux that did not get measured correctly? Primarily due to the issues with the inaccurate apertures discussed previously, a very small number of sources appeared to have optical flux when visually inspected that was not measured to be significant in our catalog (i.e., in the case of a too-large aperture, the flux is concentrated in a small number of pixels, whereas the flux error comes from the full aperture, so the signal-to-noise is low). In these cases, objects were removed from our sample. This step is somewhat qualitative, as there are cases of objects where the aperture appears correct, yet there is still a ∼1-2σ detection in a single optical band. In the majority of these cases, we left these objects in the sample because we are confident in our photometric redshift analysis. During this step, we also examined I₈₁₄ photometry for each source in the CANDELS fields; this primarily benefits the selection of z = 8 galaxies, which should not be visible at this wavelength. Three z = 8 candidates with observable I₈₁₄ flux were removed from our sample.

The most crucial step in our visual inspection is the classification and removal of stellar sources, because stellar contamination would dominate the bright end of the luminosity function if these contaminants were not considered. In particular, M-type stars and L- and T-brown type dwarf stars can have similar colors (including optical non-detections) as our high-redshift galaxies of interest, especially at $z\;\geqslant$ 6. While some studies use dwarf star colors during their selection (e.g., Bowler et al. 2012, 2014; McLure et al. 2013; Bouwens et al. 2015), many primarily use the Source Extractor "stellarity" parameter to diagnose whether a compact object is a star or a galaxy (e.g., Bouwens et al. 2011b, 2015). However, the stellarity parameter loses its ability to discern between a point source and a resolved source for faint objects. To test this further, we examined the stellarity of sources in the CANDELS GOODS catalogs. At very bright magnitudes (${J}_{125}\;\lt$ 24), there is a clear separation between stars and galaxies, with objects either having a stellarity near unity (i.e., stars) or a stellarity near zero (i.e., galaxies). However, this separation becomes less clear at ${J}_{125}\;\gt$ 25, where the stellar and galaxy sequences begin to blend together. Therefore, stellarity can be an unreliable star-galaxy separator at ${J}_{125}\;\gt$ 25, which is similar to the brightness of our brightest $z\;\geqslant$ 7 galaxies.

While the GOODS fields cover relatively small regions on the sky, the potential number of brown dwarf contaminants, even at ${J}_{125}\;\gt$ 25, is significant. The Galactic structure model of R. Ryan & N. Reid (2015, in preparation) predicts the surface density of brown dwarfs in our covered fields. In the GOODS-S region, using the area covered by the CANDELS, ERS, and HUDF09 observations, we would expect ∼6 stars of spectral type M6–T9 with J₁₂₅-band magnitudes between 25 and 27. The surface density of M6-T9 stars in GOODS-N is similar, with an expected number of stars in the field of ∼5. Thus, the expected number of 25 $\lt \;{J}_{125}\;\lt$ 27 stars of spectral type M6-T9 in our whole surveyed region is ∼11. While this number is small, the number of brown dwarfs is expected to fall off toward fainter magnitudes, thus the majority of these likely have J₁₂₅ close to 25. This magnitude is similar to those of the brightest galaxies in our sample, which dominate the shape of the bright end of the galaxy luminosity function. As stellarity is an unreliable method of identifying these sources, we must find an alternative method.

Although brown dwarfs can have similar colors to $z\;\gt$ 6 galaxies, and can be included in the initial sample, they fall on well-defined color sequences, and can thus be distinguished from true galaxies. Figure 3 shows the two color–color plots,¹⁸ which we used in tandem with the size information—examining not only stellarity, but also the FWHM and half-light radius as measured by Source Extractor—to identify stars lurking our sample (similar plots were used at z = 4 and 5). If a galaxy appeared un-resolved (defined as having a stellarity >0.8, or a half-light radius and/or FWHM similar to that of stars in the field), then we examined that object in the color–color plots as shown in Figure 3. If the object also had colors similar to a dwarf star, we removed it from our sample. Over all of our fields, a total of 23 objects were flagged as stars (many with $J\;\lt$ 25) in our $z\;\geqslant$ 6 samples: 18 from our initial $z\;\sim$ 6 galaxy sample and five from our initial $z\;\sim$ 7 galaxy sample. These objects were removed from our sample. One of the stars removed from our $z\;\sim$ 7 sample was previously flagged as a probable T-dwarf by Castellano et al. (2010). We examined a subset of eight of these stars, which were detected in the FourStar Galaxy Evolution (zFourGE) medium band imaging survey of a portion of GOODS-S, and found that all eight have $z-J1$ and $J1-J3$ colors consistent with brown dwarf stars (where J1 and J3 refer to two of the three medium bands that comprise the J band; Tilvi et al. 2013). Of these, six stars have ${J}_{125}\;\gt$ 25, meaning that our high-redshift galaxy selection criteria also originally selected about half of the expected number of faint brown dwarfs in this field. Four of these six stars have Source Extractor stellarity measurements <0.8, thus a stellarity-only rejection method would have failed to remove them. We conclude that our visual inspection step efficiently removed stellar contaminants from our sample, but we emphasize that the color examination portion was crucial to exclude the faintest stars from our sample.

Zoom In Zoom Out Reset image size

We screened for the presence of bright AGN in our sample by searching for counterparts in Chandra X-ray Observatory point source catalogs. In the GOODS-S field, we used the 4 Msec Chandra Deep Field—South (CDF-S) catalog of Xue et al. (2011); in GOODS-N, we used the 2Msec Chandra Deep Field—North catalog of Alexander et al. (2003). These catalogs have average positional accuracies of 042 and 03, respectively. To be conservative, we searched for matches in each catalog out to a radius of 1${}^{\prime\prime }$. We then visually inspected each of the 34 galaxies in our sample with a match. Seven objects, all with Chandra catalog separations >06, had nearby HST counterparts with positions that were consistent with the Chandra catalog. These seven sources remained in our sample, because these interlopers were likely providing the X-ray emission, and none had spectroscopic redshifts in the CDF-S catalog. The remaining 27 sources, all with Chandra catalog separations $\leqslant 0\buildrel{\prime\prime}\over{.} 6$ from our candidate galaxies had HST counterparts with positions consistent with the X-ray emission coming from the galaxies in our sample. Secure spectroscopic redshifts were available for four of these 27 galaxies in the CDF-S catalog: z = 3.06, 3.66, 3.70, and 4.76. We removed these 27 galaxies (25 from our z = 4 sample, and two from our z = 5 sample) from our galaxy sample. This removal is conservative, because, although the X-ray detections imply the presence of an AGN, they do not prove that the AGN dominates the UV luminosity.

As we will discuss later, one of the main results of this work is an apparent constant value of ${M}_{\mathrm{UV}}^{*}$ at $z\;\gt$ 5, which is brighter than many previous works. It is thus imperative that we have high confidence that our bright galaxies are all actually at high redshift, and not lower-redshift contaminants. To provide a further check on our bright sources, we re-examined the photometric redshifts of our bright galaxies with the addition of Spitzer Space Telescope Infrared Array Camera (IRAC; Fazio et al. 2004) imaging over our fields. This imaging probes the rest-frame optical at these wavelengths, and thus provides significant constraining power because the most likely contaminants are red, lower-redshift galaxies, which would have very different fluxes in the mid-infrared than true high-redshift galaxies. We examined sources with ${M}_{1500}\lt -21$, which is approximately the value of ${M}_{\mathrm{UV}}^{*}$ at these redshifts, and provides samples of 164, 85, 29, 18, and 3 bright galaxies at z = 4, 5, 6, 7, and 8, respectively.

During the cryogenic mission, the GOODS fields were observed by the GOODS team (M. Dickinson et al. 2015, in preparation) at 3.6, 4.5, 5.8, and 8.0 μm. Later, during Cycle 6 of the warm mission, broader regions encompassing the GOODS footprints were covered by the Spitzer Extended Deep Survey (SEDS Ashby et al. 2013) to 3σ depths of 26 AB mag at both 3.6 and 4.5 μm. A somewhat narrower subset of both fields was subsequently covered by Spitzer-CANDELS (S-CANDELS Ashby et al. 2015) to even fainter levels; reaching ∼0.5 mag deeper than SEDS in both of the warm IRAC bandpasses. The HUDF09 fields were observed by Spitzer program 70145 (the IRAC Ultra-Deep Field; Labbé et al. 2013), reaching 120, 50, and 100 hr in the HUDF Main, PAR1, and PAR2 fields, respectively. Finally, program 70204 (PI Fazio) observed a region in the ERS field to 100 hr depth. The present work is based on mosaics constructed by coadding all the above data, following the procedures described by Ashby et al. (2013). The combined data have a depth of $\gtrsim 50$ hr over both CANDELS GOODS fields and >100 hr over the HUDF main field.

Because the IRAC PSF is much broader than that of HST, our galaxies may be blended with other nearby sources. We measure Spitzer/IRAC 3.6 and 4.5 μm photometry by performing PSF-matched photometry on the combined IRAC data, which reach at least 26.5 mag (3σ) at 3.6 and 4.5 μm (Ashby et al. 2015). We utilized the TPHOT software (Merlin et al. 2015), an updated version of TFIT (Laidler et al. 2007), to model low-resolution images (IRAC images) by convolving HST imaging with empirically derived IRAC PSFs and simultaneously fitting all IRAC sources. Specifically, we used the light profiles and isophotes in the detection ($J+H$) image obtained by Source Extractor, and convolved them with a transfer kernel to generate model images for the low-resolution data. These models were then fit to the real low-resolution images, dilating the segmentation maps of the model images to account for missing flux on the edges of galaxies (Galametz et al. 2013). The fluxes of sources are determined by the model that best represents the real data. As the PSF FWHM of the high-resolution image (H-band) is negligible (∼019) when compared to those of the low-resolution IRAC images (∼17), we use the IRAC PSFs as transfer kernels. We derive empirical IRAC PSFs by stacking isolated and moderately bright stars in each field. Because our own WFC3 catalog was used as the input for TPHOT, all of our galaxies have IRAC measurements in the TPHOT catalogs. We visually inspected the positions of each of our high-redshift galaxy candidates in the IRAC images to ensure no significant contamination from the residuals of nearby bright galaxies. If an object was on or near a strong residual, we ignored the IRAC photometry in the subsequent analysis. This was the case for 18 of the 164 galaxies at z = 4, 23 of the 85 galaxies at z = 5, 3 of the 29 galaxies at z = 6, 6 of the 18 galaxies at z = 7, and 1 of the 3 galaxies at z = 8. With the contaminated fluxes removed, we found that all remaining ${M}_{1500}\lt -21$ galaxies at z = 4–8 had 3.6 μm detections of at least 3σ significance, with a magnitude range at $z\;\geqslant$ 6 of 22.7 ≤ ${m}_{3.6}\;\leqslant$ 25.8. The full description of our TPHOT IRAC photometry catalog will be presented by Song et al. (2015).

We reran EAZY for this subsample of bright galaxies, including the Spitzer/IRAC fluxes, as well as photometry from the ACS F814W filter, which was not included in the original photometric redshift calculation (see Section 2.1). We examined these updated photometric redshift results, searching for galaxies in our z = 4 and 5 samples with ${z}_{\mathrm{new}}\;\lt$ 2.5, and in our z = 6, 7, and 8 samples with ${z}_{\mathrm{new}}\;\lt$ 4. We found 14 out of 164 galaxies in our z = 4 sample, and 14 out of 85 galaxies in our z = 5 sample with ${z}_{\mathrm{new}}\;\lt$ 2.5. We found 1 galaxy out of 29 at z = 6 that appeared to be better fit with a low-redshift solution of ${z}_{\mathrm{new}}\;=$ 0.9, but no galaxies in our z = 7 or 8 samples had preferred low-redshift solutions with the inclusion of IRAC photometry.

Examining these results, out of the 28 z = 4 or 5 galaxies with preferred low-redshift solutions, 23 had photometry consistent with a true low-redshift galaxy. Four galaxies, however, had photometry that appeared to be consistent with a high-redshift galaxy with a strong emission line (Hα or [O iii]) in one IRAC band. Systems with lines such as these (i.e., EW_{[O iii]} > 500 Å) are rare locally, but appear to be more common at high redshift (e.g., van der Wel et al. 2011; Finkelstein et al. 2013; Smit et al. 2014). Although typical emission line strengths are now included in the EAZY templates, these do not account for extreme emission lines, so it is not surprising that EAZY does not return a high-redshift solution. We elect to keep these four galaxies in our sample, noting that the lack of strong rest-frame optical lines in the EAZY templates does not affect our initial sample selection, which does not make use of the IRAC photometry. A fifth galaxy (z5_GNW_13415) has a high-quality published spectroscopic redshift of z = 5.45, thus we also keep it in our sample. Perhaps unsurprisingly, the low-redshift fit for this galaxy had a very poor quality-of-fit, with ${\chi }^{2}\;=$ 127, implying that the EAZY templates are a poor match for this galaxy. The remaining 23 galaxies at z = 4 and 5 were removed from our sample. The sole $z\;\geqslant$ 6 galaxy with a preferred low-redshift solution with the inclusion of IRAC photometry, z6_GSW_3089, is shown in Figure 4. The red HST colors imply a much brighter IRAC flux than is seen. A solution at z = 0.93 yields a better fit, because the peak of stellar emission at that redshift better matches the observed IRAC fluxes. We thus removed this galaxy from our sample. This galaxy also has a spectroscopic redshift of z = 5.59 from the observations of Vanzella et al. (2009); however, it has a spectroscopic quality flag of "C," which indicates that the spectroscopic redshift is unreliable. This combined with the ∼4σ detection in the V₆₀₆ band leaves us confident that a low-redshift solution is more likely.

Zoom In Zoom Out Reset image size

With the removal of these likely contaminants, we retain a total sample of 150, 77, 28, 18, and 3 ${M}_{1500}\lt -21$ galaxies at z = 4, 5, 6, 7, and 8, respectively. The fraction of contaminants at $z\;\geqslant$ 6 (1 out of 50 z = 6, 7, and 8 galaxies, or 2%) is consistent with (albeit somewhat less than) the expected low value calculated in Section 3.8.

Our final galaxy sample is summarized in Table 2, and a catalog of all galaxies in our sample is provided in Table 3. In Figure 1 we show the absolute UV magnitude distribution of our final samples, highlighting that we cover a dynamic range of five magnitudes. In particular, the CANDELS data are crucial, because galaxies from these data dominate the total number of galaxies in our sample, and approximately double the luminosity dynamic range that we can probe. This is highlighted in the left panel of Figure 2, which shows that galaxies discovered in the CANDELS GOODS fields dominate the total number at all redshifts in our sample.

Table 2.  Summary of Final High-redshift Galaxy Samples

Redshift	Nall	N${}_{M\lt -21}$	Veff (M${}_{1500}=-$22)	Veff (M${}_{1500}=-$19)
			(105 Mpc3)	(105 Mpc3)
4 (3.5–4.5)	4156	150	12.2	4.11
5 (4.5–5.5)	2204	77	8.98	3.36
6 (5.5–6.5)	706	28	7.93	2.50
7 (6.5–7.5)	300	18	6.99	0.30
8 (7.5–8.5)	80	3	5.88	0.16

Note. The total number of sources in our final galaxy sample, after all contaminants were removed. The final two columns give the total effective volume at each redshift for two different values of the UV absolute magnitude.

Download table as:  ASCIITypeset image

To derive the rest-frame absolute magnitude at 1500 Å (M₁₅₀₀), as well as the UV spectral slope β (f ${}_{\lambda }\propto {\lambda }^{\beta }$ Calzetti et al. 1994), we fit spectral energy distributions (SEDs) from synthetic stellar population models to the observed HST photometry of our high-redshift candidate galaxies. The technique used here is similar to our previous works (Finkelstein et al. 2010, 2012a, 2012b, 2013). We used the updated (2007) stellar population synthesis models of Bruzual & Charlot (2003) to generate a grid of spectra, varying the stellar population metallicity, age, and star-formation history (SFH).¹⁹ Metallicities spanned 0.02–1 times ${z}_{\odot }$, and ages spanned 1 Myr to the age of the universe at a given redshift. We allowed several different types of SFHs, including a single burst and continuous, as well as both exponentially decaying and rising (so-called "tau" and "inverted-tau" models). To these spectra, we added dust attenuation using the starburst attenuation curve of Calzetti et al. (2000), with a range of 0 ≤ $E(B-V)$ ≤ 0.8 (0 ≤ A_V ≤ 3.2 mag). We also included nebular emission lines using the prescription of Salmon et al. (2015)—which uses the line ratios from Inoue (2011), based on the number of ionizing photons from a given model—and assuming the ionizing photon escape fraction is ≈ zero. We then redshifted these models to 0 $\lt \;z\;\lt$ 11 and added IGM attenuation (Madau 1995). These model spectra were integrated through our HST filter bandpasses to derive synthetic photometry for comparison with our observations. For each model, we computed the value of M₁₅₀₀ by fitting a 100 Å-wide synthetic top-hat filter to the spectrum centered at rest-frame 1500 Å. Likewise, for each model we measured the value of β by fitting a power law to each model spectrum using the wavelength windows specified by Calzetti et al. (1994), similar to Finkelstein et al. (2012b).

The best-fit model was found via ${\chi }^{2}$ minimization, including an extra systematic error term of 5% of the object flux for each band to account for such items as residual uncertainties in the zeropoint correction and PSF-matching process. The stellar mass was computed as the normalization between the best-fit model (which was normalized to 1 M${}_{\odot }$) and the observed fluxes, weighted by the signal-to-noise in each band. These best-fit values of M₁₅₀₀ are used in our luminosity function analysis, while β is used to correct for incompleteness in a color-dependent fashion. The uncertainties in the best-fit parameters were derived via Monte Carlo simulations, perturbing the observed flux of each object by a number drawn from a Gaussian distribution with a standard deviation equal to the flux uncertainty in a given filter. For each galaxy, 10² Monte Carlo simulations were run, providing a distribution of 10² values for each physical parameter. The 68% confidence range for each parameter was calculated as the range of the central 68% of results from these simulations. The best-fit redshift in these simulations was allowed to vary following the redshift PDF, thus folding in the uncertainty in redshift into the uncertainty in the physical parameters (most notably, the stellar mass and M₁₅₀₀; Finkelstein et al. 2012a). During this process, we only allowed the redshift to vary within ${\rm{\Delta }}z=\pm 1$ of the best-fit photometric redshift. This excludes any low-redshift solution from biasing the uncertainties on a given parameter. The amount of the integrated P(z) at $z\;\gt$ 3 excluded via this step was typically ≤10% (at z = 6).

A key issue in any study of high-redshift galaxies is the risk of sample contamination, either by spurious sources or by lower-redshift interlopers. The gold standard for eliminating contamination is to obtain spectroscopic redshifts. This is clearly unfeasible for all galaxies in our sample (until the next generation of space and ground-based telescopes), but there is significant archival spectroscopic data. As discussed in Section 3.2, we find excellent agreement between available spectroscopic redshifts and our photometric redshifts, with ${\sigma }_{z/(1+z)}\;=$ 0.031. In particular, the four bright²⁰ galaxies ($24.9\lt {J}_{125}\lt 25.7$) with confirmed ${z}_{\mathrm{spec}}\gt 6.5$ have an excellent agreement with spectroscopic redshifts: z7_GNW_24443 with ${z}_{\mathrm{phot}}\;=$ 6.66 and ${z}_{\mathrm{spec}}\;=$ 6.573 (Rhoads et al. 2013); z7_GSD_21172 with ${z}_{\mathrm{phot}}\;=$ 6.73 and ${z}_{\mathrm{spec}}\;=$ 6.70 (Hathi et al. 2008); z7_GNW_4703 with ${z}_{\mathrm{phot}}\;=$ 7.19 and ${z}_{\mathrm{spec}}\;=$ 7.213 (Ono et al. 2012); and z7_GND_42912²¹ with ${z}_{\mathrm{phot}}\;=$ 7.45 and ${z}_{\mathrm{spec}}\;=$ 7.51 (Finkelstein et al. 2013). The brightest source in our z = 6 sample, z6_GSW_12831 with ${M}_{1500}=-22.1$ and ${z}_{\mathrm{phot}}\;=$ 5.77, is confirmed with ${z}_{\mathrm{spec}}\;=$ 5.79 (Bunker et al. 2003). This galaxy has a 3σ detection in the V₆₀₆-band, which could have resulted in its exclusion from a typical LBG color–color selection sample, although the observed flux can be explained by non-ionizing UV photons transmitted through the Lyα forest.

In general, spectroscopic follow up of sources selected on the basis of their Lyman breaks (either color–color selection, or photometric redshift selection) finds a very small contamination rate by low-redshift sources (e.g., Pentericci et al. 2011). However, given the apparent difficulty in detecting Lyα emission at $z\;\gt$ 6.5 (e.g., Pentericci et al. 2011; Ono et al. 2012; Finkelstein et al. 2013), the true effect of contamination at these higher redshifts is not empirically known. In this subsection, we will attempt to estimate our contamination fraction by other means.

3.8.1. Properties of the Image Noise

3.8.2. Estimates from Stacked Redshift PDFs

Contamination by low-redshift interlopers is a more complicated issue. While extreme emission line galaxies at lower redshift could theoretically be an issue (Atek et al. 2011), our requirement of detections in two bands (as well as the frequent detections in more than two bands for all but the highest redshift objects in the z = 8 sample) makes  a significant contamination by these sources unlikely. The most likely possible contaminants are faint red galaxies at $z\;\leqslant$ 2 (e.g., Dickinson et al. 2000). These galaxies can be too faint to be detected in our optical imaging, but their red SEDs yield detections in the WFC3/IR bands. Although faint sources that are very red will have a disfavored high-redshift solution with our current photometric selection, we have information on their likelihoods encoded in our redshift PDFs. Figure 5 shows the redshift PDFs of galaxies in each of our three highest redshift samples, stacked in magnitude bins of ΔM = 1 mag. At all redshifts and all magnitudes, ≳85% of the redshift PDF is at $z\;\gt$ 4, implying that there is not significant contamination by lower-redshift galaxies.

Zoom In Zoom Out Reset image size

The position of the secondary redshift peak is consistent with the detection of a 4000 Å break, rather than the Lyman break (at z = 6, 7, and 8, this gives ${z}_{\mathrm{secondary}}\;=$ 1.1, 1.4, and 1.7, respectively). At z = 8, the possible contamination from galaxies at $z\;\lt$ 4 is <10.5%, which is primarily due to the fact that at z = 8, a galaxy will have to be undetected in most of the filters we consider here (there is an additional ∼8% chance the true redshift is in the range $4\lt \;z\;\lt 6$). The worst case is for faint galaxies at z = 6, because z = 6 galaxies are typically detected in all but two filters, although even here, the indicated contamination by $z\;\lt$ 4 galaxies is ≲15%.

3.8.3. Stacking Imaging

The limits from the previous subsection are likely upper limits on the contamination fraction. When fitting photometric redshifts to rule out all low-redshift solutions, the Lyman break needs to be detected at high significance, which is the case for only the brightest galaxies (e.g., at z = 6, the brightest bin has a contamination of <2%). Additionally, these results are dependent on the templates used, which by definition do not account for unknown galaxy populations. We therefore consider two empirical tests of contamination. The first is to stack all galaxies in a sample to search for detections below the Lyman break; the results from this test for z = 6, 7, and 8 are shown in Figure 6. As expected for galaxies at the expected redshifts, there is no visible signal in the B₄₃₅-band at z = 6, i₇₇₅-band and blueward at z = 7, and I₈₁₄-band and blueward at z = 8. This confirms that the majority of the flux in our sample does not arise from lower-redshift sources.

Zoom In Zoom Out Reset image size

3.8.4. Estimates from Dimmed Real Sources

As a final test, we estimated the contamination by artificially dimming real lower redshift sources in our catalog, to see if the increased photometric scatter allows them to be selected as high-redshift candidates. This empirical test is useful because it does not rely on known spectral templates to derive the contamination, although it does assume that the fainter objects that could potentially contaminate our sample have similar SEDs to the bright objects that we dim. We performed this exercise twice, once using the combined catalog of the GOODS-S and GOODS-N Deep fields, and once in the HUDF main field, to probe fainter magnitudes. In the GOODS Deep fields, we selected all real sources with 21 $\lt \;{H}_{160}\;\lt$ 24 and ${z}_{\mathrm{phot}}\lt 3$, and reduced their observed fluxes by a factor of 20. The same was done for sources drawn from the HUDF Main field, here extending the magnitude range to be 21 $\lt \;{H}_{160}\;\lt$ 26. The limits on these magnitudes were chosen to exclude any real high-redshift sources. We replaced the true flux uncertainties of these objects with flux uncertainties from a randomly drawn real source from the full catalog from a given field with a similar magnitude as the dimmed source. We then added scatter to the dimmed fluxes, perturbing them by a random amount drawn from a Gaussian distribution with a standard deviation equal to the flux uncertainties of the object. We included two realizations of the HUDF field to increase the number of dimmed objects.

The total number of sources in our artificially dimmed catalog was 4066 in the Deep fields and 1254 in the HUDF field. We measured the photometric redshifts of these sources with EAZY in an identical manner as on our real catalogs, and then applied our sample selection to this dimmed catalog. In the Deep fields, we found that a total of 149, 134, 54, 23, and 8 dimmed objects satisfied our z = 4, 5, 6, 7, and 8 selection criteria. Investigating the original (not dimmed or perturbed) colors of these sources, we found that they are unsurprisingly red, with the bulk of sources having $V-H\;\gt$ 2 mag. Thus this parent population of red sources is responsible for the majority of the possible contamination. The contamination fraction in our high-redshift sample is then defined as

$$\begin{eqnarray}&&{\mathcal{F}}=\displaystyle \frac{\displaystyle \frac{{N}_{\mathrm{dimmed},\mathrm{select}}}{{N}_{\mathrm{dimmed},\mathrm{red}}}\ast {N}_{\mathrm{total},\mathrm{red}}}{{N}_{z}}\end{eqnarray} \tag{ 1 }$$

where ${N}_{\mathrm{dimmed},\mathrm{select}}$ was the number of dimmed sources satisfying our high-redshift sample selection; ${N}_{\mathrm{dimmed},\mathrm{red}}$ was the total number of sources in the dimmed catalog with original colors of $V-H\;\lt$ 2; ${N}_{\mathrm{total},\mathrm{red}}$ was the number of sources in the full object catalog with 25 $\lt \;H\;\lt$ 27, ${z}_{\mathrm{phot}}\;\lt$ 3, and $V-H\;\lt$ 2; and N_z was the number of true galaxy candidates in a given redshift bin. For example, at z = 6, where we found 54 dimmed galaxies that satisfied our selection criteria (${N}_{\mathrm{dimmed},\mathrm{select}}=54$), ${N}_{\mathrm{dimmed},\mathrm{red}}\;=$ 1023, ${N}_{\mathrm{total},\mathrm{red}}\;=$ 695, and ${N}_{z}\;=$ 322, giving an estimated contamination fraction ${\mathcal{F}}\;=$ 11.4%. Thus for sources with 25 $\lt \;H\;\lt$ 27, we found an estimated contamination fraction of ${\mathcal{F}}\;=$ 4.5%, 8.1%, 11.4%, 11.1%, and 16.0% at z = 4, 5, 6, 7, and 8, respectively. We performed the same exercise in the HUDF, for fainter sources with 26 $\lt \;H\;\lt$ 29, finding that 30, 21, 8, 8, and 0 sources satisfied our z = 4, 5, 6, 7, and 8 selection criteria, respectively, giving a contamination fraction of 9.1%, 11.6%, 6.2%, 14.7%, and <4.9% at z = 4, 5, 6, 7, and 8, respectively.

Broadly speaking, we estimate relatively small contamination fractions of ∼5%–15%, which is in line with the estimates from the stacked P(z) curves. As the bulk of contaminants appear to be red galaxies, it is interesting to compare to the space density of these potentially contaminating sources. This was recently estimated by Casey et al. (2014b), who found that dusty star-forming galaxies at $z\;\lt$ 5 will contaminate $z\;\gt$ 5 galaxy samples at a rate of <1%. This is much less than our contamination estimates, thus we may have overestimated the contamination rate; although it may not be inconsistent once photometric scatter is applied to faint, red galaxies, making it easier for them to scatter into our sample. In any case, the expected contamination rate is quite small, therefore we do not reduce our observed number densities for the expected minimal contamination.

Now that we have a final galaxy sample with measured values of M₁₅₀₀, as well as the effective volumes from our completeness simulations in the previous section, we can proceed to measure the rest-frame UV luminosity function at z = 4, 5, 6, 7, and 8. We calculate the luminosity function in two ways: a parametric version assuming that the luminosity function takes the form of a Schechter (1976) function, and a non-parametric stepwise maximum likelihood (SWML) calculation.

The fitting of a Schechter function is well motivated, as it successfully matches the observed rest-UV luminosity functions at lower redshifts (e.g., Bouwens et al. 2006; Reddy & Steidel 2009). This function is characterized by a power law at the faint end with slope α, and an exponential cutoff at the bright end, transitioning between the two regimes at the characteristic magnitude $M*$. The parameter $\varphi *$ sets the normalization of this function. The number density at a given magnitude is then given by

$$\begin{eqnarray}&&\varphi (M)=0.4\;\mathrm{ln}\;(10)\;\varphi *\;{10}^{-0.4(M-M*)(\alpha +1)}\;{e}^{-{10}^{-0.4(M-M*)}}.\end{eqnarray} \tag{ 4 }$$

For the measurement of the luminosity function assuming a Schechter functional form, we calculated the likelihood that the number of observed galaxies in a given magnitude bin is equal to that for an assumed value of the Schechter parameters M* and α. Rather than performing a grid-based search, we performed a Markov Chain Monte Carlo (MCMC) search algorithm to better span the parameter space, as well as to better characterize the uncertainties on the Schechter parameters. We performed this calculation in bins of absolute magnitude with ΔM = 0.5 mag, ranging from −24 ≤ ${M}_{1500}\leqslant -17$. At the bright end we are in the limit of small numbers, therefore we model the probability distribution as a Poissonian distribution (e.g., Cash 1979; Ryan et al. 2011), with:

$$\begin{eqnarray}&&{{\mathcal{C}}}^{2}(\varphi )=-2\;\mathrm{ln}\;{\mathcal{L}}(\varphi )\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}{C}^{2}(\varphi ) & = & -2\displaystyle \sum _{i}\displaystyle \sum _{j}\;{N}_{j,\mathrm{obs}}\;\mathrm{ln}({N}_{j,\mathrm{model}})\\ & & -{N}_{j,\mathrm{model}}-\mathrm{ln}({N}_{j,\mathrm{obs}}!)\end{eqnarray} \tag{ 6 }$$

where ${\mathcal{L}}(\varphi )$ is the likelihood that the expected number of galaxies (N_model) matches that observed (N_obs) for a given value of M* and α, and C² is the goodness-of-fit statistic. The subscripts i and j represent the subfields and magnitude bins, respectively. The final goodness-of-fit is the sum over all fields and magnitudes in a given redshift bin. We use the effective volume results for a given redshift, magnitude bin, and field to convert from the model number density to the expected number, calculating $\varphi *$ as the normalization, such that the total expected number of galaxies over all magnitude bins matches the total number of observed galaxies.

For each magnitude bin, we performed 10 independent MCMC chains utilizing a Metropolis–Hastings algorithm, each of 10⁵ steps, building a distribution of M*, α, and $\varphi *$ values for each field. During each step of the chain, the likelihood of a given model was computed for each of our observed fields, and then added together to compute the likelihood for the sample as a whole (we also recorded the individual field values, see Section 6.5). Prior to each recorded chain, we performed a burn-in run with a number of steps equal to 10% of the number of steps in each chain. The starting point for the burn is a brute-force ${\chi }^{2}$ fit of a grid of α and M* values to our data. At the end of the burn, the final values of the parameters from the last step were then the starting points for each chain. The rest of the burn-in results were not recorded. During each step, new values of M* and α were chosen from a random Gaussian distribution, with the Gaussian width tuned to generate an approximate acceptance rate of 23%. During each step $\varphi *$ was calculated as the normalization. If the difference between the likelihood of the model for the current step exceeds that from the previous step by more than a randomly drawn value (i.e., ≡2 ln (n); where n is a uniform random number between zero and unity), then the current values of the Schechter function parameters were recorded. If not, the chain reverted to the value from the previous step.

By running 10 independent chains, we mitigate against being trapped by local minima in the parameter space. Our final result links these 10 chains together, giving a distribution of 10⁶ values of the Schechter function parameters at each redshift. The results were visually inspected to confirm that the chains reached convergence. For each Schechter function parameter, the best-fit values were taken to be the median of the distribution, with the uncertainty being the central 68% of the distribution. These results are given in Table 4. For the z = 8 Schechter function fit, we imposed a top-hat prior, forcing ${M}_{\mathrm{UV}}^{*}$ to be fainter than −23. Without this prior, the fit preferred a much brighter value of ${M}_{\mathrm{UV}}^{*}$, such that the observed data points all lay on the faint-end slope (i.e., a single power law). We discuss the implications of this in Section 6.6.

Table 4.  Schechter Function Fits to the Luminosity Function

Redshift	${M}_{\mathrm{UV}}^{*}$	α	$\varphi *$
			(Mpc−3)
4	−20.73${}_{-0.09}^{+0.09}$	−1.56${}_{-0.05}^{+0.06}$	(14.1${}_{-1.85}^{+2.05}$) × 10−4
5	−20.81${}_{-0.12}^{+0.13}$	−1.67${}_{-0.06}^{+0.05}$	(8.95${}_{-1.31}^{+1.92}$) × 10−4
6	−21.13${}_{-0.31}^{+0.25}$	−2.02${}_{-0.10}^{+0.10}$	(1.86${}_{-0.80}^{+0.94}$) × 10−4
7	−21.03${}_{-0.50}^{+0.37}$	−2.03${}_{-0.20}^{+0.21}$	(1.57${}_{-0.95}^{+1.49}$) × 10−4
8	−20.89${}_{-1.08}^{+0.74}$	−2.36${}_{-0.40}^{+0.54}$	(0.72${}_{-0.65}^{+2.52}$) × 10−4

Note. The final values for each parameter are the median of the parameter distribution from the MCMC analysis. The quoted errors represent the 68% confidence range on each parameter.

Download table as:  ASCIITypeset image

Although we computed the volumes down to very faint magnitudes, we would be highly incomplete if we included these faint galaxies and calculated the luminosity function down to ${M}_{1500}=-17$ or fainter (Figure 9). In practice, it is our deepest field (i.e., the HUDF) that determines how faint we can constrain the luminosity function. The HUDF drops below 50% completeness at magnitudes fainter than ${M}_{1500}\sim -17.5$ at z = 4, 5, and 6; −18 at z = 7; and −18.5 at z = 8. Thus, in our calculation of the luminosity function, we only include a given field's contribution at a given magnitude if it is above the 50% completeness limit for that magnitude and redshift. Extending the analysis fainter will give results that are dominated by the incompleteness correction.

As shown in Figure 9, while the volume per unit area for the different fields is very tight at z = 4, there is a progressively larger scatter apparent when moving toward higher redshift, representing a systematic uncertainty in the effective volume calculation. One likely culprit is the fact that the volumes depend on the distribution of the sizes and colors of objects in a given field. For fields with few sources (i.e., the smallest fields at the highest redshifts), there may only be a single object in a given magnitude bin. To mitigate significant variances in the effective volume at the bright end, where the numbers of sources are small, we set the effective volume in a given redshift bin and field in bright bins with less than three objects equal to the value in the brightest bin with more than three objects (i.e., if there are no magnitude bins with more than three objects at M $\lt \;-$21, the effective volumes for all brighter bins are set equal to the value at $M=-21$). This change has no discernable effect on our luminosity function results, as this is well above the 90% completeness limit for any of our fields and is thus only done to keep small numbers of galaxies from significantly affecting the volumes.

Another possible issue is the source density of simulated objects. In the smaller fields (i.e., HUDF, HUDF parallels, and HFF parallels) we input sources with twice the surface density as in the larger fields to speed up the computing time. As some sources (just like real galaxies) will inevitably fall on top of real sources, and thus not be recovered, an increased source density could result in a (slightly) lower completeness. This is just what is observed in these fields, as shown in Figures 7 and 9. To account for this uncertainty, we measured the spread in volume per unit area in each field at ${M}_{1500}=-21$ at each redshift, which we found to be ∼1.5%, 3.8%, 6.2%, 7.8%, and 13% at z = 4, 5, 6, 7, and 8, respectively. At each step in the MCMC chain, we perturbed the effective volume by this amount to account for this systematic uncertainty in our luminosity function results.

We also examined a non-parametric approach to studying evolution in the luminosity function. This is particularly warranted at very high redshift, where the effects responsible for suppressing the bright end of the luminosity function and causing the exponential decline in number density (e.g., AGN feedback or dust attenuation) may be less relevant. We thus calculated the SWML luminosity function, which is essentially the number density at a given magnitude, free from assumptions about the functional form of number density with magnitude. We also calculated the SWML luminosity function using an MCMC sampler. In this case, since the number densities in the magnitude bins are not linked by an overarching function, we calculated the number density in each magnitude bin independently.

For each magnitude bin and for each field, the likelihood was calculated using Equations (5) and (6) that a given randomly drawn value of $\varphi (M)$ will give the observed number of galaxies. The actual recorded value of $\varphi (M)$ is that which maximizes the likelihood. While in practice, this yields very similar results as one would get by simply taking the observed number and dividing by the effective volume (consistent within a few percent for bins with more than a few galaxies), our approach has two advantages. First, in the limits where numbers are small, this approach is more accurate because it properly accounts for the Poissonian likelihood. Second, this approach generates a full probability distribution for the number densities in each magnitude bin, allowing for the derivation of accurate asymmetric uncertainties. Our SWML luminosity function determinations and best-fit Schechter functions are given in Table 5 and shown in Figure 10.

Zoom In Zoom Out Reset image size

Table 5.  Rest-frame Ultraviolet Luminosity Function: Stepwise Maximum Likelihood Method

M1500	$\varphi \;(z\;\approx$ 4)	$\varphi \;(z\;\approx$ 5)	$\varphi \;(z\;\approx$ 6)	$\varphi \;(z\;\approx$ 7)	$\varphi \;(z\;\approx$ 8)
	(10−3 Mpc−3 mag−1)	(10−3 Mpc−3 mag−1)	(10−3 Mpc−3 mag−1)	(10−3 Mpc−3 mag−1)	(10−3 Mpc−3 mag−1)
−23.0	<0.0016	<0.0023	<0.0025	<0.0029	<0.0035
−22.5	0.0093${}_{-0.0033}^{+0.0045}$	0.0082${}_{-0.0035}^{+0.0050}$	<0.0025	<0.0029	<0.0035
−22.0	0.0276${}_{-0.0062}^{+0.0074}$	0.0082${}_{-0.0036}^{+0.0051}$	0.0091${}_{-0.0039}^{+0.0057}$	0.0046${}_{-0.0028}^{+0.0049}$	<0.0035
−21.5	0.1192${}_{-0.0132}^{+0.0145}$	0.0758${}_{-0.0125}^{+0.0137}$	0.0338${}_{-0.0085}^{+0.0105}$	0.0187${}_{-0.0067}^{+0.0085}$	0.0079${}_{-0.0046}^{+0.0068}$
−21.0	0.2968${}_{-0.0219}^{+0.0230}$	0.2564${}_{-0.0240}^{+0.0255}$	0.0703${}_{-0.0128}^{+0.0148}$	0.0690${}_{-0.0144}^{+0.0156}$	0.0150${}_{-0.0070}^{+0.0094}$
−20.5	0.6491${}_{-0.0347}^{+0.0361}$	0.5181${}_{-0.0338}^{+0.0365}$	0.1910${}_{-0.0229}^{+0.0249}$	0.1301${}_{-0.0200}^{+0.0239}$	0.0615${}_{-0.0165}^{+0.0197}$
−20.0	1.2637${}_{-0.0474}^{+0.0494}$	0.9315${}_{-0.0482}^{+0.0477}$	0.3970${}_{-0.0357}^{+0.0394}$	0.2742${}_{-0.0329}^{+0.0379}$	0.1097${}_{-0.0309}^{+0.0356}$
−19.5	1.6645${}_{-0.0618}^{+0.0630}$	1.2086${}_{-0.0666}^{+0.0488}$	0.5858${}_{-0.0437}^{+0.0527}$	0.3848${}_{-0.0586}^{+0.0633}$	0.2174${}_{-0.1250}^{+0.1805}$
−19.0	2.6392${}_{-0.1165}^{+0.1192}$	2.0874${}_{-0.1147}^{+0.1212}$	0.8375${}_{-0.0824}^{+0.0916}$	0.5699${}_{-0.1817}^{+0.2229}$	0.6073${}_{-0.2616}^{+0.3501}$
−18.5	3.6169${}_{-0.6091}^{+0.6799}$	3.6886${}_{-0.3725}^{+0.3864}$	2.4450${}_{-0.3515}^{+0.3887}$	2.5650${}_{-0.7161}^{+0.8735}$	1.5110${}_{-0.7718}^{+1.0726}$
−18.0	5.8343${}_{-0.8204}^{+0.8836}$	4.7361${}_{-0.4413}^{+0.4823}$	3.6662${}_{-0.8401}^{+1.0076}$	3.0780${}_{-0.8845}^{+1.0837}$	⋯
−17.5	6.4858${}_{-0.9467}^{+1.0166}$	7.0842${}_{-1.1364}^{+1.2829}$	5.9126${}_{-1.2338}^{+1.4481}$	⋯	⋯

Note. Magnitude bins with zero objects are shown as upper limits, calculated as 1/V_eff/ΔM in that magnitude bin for that redshift.

Download table as:  ASCIITypeset image

As shown in Figure 10, the qualitative shape of the SWML luminosity functions at all redshifts we consider here are similar, in that bright galaxies are rare and faint galaxies are relatively common. Additionally, when examining the Schechter fits (solid line), we see that they are consistent with the SWML determinations. Surprisingly, the best-fit Schechter function parameters (Table 4) show little evolution in ${M}_{\mathrm{UV}}^{*}$. However, from z = 4 to 8, the uncertainty on ${M}_{\mathrm{UV}}^{*}$ gets progressively larger, to 0.4 (0.9) mag at z = 7 (8). This is easy to understand, as at all redshifts, our data set contains galaxies in only 1–2 bins brightward of ${M}_{\mathrm{UV}}^{*}$. Ideally, one would prefer to have multiple bins in magnitude on either side of ${M}_{\mathrm{UV}}^{*}$ to obtain robust constraints. However, that requires a larger volume than we consider in this analysis. In Figure 11, we fit the evolution of ${M}_{\mathrm{UV}}^{*}$ with redshift with a linear function, using our results at z = 4–8. We find that ${dM}*/{dz}=-$0.12 ± 0.09; thus, our data do not support a significant evolution of ${M}_{\mathrm{UV}}^{*}$ with redshift.

Zoom In Zoom Out Reset image size

We also fit similar functions to see if we detect evolution in α and $\varphi *$. As shown in Figure 11, we do see significant evolution in the faint-end slope α, with it becoming steeper at higher redshift, as $d\alpha /{dz}=-$0.19 ± 0.04 (4.8σ significance). We see a similar significance in the evolution of the characteristic number density $\varphi *$, which evolves as $d\mathrm{log}\varphi */{dz}=-$0.31 ± 0.07 (4.4σ significance). Thus, while ${M}_{\mathrm{UV}}^{*}$ does not significantly evolve with redshift from z = 4 to 8, both α and $\varphi *$ do, in that the number density decreases and the faint-end slope becomes steeper with increasing redshift. In particular, this decline in characteristic number density is by a factor of ∼20 over a period of time of less than 1 Gyr. Although the steepening of the faint end is consistent with previous studies (e.g., Bouwens et al. 2012), the un-evolving ${M}_{\mathrm{UV}}^{*}$ and strong number density evolution are the opposite of what was presented in the literature just a year ago (e.g., Bouwens et al. 2007, 2011b; McLure et al. 2013). This updated evolutionary picture will be crucial when projecting number counts for future HST and JWST surveys.

In Figure 12, we show our determinations of the luminosity functions together at all five redshifts, along with the joint confidence contours on ${M}_{\mathrm{UV}}^{*}$, α, and $\varphi *$. It is apparent that there is significant evolution in the luminosity function, with a drop in number density from z = 4 to 8, as well as a gradual steepening of the faint-end slope. The apparent non-evolution of the characteristic magnitude is visible as the roughly constant magnitude of the "knee" of the luminosity function.

Zoom In Zoom Out Reset image size

By definition, our method of computing the luminosity function is dependent on magnitude binning, because we compare the observed number to that expected based on a given model in magnitude bins of width 0.5 mag. While galaxies close to one side of a magnitude bin have the potential to scatter to another bin, the typical uncertainties on the UV absolute magnitudes of galaxies in our sample are ∼0.2 mag. Additionally, galaxies can shift both ways; thus, while one galaxy moves out of a bin, another may move in (although this effect will not be symmetric given the shape of the luminosity function). In our results in the previous section, we assumed that magnitude scatter does not significantly impact our results.

To investigate the impact of this assumption, we performed another iteration of MCMC fitting to our data, here allowing galaxies to scatter between magnitude bins. At each step in the MCMC chain, a new value of M₁₅₀₀ was drawn for each galaxy from the 100 SED-fitting Monte Carlo simulation results. The spread in these values encompassed both the photometric scatter in the observed filters and the uncertainty in the photometric redshift (see Section 3.7). To compare to our fiducial luminosity function values, we recorded both the median Schechter parameter results and the median number density in each magnitude bin, because, unlike our fiducial MCMC run, these varied during each step as the magnitudes changed. At all redshifts, our fiducial values of the stepwise luminosity functions are consistent with these "magnitude-scatter" values within 10% at M $\geqslant -$21.5, and typically within 2%–3%. The sole exception is in the brightest bin (−22 at z = 4–6, and −21.5 at z = 7 and 8), where our fiducial number density values are higher by ∼15%–20% (60% at z = 7, where there is only a single galaxy in this bin). We examined the Schechter fit to see whether this bright-end difference affects our results. Values of both ${M}_{\mathrm{UV}}^{*}$ and α derived when allowing galaxies to shift between bins are consistent with our fiducial values within 0.1 mag and <3%, respectively. We conclude that the relatively small (∼20%) uncertainties in the absolute magnitudes of our galaxies do not have a significant impact on our luminosity function results.

Given that our results show that the Schechter functional parameters may not be a robust method of tracking galaxy evolution (e.g., a non-evolving value of ${M}_{\mathrm{UV}}^{*}$ does not mean that the galaxy populations are not evolving), we examine the evolution in a non-parametric way. In Figure 13 we show the evolution of the stepwise luminosity function, plotting the number density corresponding to galaxies at ${M}_{\mathrm{UV}}=-21$ and −19 versus redshift. From z = 8 to 4, the abundance of brighter galaxies increases faster than faint galaxies. This trend halts at z = 4, where bright galaxies have an approximately constant abundance down to z = 2, and then turn over. Faint galaxies, however, continue increasing in abundance down to z = 2, where they also turn over. This figure highlights the phenomenon of downsizing, where bright/large galaxies grow faster at early times (e.g., Cowie et al. 1996, see also Lundgren et al. 2014). This is different from the expectation one would get simply from examining Schechter fits, because the luminosity functions don't evolve much over the range 2 $\lt \;z\;\lt$ 4 (e.g., Reddy & Steidel 2009). Given that the trends here mimic the evolution of the cosmic SFR density, we fit the function provided by Madau & Dickinson (2014) to our data for both number densities, given by

$$\begin{eqnarray}&&\varphi (z)=A\;\displaystyle \frac{{(1+z)}^{\alpha }}{1+{[(1+z)/B]}^{\gamma }}\;{\mathrm{mag}}^{-1}\;{\mathrm{Mpc}}^{-3}.\end{eqnarray} \tag{ 7 }$$

The evolution with redshift is thus proportional to ($1+z$)${}^{\alpha }$ at low redshift, and ($1+z$)${}^{\alpha -\gamma }$ at high redshift, with A and B as dimensionless coefficients. Fitting all of the data, including the literature data shown in Figure 13, in this way, we confirm that at $z\;\gt$ 3, bright galaxies change in abundance faster, as ($1+z$)${}^{-5.9\pm 0.4}$, than faint galaxies, which go as ($1+z$)${}^{-3.3\pm 0.3}$.

Zoom In Zoom Out Reset image size

Another interesting aspect is to compare the trends observed to our predicted abundance of bright z = 9 galaxies (see Section 9). The trend observed here slightly overestimates our predicted z = 9 abundance, although if we assume the uncertainties on our z = 8 number density applies to z = 9, our trend is consistent with this prediction. In any case, this trend of abundance with redshift lends more weight to our expectation of a significant abundance of bright z = 9 galaxies. This figure neglects the impact of dust attenuation, because we are only looking at the observed UV magnitudes. The dust attenuation appears to be luminosity dependent (bright galaxies are dustier than faint galaxies; e.g., Bouwens et al. 2014), as well as being higher at lower redshift. So, correcting for dust would not only increase the abundance of bright galaxies more than that of faint galaxies, it would increase it by more at lower redshift, thus enhancing the differences between faint and bright galaxies at $z\;\gt$ 4.

Our result of a similarly bright value of ${M}_{\mathrm{UV}}^{*}$ at z = 6, 7, and 8 is a dramatic change from previously published results. In Figure 10, we show the stepwise luminosity function results from several relevant studies from the literature. Figure 15 shows our uncertainty results, highlighting both the distribution of ${M}_{\mathrm{UV}}^{*}$ and α from the MCMC chains, as well as the covariance between the two parameters, along with previous determinations of ${M}_{\mathrm{UV}}^{*}$ and α (Bouwens et al. 2007, 2011b; McLure et al. 2009, 2013; Ouchi et al. 2009; van der Burg et al. 2010; Schenker et al. 2013; Willott et al. 2013; Bowler et al. 2014). In this subsection we compare solely to previous work—we reserve the comparison to the contemporaneous work by Bouwens et al. (2015) to Section 6.4.1 below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At z = 4 and z = 5, both our binned luminosity function data points, as well as our Schechter function parameters, are in excellent agreement with the ground-based study of van der Burg et al. (2010). We are also in excellent agreement with the ground-based study of McLure et al. (2009) at z = 5. We find good agreement with the space-based study of Bouwens et al. (2007) at z = 5, but at z = 4 the Bouwens et al. (2007) result lies outside our 2σ confidence region on the Schechter function parameters, in that we prefer a shallower faint-end slope and a fainter value for ${M}_{\mathrm{UV}}^{*}$.

At z = 6, our binned luminosity function data points are consistent within 1-2σ with the Bouwens et al. (2007) results at the faint end (Figure 10). At the bright end, our data are higher than those from both ground-based studies (although typically only different at the 1-2σ level). This is somewhat counterintuitive, as one may expect the ground-based studies to suffer a higher contamination rate, particularly for relatively fainter sources at higher redshift, due to their inability to resolve stars from galaxies; however, it may also be explained due to an aggressive sample selection required to minimize contamination. In any case, the differences are not highly significant, with the exception of the brightest data point from Willott et al. (2013), which gives a number density at $M=-22.5$ of 2.7 × 10⁻⁸ Mpc⁻³. While this is consistent with our upper limit at that magnitude, it is a factor of ∼250 lower than our number density only 0.5 mag fainter at $M=-22$ (see Table 5). Given the results at similar magnitudes at lower redshifts, it is highly unlikely that there is such a steep drop in number density over only a 0.5 mag interval, although future large area studies can better investigate the difference (Bowler et al. 2015). The larger discrepancy comes when comparing the Schechter function parameters. Specifically, Bouwens et al. (2007) found $M*=-$20.29 ± 0.19, and McLure et al. (2009) found $M*=-$20.04 ± 0.12. Both values are significantly (2–3σ) fainter than our derivation of $M*=-$21.13${}_{-0.31}^{+0.25}$. For the space-based study of Bouwens et al. (2007), this is understandable, because at that time only optical data were available; thus, z = 6 galaxies were selected via detections in only one band, and a robust determination of their UV absolute magnitudes was difficult. For the ground-based study of McLure et al. (2009), a cause for the difference is less clear, although the different data being used certainly plays a role.

Comparing to previous works at z = 7, we find broadly similar results, in that our results are consistent with the derived number densities from previous studies, yet our Schechter fit prefers a much brighter value of ${M}_{\mathrm{UV}}^{*}$. This is easier to understand, because a number of previous studies had less data available, and thus, utilizing smaller volumes, were unable to constrain the bright end (e.g., Bouwens et al. 2011b; Schenker et al. 2013). The exception is the recent work by McLure et al. (2013), which used a similar volume as our study, except they used the CANDELS UDS field rather than the CANDELS GOODS-N field. Examining the brightest data point from McLure et al. (2013) at $M=-$21, the number density is about a factor of two below our data point. However, the discrepancy is mitigated by two factors. First, as discussed by Bouwens et al. (2015), the use of fixed diameter circular apertures by McLure et al. (2013) systematically underestimates the fluxes for bright, more extended, galaxies. Bouwens et al. (2015) estimated the amplitude of this effect to be ∼0.25 mag. Shifting the brightest McLure et al. (2013) data point by 0.25 mag brings it into agreement with our results. Second, the CANDELS GOODS-N field appears to have an overdensity of z = 7 galaxies. Specifically, when comparing the number density of z = 7 galaxies in the GOODS-N Deep field to the GOODS-S Deep field in Figure 15, GOODS-N has a higher number density at all magnitudes. While we have not selected galaxy samples in the UDS, we can examine this further by recomputing our z = 7 stepwise luminosity function using only the GOODS-S and HUDF fields. At magnitudes fainter than −21, the results do not change appreciably, as the GOODS-N Deep and Wide fields lie on either side of our Schechter fit at those magnitudes. However, the results using only GOODS-S provide a number density ∼33% lower at $M=-21.5$ than our fiducial luminosity function. This difference is at the 1σ level, due to the large Poisson noise contribution in this bin, and thus is not highly significant.

We also compare to several ground-based studies at z = 7. Ouchi et al. (2009) identified 22 bright $z\;\sim$ 7 candidate galaxies over ∼0.4 deg². Their data points based on detected galaxies are consistent with our own, though their strict upper limits at $M\sim -22$ push their Schechter fit to a fainter value of ${M}_{\mathrm{UV}}^{*}=-$20.1, although the large uncertainty (0.76) leaves ${M}_{\mathrm{UV}}^{*}$ consistent with our fit at only slightly more than the 1σ level. The stepwise luminosity function from Castellano et al. (2010) based on deep HAWK-I data agrees well with our results, while the results from the zFourGE medium band survey of Tilvi et al. (2013) agree at $M=-$21.5, but differ by ∼2σ at $M=-20.5$.

Recently, Bowler et al. (2014) made a significant improvement in search volume from the ground, discovering 34 luminous $z\;\sim$ 7 galaxy candidates over 1.65 deg², from the UltraVISTA survey data over the COSMOS field (McCracken et al. 2012) and the UKIDSS survey over the UDS field (Lawrence et al. 2007). Broadly speaking, they are consistent with our results, and are highly inconsistent with the previous determinations of ${M}_{\mathrm{UV}}^{*}\sim -20$ (Figure 10). There is a mild tension at $M=-$21.75, where the value of our Schechter fit at that point is 2σ higher than their derived number density. However, this is their faintest magnitude bin, and is only ∼50% complete, thus this data point relies the most on the completeness correction. In any case, the fact that Bowler et al. (2014) found $z\;\sim$ 7 candidates out to very bright magnitudes gives us confidence that our brighter determination of ${M}_{\mathrm{UV}}^{*}$ is not necessarily dominated by cosmic variance in our fields, but is a true feature of the z = 7 universe. However, this uncertainty on ${M}_{\mathrm{UV}}^{*}$ of ∼0.4 mag makes it apparent that more data are needed to constrain this parameter further.

At z = 8, we again find consistent number densities with previous studies, although our larger volume allows us to find more rare, bright ($M=-$21.5) galaxies than observed in some previous surveys, pushing them to lower values of ${M}_{\mathrm{UV}}^{*}$ (although our uncertainty on ${M}_{\mathrm{UV}}^{*}$ is large, so the difference in our determination is not significantly different from previous studies). As noted, in our fit of the z = 8 luminosity function we constrained ${M}_{\mathrm{UV}}^{*}$ to be fainter than −23, to avoid un-physically bright values, which tended to be preferred in an unconstrained fit. We note two important points when comparing to previous studies. First, while our study did not utilize the pure parallel BoRG (Trenti et al. 2011) and HIPPIES (Yan et al. 2011) programs, our bright end is consistent with that from Schmidt et al. (2014), based on a determination of the z = 8 luminosity function over 350 arcmin² of pure parallel data (for comparison, our search area at z = 8 comprised ∼300 arcmin²; Table 1). The multiple sight-lines of BoRG and HIPPIES leave their results less susceptible to cosmic variance effects, so the agreement implies that cosmic variance may not be strongly affecting our bright end, although we explore this in Section 6.5.

A potentially larger difference between our results and those of previous studies is also seen at the faint end, in that our faint-end slope is possibly steeper than previously found. However, our uncertainty is large, such that our result of $\alpha =-$2.36${}_{-0.40}^{+0.54}$ is consistent with previous results of $\alpha \approx -$2 (e.g., Bouwens et al. 2011b; McLure et al. 2013; Schenker et al. 2013). Previous studies use galaxies as faint as M $=-17.5$ in their determination of the faint-end slope at z = 8. As discussed, and shown in Figure 8, we find that we fall below 50% completeness at $M\gt -$18.5, and thus we do not use galaxies fainter than that in our luminosity function determinations. While robust estimates of the number densities of galaxies at −18.5 ≤ M $\leqslant \;-17.5$ would certainly improve the confidence on the faint-end slope, we use the same deep data sets as the other referenced studies (HUDF). We would expect the incompleteness to be similar between all studies, although it does depend on sample selection and the exact details of the incompleteness simulations. In any case, constraints on the faint-end slope at z = 8 should improve in the near future with further data from the HFF program. Our inclusion of the HFF parallel imaging, even contributing only a small number of galaxies at $z\;\gt$ 7, did improve the fractional error on the faint-end slope (${\sigma }_{\alpha }/\alpha$) by 2% and 8% at z = 7 and z = 8, respectively.

Finally, there have also been theoretical estimates of the luminosity functions at these redshifts, most prominently from Jaacks et al. (2012a), who made predictions in good agreement with our observed luminosity functions. Specifically, their simulations also predict bright values of ${M}_{\mathrm{UV}}^{*}$ of −21.15, −20.82, and −21.00 at z = 6, 7, and 8, respectively. They also found quite steep faint-end slopes of −2.15${}_{-0.15}^{+0.24}$, −2.30${}_{-0.18}^{+0.28}$, and −2.51${}_{-0.17}^{+0.27}$ at z = 6, 7, and 8, respectively. Within the uncertainties, these faint-end slopes are consistent with our own, although the apparent agreement at z = 8 is tantalizing (as mentioned above, we cannot constrain the slope to be so steep). Steep faint-end slopes of $\alpha \sim -$2 at these redshifts were also seen by Salvaterra et al. (2011) and Dayal et al. (2013), but both studies also predict a brightening in ${M}_{\mathrm{UV}}^{*}$ toward lower redshift, which is contrary to our observations.

6.4.1. Comparison to Bouwens et al. (2015)

6.4.2. Previously Published Measurement Uncertainties

The impact of cosmic variance on our measurement of the luminosity function is minimized due to our use of several fields, which are split into four widely separated regions of the sky. However, as shown in Figure 15, there is significant variance between the different fields, particularly at $z\;\geqslant$ 5. To estimate the effect of cosmic variance on our derived number densities, we used two techniques: a semi-empirical technique, combining the QUICKCV calculator provided by Newman & Davis (2002; using the updated version provided by Moster et al. 2011) with the recent clustering-based bias measurements from Barone-Nugent et al. (2014); and a semi-analytic technique, based on the semi-analytic models (SAMs) of Somerville et al. (2012). In Section 7 we discuss how these models, with a modification to the redshift dependance of the normalization of the dust attenuation, provide an excellent match to our measured UV luminosity functions.

For a given survey geometry, the QUICKCV program returns the fractional error in a count due to cosmic variance for an unbiased tracer of matter at a given redshift. For this calculation, we estimated the fractional error separately for GOODS-S, GOODS-N, MACS-0416 parallel, and Abell 2744 parallel fields, adding the variances in quadrature to derive a final value of ${\sigma }_{\mathrm{CV}}$ for a given redshift bin. In the GOODS-S field, we included the area from the three HUDF09 fields, because even the parallel fields are separated by only a few arcmin from the GOODS-S proper. For the input survey geometries, we estimate rectangular regions of the approximate shape of the GOODS fields, with an enclosed area equal to the GOODS-S Deep+Wide+ERS+HUDF09 fields for GOODS-S, and the GOODS-N Deep+Wide for GOODS-N. The field geometries were thus $10\buildrel{\,\prime}\over{.} 2$ × 1503 for GOODS-S, $9\buildrel{\prime\prime}\over{.} 51$× 1465 for GOODS-N, and $2\buildrel{\,\prime}\over{.} 1$ × 21 for each of the HFF parallel fields. However, at faint magnitudes, our galaxy sample primarily comes from the HUDF, thus we also estimated the QUICKCV-derived cosmic variance uncertainty with the HUDF area only, with a geometry of $2\buildrel{\,\prime}\over{.} 26$ × 226. To convert these unbiased estimates of cosmic variance to values appropriate to our galaxy sample, we use the recently published clustering-based bias measurements from Barone-Nugent et al. (2014), who used the galaxy sample of Bouwens et al. (2015) for their calculation (because they did not measure the clustering at z = 8, we use the z = 7 bias values for our z = 8 cosmic variance estimate). They estimated the bias for both bright and faint galaxies, splitting their sample at ${M}_{\mathrm{UV}}=-19.4$ at each redshift. Our estimates of the fractional uncertainty on galaxy counts due to cosmic variance from this method are shown as the gray bars in Figure 17, where we show values of this quantity for both bright and faint galaxies.

Zoom In Zoom Out Reset image size

For our semi-analytic cosmic variance estimate, we used mock catalogs of the Somerville et al. (2012) SAMs, which cover an area ∼40× larger than that of the combined CANDELS/GOODS fields. We extract independant, GOODS-sided volumes from these catalogs, exploring the variation in number counts in the independant volumes as a function of UV absolute magnitude. At magnitudes brighter than −18.5 at z = 4–6 (−19 at z = 7; −19.5 at z = 8) we estimated our survey as being two $10^{\prime} \times 16^{\prime}$ fields, representing a combination of the CANDELS/GOODS fields with the five single-WFC3 pointing fields (HUDF09 and HFF). At fainter magnitudes, where the majority of our objects come from the HUDF, we assume a single 226 × 226 field. We calculate the 1σ fractional uncertainty on the number density, ${\sigma }_{\mathrm{cv}}/N$, by bootstrap resampling galaxies in a given M_UV bin at each redshift. This 1-sigma fractional uncertainty includes the Poisson noise, thus we subtract the Poisson errors using the recipe of Gehrels (1986) to calculate the fractional uncertainty on the number density due to cosmic variance only. The uncertainty for the total survey volume is calculated by adding the variance for two GOODS-sized fields in quadrature for bright bins (and the HUDF-only for faint bins), and is shown in Figure 17.

Comparing the SAM-derived cosmic variance values to those from QUICKCV, we find generally excellent agreement. The SAM method predicts, as expected, a larger cosmic variance uncertainty for the brightest galaxies, although this is understood as the Barone-Nugent et al. (2014) "bright" sample encompassed galaxies down to ${M}_{\mathrm{UV}}=-$19.4, and thus likely has a relatively faint median magnitude. Future measurements of the bias in finer-resolution bins of bright galaxies can probe this effect further. The agreement at the faint end is generally good as well, with the exception of z = 6, where the semi-empirical method predicts a somewhat low uncertainty due to the Barone-Nugent et al. (2014) bias measure at z = 6, which is slightly lower than at z = 5. While both of these methods are estimates, and thus may not be extremely accurate, the general agreement between these two different techniques implies that our estimates of the cosmic variance uncertainty are not highly inaccurate.

To assess the impact of cosmic variance on our measured luminosity functions, we compared both of our estimates of the cosmic variance uncertainty to the Poisson noise from our stepwise luminosity functions, as shown in Figure 17. At all redshifts, the data at M $\lt \;-21$ are dominated by Poisson noise (M $\lt \;-$20.5 at $z\;\geqslant$ 6), thus we do not expect cosmic variance to be dominating the uncertainties on the bright end of the luminosity functions derived here. However, cosmic variance may play some role at the faint end, where we are restricted to small fields. At z = 7, there does appear to be a step in the stepwise luminosity function at $M\geqslant -$18.5, where the $M=-19$ point is below our best-fit Schechter function, and the $M=-18.5$ point is above. At M $\geqslant \;-18.5$ our data come from only the HUDF main field, thus this break represents the point where we become reliant on a single small field. As shown in Figure 17, the cosmic variance uncertainty at faint magnitudes at z = 7 is ∼40%, which is somewhat larger than the Poisson uncertainty; thus, cosmic variance may bear some responsibility for this discontinuity in the luminosity function at the faint end. Future measures of the luminosity function at $M\geqslant -18.5$ from the HFF lensing program may improve these constraints, but while faint galaxies may be found, the volumes will still be incredibly small (e.g., Robertson et al. 2014). Thus, robust constraints on the number densities at this faint level at $z\;\geqslant$ 7 may need to wait until the JWST. We conclude that while the effects of cosmic variance are not negligible, they are not the dominant source of uncertainty on the abundances of the bright objects we have discovered at very high redshift.

When allowed to choose any value of ${M}_{\mathrm{UV}}^{*}$, our z = 8 Schechter function fit preferred very bright values of ${M}_{\mathrm{UV}}^{*}$, such that all observed data points lay on the faint-end slope part of the function. This implies that the z = 8 luminosity function is consistent with a single power law. Such a functional form is what one might expect when the feedback effects that govern the bright end at lower redshift (i.e., mainly feedback due to accreting supermassive black holes) disappear, or if dust attenuation ceases to be a factor. Bowler et al. (2014) recently postulated that the z = 7 luminosity function is better fit by a double power law, rather than a Schechter form. At z = 7, our stepwise data appear consistent with the Schechter fit out to the brightest magnitudes we cover. To see whether our data show a preference for a Schechter functional form at all redshifts, we performed three fits to the data: a Schechter fit, a single power law, and a double power law. To place these fits on equal ground, we found the best-fit parameters for each function using a simple maximum likelihood routine. For the Schechter fit, we used the function shown in Equation (4), investigating a range of ${M}_{\mathrm{UV}}^{*}$ with ${\rm{\Delta }}M\;=$ 0.1 mag, and α with ${\rm{\Delta }}\alpha \;=$ 0.02. We approximated a single power law using the Schechter functional form with ${M}_{\mathrm{UV}}^{*}$ fixed at −30. For the double power law, we used the form given in Equation (2) of Bowler et al. (2014), which is similar to the Schechter function at the faint end, but replaces the bright end with a second power law with slope β. In all cases, $\varphi *$ is found as the normalization, such that the total number of expected objects for a given function is equal to the number observed. The likelihood that a given functional form represents our data was calculated in an identical manner as in Section 5.1, using Equations (5) and (6).

To compare the results from these fits at each redshift, we used the Bayesian information criterion (BIC). This is similar to a ${\chi }^{2}$ statistic, except that it takes into account both the number of data points and the number of free parameters. For a model to be preferred over a competing model, it must have a BIC lower by at least 2. This is sensible, because adding a free parameter must yield a better fit for that model to be preferred. The BIC is calculated as

$$\begin{eqnarray}&&\mathrm{BIC}=-2\ \mathrm{ln}({\mathcal{L}})+k\ \mathrm{ln}(N)\end{eqnarray} \tag{ 8 }$$

where N is the number of data points and k is the number of free parameters (Liddle 2004). For the Schechter, double power law, and single power law fits, the number of free parameters are 2 (${M}_{\mathrm{UV}}^{*}$, α), 3 (${M}_{\mathrm{UV}}^{*}$, α, β), and 1 (α), respectively (we do not count $\varphi *$ as a free parameter because it is a normalization). The number of data points is the number of galaxies from our sample used in the fit, which is restricted to those brighter than the 50% completeness limits discussed previously. This gives N = 2788, 1812, 605, 221, and 47 galaxies at z = 4, 5, 6, 7, and 8, respectively.

The results from this analysis are shown in Table 6. A difference in the absolute value of the BIC of 2 is interpreted as "positive" evidence, whereas a difference of 6 or higher is "strong" evidence, both in favor of the model with the smaller value. In Table 6, in addition to the value of the BIC, we show the difference between the BIC values for the Schechter versus double power law, and Schechter versus single power law. In this formalism, a negative difference is in favor of the Schechter function. Comparing the Schechter function versus the double power law, we find that a Schechter form is strongly preferred to either a double or single power law at z = 4–7. This is not surprising; there is clearly a deficit of observed galaxies at the bright end when compared to the best-fit power law (Figure 10). However, no such deficit is visible at z = 8, and this is confirmed, as both the Schechter fit and the single power law fit have effectively identical values as the BIC. We conclude that our data support an exponential decline at the bright end of the luminosity function at z = 4, 5, 6, and 7. At z = 8, we do not see any evidence for a decline in the bright end, at least out to $M=-$21.5. Further data are needed to show whether one can either detect, or rule out, a decline at the bright end at z = 8. If the latter is true, it could indicate a significant change in the halo masses of bright galaxies, a drop in dust attenuation in bright galaxies, or a change in the physics governing the feedback in bright galaxies in the distant universe.

Although it is presently generally assumed that galaxies dominated the ionizing photon budget for the reionization of the IGM, it has been difficult for observations to obtain robust proof. Analyses of the IGM via line-of-sight quasar observations have been able to show that reionization was likely complete by $z\;\sim$ 6 (e.g., Fan et al. 2006; although see Mesinger 2010 and Becker et al. 2015). Additionally, observations of the cosmic microwave background (CMB) radiation constrain the total optical depth due to electron scattering, which, while it cannot directly inform the duration of reionization, can give an estimate of the reionization redshift (z_reion) if reionization was instantaneous. The results from the WMAP nine-year data set give ${\tau }_{\mathrm{es}}\;=$ 0.088 ± 0.014, which corresponds to ${z}_{\mathrm{reion}}\;=$ 10.6 ± 1.2 (Hinshaw et al. 2013), while the recent Planck results suggest ${z}_{\mathrm{reion}}\;=$ ${8.8}_{-1.1}^{+1.2}$ (Planck Collaboration et al. 2015).

The primary reason for the current uncertainties in the contribution of galaxies to reionization lies in the uncertainty in the faint-end slope measurements, as well as in the assumptions of the escape fraction of ionizing photons (f_esc) and the clumping factor in the IGM (C). The clumping factor is primarily constrained theoretically, but most studies agree that it is low (<6) at high redshift (e.g., Faucher-Giguère et al. 2008; Pawlik et al. 2009; McQuinn et al. 2011; Finlator et al. 2012). To infer a number of escaping ionizing photons from observations of galaxies, one first needs to take the observed UV light and assume an IMF and a metallicity. Then, to calculate the number of these ionizing photons available for reionization, multiply by an assumed value of f_esc. However, it is difficult to constrain f_esc directly at high redshift, because the correction for intervening IGM absorption systems is extremely high at $z\;\gt$ 4. Significant effort is being expended on observationally constraining f_esc at $z\;\lt$ 4. Although bright galaxies at $z\;\sim$ 1 have very low escape fractions (relative to the UV emission) of ${f}_{\mathrm{esc},\mathrm{rel}}\lt$ 2% (Siana et al. 2010), escaping ionizing emission has been observed from small fractions of galaxies probed by studies at $z\;\sim$ 3–4 (e.g., Steidel et al. 2001; Shapley et al. 2006; Iwata et al. 2009; Vanzella et al. 2010; Nestor et al. 2011); though some ground-based studies may suffer from contamination by intervening sources (e.g., Vanzella et al. 2012). Recent results imply that escape fractions from star-forming galaxies at $z\;\sim$ 2–3 range from 5% to 20%, with lower-mass galaxies, especially those with Lyα in emission, having a greater likelihood of having detectable escaping ionizing emission (e.g., Mostardi et al. 2013; Nestor et al. 2013).

Finkelstein et al. (2012a) used measurements of the emission rate of ionizing photons from observations of the Lyα forest in quasar spectra to place an upper limit on f_esc from galaxies. Assuming that the rest-frame UV luminosity function extended down to ${M}_{\mathrm{UV}}=-$13, the escape fraction must be ${f}_{\mathrm{esc}}\;\lt$ 13% to avoid violating the Lyα forest measurements of Bolton & Haehnelt (2007). Using this value, and assuming C = 3, the luminosity functions available at the time were consistent with a wide range of reionization histories, including an end redshift as late as $z\;\lesssim$ 5, and an ionized fraction at $z\;\sim$ 7 from 30% to 100%. Kuhlen & Faucher-Giguère (2012) and Robertson et al. (2013) did similar analyses, folding in additional observables (e.g., the Lyα forest and CMB), and found that in order to complete reionization by $z\;\sim$ 6, the luminosity function must extend much deeper than can presently be observed, and/or the average escape fraction must be higher at higher redshift.

Here, we use our updated luminosity functions to re-examine the contribution of galaxies to reionization. Figure 20 shows both the observable specific UV luminosity density (${\rho }_{\mathrm{UV}}$), which we define to be that above our 50% completeness limit, as well as the total ${\rho }_{\mathrm{UV}}$, which we define as the integrated luminosity function down to ${M}_{1500}=-$13. We then compare these values to the critical number of UV photons necessary to sustain an ionized IGM at a given redshift, taken from Madau et al. (1999). This figure is similar to Figure 3 from Finkelstein et al. (2012a), thus we refer the reader there for more details. Effectively, these critical curves depend on assumptions about the stellar IMF, metallicity, f_esc, and clumping factor. The first two are responsible for the conversion from observed UV photons to intrinsic ionizing photons. We assumed a Salpeter IMF, and the width of the curves denote the impact of changing the metallicity from 0.2 $\leqslant \;Z$/Z${}_{\odot }$ ≤ 1.0. We show several curves for the reader's choice of the ratio of $C/{f}_{\mathrm{esc}}$. Here, we use a fiducial value of C = 3 and ${f}_{\mathrm{esc}}\;=$ 13%, which are consistent with Finkelstein et al. (2012a).

Zoom In Zoom Out Reset image size

The right panel of Figure 20 shows the ionization history of the IGM, comparing our derived value for the total specific UV luminosity density to our fiducial model of C = 3 and ${f}_{\mathrm{esc}}\;=$ 13%, folding in the values at z = 10.4 from Bouwens et al. (2015) to extend our analysis beyond z = 8. Our luminosity functions are consistent with a reionization history that starts at $z\;\sim$ 11, and ends by $z\;\gt$ 5. Although the exact value of the volume ionized fraction in the IGM is uncertain between these redshifts, due to the persistent uncertainty in the faint-end slope, our results imply the following constraints (given the caveat of our assumptions). At z = 6, we can constrain ${x}_{{\rm{H}}\;{\rm{II}}}\;\gt$ 0.85 (1σ), while out to the limit of our observations at z = 8 the data are still consistent with a fully ionized IGM (68% C.L. of 0.15 $\lt \;{x}_{{\rm{H}}\;{\rm{II}}}=\lt$1.0). We find a midpoint of reionization (x${}_{{\rm{H}}\;{\rm{II}}}\;=$ 0.5) of 6.7 $\lt \;z\;\lt$ 9.4 (68% C.L.).

Broadly speaking, measurements from quasar spectra as well as from Lyα emission from galaxies support a reionization scenario consistent with what we derive (Figure 20). The constraints from Lyα emission are heavily model dependent, and studies claiming a very low value of x_{H ii} may be assuming a velocity offset of Lyα from systemic that is too high (e.g., Stark et al. 2014). While our results are in slight tension (1.3σ) with the results from WMAP9, they are in excellent agreement with the more recent results from Planck. From our fiducial reionization history, we find ${\tau }_{\mathrm{es}}\;=$ 0.063 ± 0.013, which is highly consistent with the measurement from Planck of ${\tau }_{\mathrm{es}}\;=$ 0.066 ± 0.012. We remind the reader that our results did not use the Planck results as a constraint; rather, the inferred reionization history from our observed luminosity functions (with our assumptions on C and f_esc) are in remarkable agreement with the Planck observations.

Future observations are necessary to improve the constraints on reionization from galaxies. Specifically, more robust measurements of the faint-end slope α at z = 6–8 can dramatically shrink the uncertainties on ${\rho }_{\mathrm{UV}}$, subsequently reducing the width of our plausible values of x_{H ii}. Likewise, improving the measurements at $z\;\geqslant$ 9 will inform us on whether the ionization fraction of the IGM at that early time was significantly non-zero. The HFF program will improve both of these areas, although definitive results will likely not be obtained until the JWST era.