Bridging Star-forming Galaxy and AGN Ultraviolet Luminosity Functions at z = 4 with the SHELA Wide-field Survey

In this study, we use imaging in nine photometric bands spanning the optical to mid-IR in the SHELA Field. The SHELA Field is centered at R.A. = 1^h22^m00^s, decl. = +0°00'00'' (J2000) and extends approximately ±65 in R.A. and ±125 in decl. The optical bands consist of u', g', r', i', and z' from the DECam and cover ∼17.8 deg² of the SHELA footprint (I. Wold et al. 2018, in preparation). The mid-IR bands include the 3.6 and 4.5 μm from the Infrared Array Camera (IRAC) aboard the Spitzer Space Telescope and cover 24 deg² (Papovich et al. 2016). In addition, we include near-IR photometry in J and K_s from the February 2017 version of the VISTA-CFHT Stripe 82 Near-infrared Survey⁹ (VICS82; Geach et al. 2017), which covers ∼85% of the optical imaging footprint and has 5σ depths of J = 21.3 mag and K_s = 20.9 mag. Figure 1 shows the filter transmission curves for the nine photometric bands used, overplotted with model high-z galaxy spectra illustrating the wavelength coverage of our data set.

Zoom In Zoom Out Reset image size

The DECam images were processed by the NOAO Community Pipeline (CP). A detailed description of the CP reduction procedure can be found in the DECam Data Handbook on the NOAO website;¹⁰ however, we provide a brief summary of the procedure here. First, the DECam images were calibrated using calibration exposures from the observing run. The main calibration steps included an electronic bias calibration, saturation masking, bad pixel masking and interpolation, dark count calibration, linearity correction, and flat-field calibration. Next, the images were astrometrically calibrated with Two Micron All Sky Survey reference images. Finally, the images were remapped to a grid where each pixel is a square with a side length of 027. Observations taken on the same night were then co-added.

The CP data products for the SHELA field were downloaded from the NOAO Science Archive.¹¹ The data products include the co-added images, remapped images, data quality maps (DQMs), exposure time maps (ETMs), and weight maps (WMs). The co-added images from the CP were not intended for scientific use, so we opted to co-add the remapped images. We followed the co-adding procedure of I. Wold et al. (2018, in preparation), which we summarize here. Using the software package SWARP (Bertin et al. 2002), the subimages stored in the FITS files of the remapped images were stitched together and background subtracted. The remapped images were combined using a weighted mean procedure optimized for point sources. The weighting of each image is a function of the seeing, transparency, and sky brightness and is defined by Equation (A3) in Gawiser et al. (2006) as

$$\begin{eqnarray}&&{w}_{i}^{\mathrm{PS}}={\left(\displaystyle \frac{{\mathrm{factor}}_{i}}{{\mathrm{scale}}_{i}\times {\mathrm{rms}}_{i}}\right)}^{2},\end{eqnarray} \tag{ 1 }$$

where scale_i is the image transparency (defined as the median brightness of the bright unsaturated stars after normalizing the brightness measurement of each star by its median brightness across all exposures), rms_i is the root mean square of the fluctuations in background pixels, and factor_i is defined as

$$\begin{eqnarray}&&{\mathrm{factor}}_{i}=1-\exp \left(-1.3\displaystyle \frac{{\mathrm{FWHM}}_{\mathrm{stack}}^{2}}{{\mathrm{FWHM}}_{i}^{2}}\right),\end{eqnarray} \tag{ 2 }$$

where FWHM_stack is the median full width at half maximum (FWMH) of bright unsaturated stars in an unweighted stacked image and FWHM_i is the median FWHM of bright unsaturated stars in each individual exposure.

The seeing and transparency measurements were determined using a preliminary source catalog generated for each resampled image using the Source Extractor (SExtractor) software package (Bertin & Arnouts 1996). The seeing measurements from the final stacked images are listed in Table 1.

Table 1.  DECam Imaging Summary—Seeing and Limiting Magnitude

			u'		g'		r'		i'		z'
Subfield	R.A.	Decl.	FWHM	Depth	FWHM	Depth	FWHM	Depth	FWHM	Depth	FWHM	Depth
ID No.	(J2000)	(J2000)	(arcsec)	(mag)	(arcsec)	(mag)	(arcsec)	(mag)	(arcsec)	(mag)	(arcsec)	(mag)
SHELA-1	1h00m528	−0°00'36''	1.12	25.2	1.06	24.8	1.0	24.8	0.87	24.5	0.86	23.8
SHELA-2	${1}^{{\rm{h}}}{07}^{{\rm{m}}}{02.4}^{{\rm{s}}}$	−0°00'36''	1.2	25.2	1.26	24.8	1.28	24.6	1.39	23.9	0.96	23.6
SHELA-3	1h13m120	−0°00'36''	1.22	25.4	1.36	25.0	1.14	24.7	1.04	24.5	1.21	23.5
SHELA-4	1h19m216	−0°00'36''	1.15	25.3	1.4	24.4	1.05	24.3	1.0	22.1	1.13	23.7
SHELA-5	1h25m312	−0°00'36''	1.21	25.1	1.07	24.9	1.02	24.3	0.93	23.9	0.85	23.6
SHELA-6	1h31m408	−0°00'36''	1.26	25.4	1.37	24.9	1.26	24.5	1.27	24.2	0.86	23.5

Note. The FWHM values are for the stacked DECam images before PSF matching and have units of arcseconds. The magnitudes quoted are the 5σ limits measured in 189-diameter apertures on the PSF-matched images (see Section 2.3).

Download table as:  ASCIITypeset image

After discovering that the original WMs from the CP had values inconsistent with the ETM, we created custom rms maps for the co-added images. The initial rms per pixel was defined as the inverse of the square root of the exposure time. The median of the rms map is scaled to the global pixel-to-pixel rms, which is defined as the standard deviation of the fluxes in good-quality, blank-sky pixels. Good-quality, blank-sky pixels are pixels not included in a source according to our initial SExtractor catalog (see Section 2.3 for discussion of our source extraction procedure) and have an exposure time greater than 0.9 times the median value.

The DECam imaging data were photometrically calibrated with photometry from the SDSS Data Release 11 (DR11; Eisenstein et al. 2011) using only F0 stars. F0 stars were used because their SEDs span all five optical filters while appearing in the sky at a sufficiently high surface density to provide statistically significant numbers in each DECam image. We began by creating a preliminary source catalog for the stacked DECam images using SExtractor and position-matching to the SDSS source catalog. Then, we selected F0 stars using SDSS colors by integrating an F0-star model spectrum from the 1993 Kurucz Stellar Atmospheres Atlas (Kurucz 1979) with each of the five optical SDSS filter curves. For sources in the catalog to be identified as an F0 star, the total color differences, using colors for all adjacent bands, added in quadrature must have been less than 0.35. We then calculate the expected magnitude offset between SDSS and DECam filters for F0 stars, which are as follows: 0.33 for u', 0.02 for g', −0.001 for r', −0.02 for i', and −0.01 for z'. The zero-point for each filter was then calculated as

$$\begin{eqnarray}&&\mathrm{ZPT}=\mathrm{median}({m}_{\mathrm{SDSS}}^{\mathrm{AB}}-{m}_{\mathrm{DECam}}-{\rm{\Delta }}{m}_{\mathrm{offset}}),\end{eqnarray} \tag{ 3 }$$

where Δm_offset is the expected magnitude offset between SDSS and DECam filters for F0 stars. After the zero-points were applied to each stacked image, the image pixel values were converted to units of nJy.

Studying galaxy properties relies on accurate measurements of galaxy colors. One way to obtain accurate colors is to perform fixed-aperture photometry where you measure source fluxes in every band using the same-sized aperture. However, since the DECam images where taken in varying seeing conditions, they have point-spread functions (PSFs) with a range of FWHM. To perform fixed-aperture photometry on these images, the PSFs of all the imaging covering a single patch of sky must be adjusted to have a similar PSF size. We divided the DECam imaging into six subfields (each defined as roughly one DECam pointing). For each subfield, we enlarged the PSFs of the stacked images to match the PSF of the stacked image with the largest PSF in that subfield. For example, in subfield SHELA-1 we matched the PSFs of the g', r', i', and z' stacked images to the PSF of the u' band. To enlarge the PSFs, we adopted the procedure of Finkelstein et al. (2015), who used the IDL deconv_tool Lucy−Richardson deconvolution routine. This routine takes as inputs two PSFs (the desired larger PSF and the starting smaller PSF) and the number of iterations to run and outputs a convolution kernel. The input image PSFs were produced by median-combining the 100 brightest stars (sources with stellar classifications in SDSS DR11) in each image. Before combining, the stars were oversampled by a factor of 10, recentered, and then binned by 10 to ensure that the star centroids aligned. We ran the deconvolution routine with an increasing number of iterations until the PSF of the convolved image (again measured from stacking stars) had a flux within a 7 pixel (189) diameter aperture matched to that of the PSF of the target image to within 5%. The total fluxes were measured in 30 pixel (81) diameter circular apertures. In Figure 2 we show the results of PSF matching in each subfield by displaying a comparison of the enlarged PSFs to the largest PSF as a percent difference.

Zoom In Zoom Out Reset image size

Due to variations in intrinsic galaxy colors and variable image depth and sky coverage, some galaxies will not appear in all bands. To get photometric measurements of all sources in every DECam image, we combined the information in the five optical band images into a single detection image. We followed the procedure of Szalay et al. (1999) and summed the square of the signal-to-noise ratio (S/N) in each band pixel by pixel as follows:

$$\begin{eqnarray}&&{D}_{i}=\sqrt{\sum \displaystyle \frac{{F}_{\mathrm{band},i}^{2}}{{\sigma }_{\mathrm{band},i}^{2}}},\end{eqnarray} \tag{ 4 }$$

where D_i is the detection image ith pixel value, F_band,i is the ith pixel flux in the band image, and σ_band,i is the rms at that ith pixel pulled from the band rms image. A weight image associated with this detection image was created where pixels associated with detection map pixels with data in at least one band have a value of unity and pixels associated with detection map pixels without data have a value of zero.

Photometry was measured on the PSF-matched images using the SExtractor software (v2.19.5; Bertin & Arnouts 1996). Catalogs were created for each of the six SHELA subfields with SExtractor in two-image mode using the detection image described above and cycling through the five DECam bands as the measurement image. In the final source catalog, we maximized the detection of faint sources while minimizing false detections by optimizing the combination of the SExtractor parameters DETECT_THRESH and DETECT_MINAREA. We did this by running SExtractor with an array of combinations of DETECT_THRESH and DETECT_MINAREA and chose the combination of 1.6 and 3, respectively, which detected all sources that appeared real by visual inspection and included the fewest false-positive detections from random noise fluctuations.

We measure source colors in 189-diameter circular apertures (which corresponds to an enclosed flux fraction of 59%–75% for unresolved sources in our subfields). To obtain the total flux, we derived an aperture correction defined as the flux in a 189-diameter aperture divided by the flux in a Kron aperture (i.e., FLUX_AUTO), using the default Kron aperture parameters of PHOT_AUTOPARAM = 2.5, 3.5, which has been shown to calculate the total flux to within ∼5% (Finkelstein et al. 2015). This correction was derived in the r' band on a per-object basis to account for different source sizes and ellipticities and was applied to the fluxes in the other DECam bands per subfield. In areas where the subfields overlapped, sources with positions that matched to within 12 in neighboring subfield catalogs had their fluxes mean-combined after being weighted by the inverse square of their uncertainties (see Section 2.4). The DECam source fluxes were corrected for Galactic extinction using the color excess E(B − V) measurements by Schlafly & Finkbeiner (2011). We obtained E(B − V) values using the Galactic Dust Reddening and Extinction application on the NASA/IPAC Infrared Science Archive (IRAS) website.¹² We queried IRAS for E(B − V) values (the mean value within a 5' radius) for a grid of points across the SHELA field with 4' spacing and assigned each source the E(B − V) value from the closest grid point. The Cardelli et al. (1989) Milky Way reddening curve parameterized by R_V = 3.1 was used to derive the corrections at each band's central wavelength. We compared the extinction-corrected DECam photometry with SDSS DR14 per subfield and band and found agreement for point sources to better than δm < 0.05−0.2 mag in terms of scatter.

We estimated photometric uncertainties in the DECam images by estimating the image noise in apertures as a function of pixels per aperture, N, following the procedure described in Section 2 of Papovich et al. (2016). There are two limiting cases for the uncertainty in apertures with N pixels, σ_N. If pixel errors are completely uncorrelated, the aperture uncertainty scales as the square root of the number of pixels, ${\sigma }_{N}={\sigma }_{1}\times \sqrt{N}$, where σ₁ is the standard deviation of sky background pixels. If pixel errors are completely correlated, then σ_N = σ₁ × N (Quadri et al. 2006). Thus, the aperture uncertainty will scale as N^β with 0.5 < β < 1.

To estimate the aperture noise as a function of N pixels, we measured the sky counts in 80,000 randomly placed apertures ranging in diameter from 027 to 81 across each stacked DECam image. We required apertures to fall in regions of the background sky, which we define as the region where the ETM has the value of at least the median exposure time (ensuring that >50% of each image was considered), excluding detected sources and pixels flagged in the DQM. We also required that the apertures do not overlap with each other. We then estimated σ_N for each aperture with N pixels by computing the standard deviation of the distribution of aperture fluxes from the normalized median absolute deviation, σ_nmad (Beers et al. 1990). We calculated σ₁ by computing σ_nmad for all pixels in the background sky as defined above. Figure 3 shows an example of the measured flux uncertainty in a given aperture, σ_N, as a function of the square root of the number of pixels in each aperture, N, for the five DECam bands in the subfield SHELA-1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Following Papovich et al. (2016), we fit a parameterized function to the noise in an aperture of N pixels, σ_N, as

$$\begin{eqnarray}&&{\sigma }_{N}={\sigma }_{1}(\alpha {N}^{\beta }+\gamma {N}^{\delta }),\end{eqnarray} \tag{ 5 }$$

where σ₁ is the pixel-to-pixel standard deviation in the sky background and α, β, γ, and δ are free parameters. The best-fitting parameters in Equation (5) for the combined DECam images are listed in Table 2. While the second term was intended to aid in fitting the data at large N values, in actuality the second term contributed significantly at all N values, resulting in β ≈ 0.9–1. Nevertheless, our functional fits reproduce the data well as can be seen, for example, in Figure 3. To estimate how correlated the pixel-to-pixel noise is, we fit σ_N with only the first term in Equation (5) and found typical values of β ≈ 0.65–0.70, suggesting slightly correlated pixel-to-pixel noise.

The photometric errors estimated by Equation (5) were scaled to apply to flux measurements outside the region with the median exposure time. The flux uncertainty for the ith source in band b in the subfield f is calculated as

$$\begin{eqnarray}&&{\sigma }_{i,f,b}^{2}={\sigma }_{N,f,b}^{2}{\left(\displaystyle \frac{{\mathrm{rms}}_{i,f,b}}{{\mathrm{rms}}_{\mathrm{med},f,b}}\right)}^{2},\end{eqnarray} \tag{ 6 }$$

where σ_N,f,b is given by Equation (5) for each band and subfield, rms_i,f,b is the value of the rms map at the central pixel location of the ith source in each band and subfield, and rms_med,f,b is the median value of the rms map in each band and subfield. The photometric error estimates exclude Poisson photon errors, which we estimate to contribute <5% uncertainty to the optical fluxes of our high-z candidates.

As described in Section 2.3, all source fluxes are measured in circular apertures of 189 diameter and scaled to total on a per-object basis. Likewise, the flux uncertainties in the apertures were scaled by the same amount to determine the total flux uncertainties, so that the S/N for the total flux is the same as for the aperture flux.

As discussed below, we wish to enhance the validity of our z = 4 galaxy sample by including IRAC photometry in our galaxy sample selection. While we could allow this by position-matching the published Spitzer/IRAC catalog from Papovich et al. (2016) to our DECam catalog, this is not optimal for two reasons. First, the Papovich et al. (2016) catalog is IRAC detected, and so it only includes sources with significant IRAC flux, while for our purposes even a nondetection in IRAC can be useful for calculating a photometric redshift. Second, this catalog uses apertures defined on the positions and shapes of the IRAC sources, while the larger PSF of the IRAC data results in significant blending, which is a larger issue at fainter magnitudes, where we expect to find the bulk of our sources of interest. For these reasons, we applied the Tractor image modeling code (Lang et al. 2016a, 2016b) to perform "forced photometry," which employs prior measurements of source positions and surface brightness profiles from a high-resolution band to model and fit the fluxes of the source in the remaining bands, splitting the flux in overlapping objects into their respective sources. We specifically used the Tractor to optimize the likelihood for the photometric properties of DECam sources in each of IRAC 3.6 and 4.5 μm bands given initial information on the source and image parameters. The input image parameters of IRAC 3.6 and 4.5 μm images included a noise mode, a PSF model, image astrometric information (WCS), and calibration information (the "sky noise" or rms of the image background). The input source parameters included the DECam source positions, brightness, and surface brightness profile shapes. The Tractor proceeds by rendering a model of a galaxy or a point source convolved with the image PSF model at each IRAC band and then performs a linear least-squares fit for source fluxes such that the sum of source fluxes is closest to the actual image pixels, with respect to the noise model. We describe how we use the Tractor to perform forced photometry on IRAC images in detail below.

We begin our source modeling procedure by selecting the fiducial band high-resolution model of each source. We use the fluxes and surface brightness profile shape parameters measured in our DECam detection image because the image combines the information of all sources in the five optical band images (as described in Section 2.3). Second, we use a one-component circular Gaussian to model the PSF. During the modeling of each source, we allow Tractor to optimize the Gaussian σ value, in addition to optimizing a source flux. We find that the median of the optimized Gaussian σ is 080 (equivalent to an FWHM of 188) for both IRAC 3.6 and 4.5 μm images, consistent with those measured from an empirical point response function for the 3.6 and 4.5 μm image (FWHM of 197; see Papovich 2016, Section 3.4).

In practice, we extracted an IRAC image cutout of each source in the input DECam catalog. We selected the cutout size of 16'' × 16''. This cutout size represents a trade-off between minimizing computational costs related to larger cutout sizes and ensuring that the sources lie well within the cutout extent. The sources of interest within the cutout are modeled as either unresolved (i.e., a point source) or resolved based on the DECam detection image. We considered a source to be resolved if an estimated radius r > 01. We define the radius as ${r}_{\mathrm{source}}=a\times \sqrt{b/a}$, where a is a semimajor axis and b/a is an axis ratio. We perform the photometry for resolved sources using a de Vaucouleurs profile (equivalent to a Sérsic profile with n = 4) with shape parameters (semimajor axis, position angle, and axis ratio) measured using our DECam detection image. We have also performed the photometry using an exponential profile (equivalent to a Sérsic profile with n = 1), but we do not find any significant difference between the IRAC flux measurements for the two galaxy profiles. Therefore, we adopt a de Vaucouleurs profile to model all resolved sources. The Tractor simultaneously modeled and optimized the sources of interest and neighboring sources within the cutout. Finally, the Tractor provided the measurement IRAC flux of each DECam source with the lowest reduced chi-squared value. We validated the Tractor-based IRAC fluxes by comparing the fluxes of isolated sources (no neighbors within 3'') to the published Spitzer/IRAC catalog from Papovich et al. (2016). For both bands, we found good agreement with a bias offset of δm < 0.05 mag and a scatter of <0.11 mag down to m = 20.5 mag and a bias offset of δm < 0.13 mag and a scatter of <0.26 down to m = 22 mag.

We selected our sample of high-redshift galaxies using a selection procedure that relies on photometric redshift (z_phot) fitting, which leverages the combined information in all photometric bands used. We obtained z_phot values and z_phot probability distribution functions (PDFs) from the EAZY software package (Brammer et al. 2008). For this analysis, we use the "z_a" redshift column from EAZY, which is produced by minimizing the χ² in the all-template linear combination mode. We also tried the "z_peak" column, and our resulting luminosity function did not change significantly. EAZY assumes the intergalactic medium (IGM) prescription of Madau (1995). We did not use any magnitude priors based on galaxy luminosity functions when running EAZY, as the existing uncertainties in the bright end would bias our results (effectively assuming a flat prior). We then applied a number of selection criteria using the z_phot PDFs from EAZY following the procedure of Finkelstein et al. (2015), which we summarize here, to construct a z ≈ 4 galaxy sample. First, we required a source to have an S/N of greater than 3.5 in the two photometric bands (r' and i' bands) that probe the UV continuum of our galaxy sample, which has been shown to limit contamination by noise to negligible amounts by Finkelstein et al. (2015). Next, we required the integral of each source's z_phot PDF for z > 2.5 to be greater than 0.8 and the integral of the z_phot PDF from z = 3.5–4.5 to be greater than the integral of the z_phot PDF in all other Δz = 1 bins centered on integer-valued redshifts. Finally, we required a source to have no photometric flags on its u', g', r', and i' flux measurements in our photometric catalog. These four bands are crucial for probing both sides of the Lyman break of galaxies in the redshift range of interest. Photometric flags indicate saturated pixels, transient sources, or bad pixels as defined by the CP (see Section 2.2).

Initially, we required sources to pass this z_phot selection procedure when only the DECam and IRAC bands were used. After inspecting candidates from this selection and finding low-z contaminants (see Section 3.2), we moved to include the VISTA J and K photometry when available. The candidate sample of high-z galaxies after including the VISTA J and K photometry significantly reduced the amount of obvious low-z contamination as described in the following section; however, a new class of contaminant appeared in the candidate sample owing to erroneous VISTA photometry for some sources. To overcome this, we required a source to pass the z_phot selection process with and without including the VISTA J and K photometry. This double z_phot selection process resulted in an initial sample of 4364 high-z objects. Next, we performed an investigation into possible contamination, which resulted in additional selection criteria and a refined sample.

Photometric studies of high-z objects can be contaminated by galactic stars and low-z galaxies whose 4000 Å break can mimic the Lyα break of high-z galaxies. The inclusion of IRAC photometry is crucially important for removing galactic stars from our sample, as galactic stars have optical colors very similar to z = 4 objects (Figure 4). While inspecting the photometry and best-fitting SED templates to confirm that our fits were robust and our high-z galaxy candidates were convincing, we found evidence of contamination in our preliminary photometrically selected sample. We explored ways of identifying and removing the contamination, including cross-matching our sample with proper-motion catalogs, SDSS spectroscopy, and X-ray catalogs. Additionally, we implemented machine learning methods.

Zoom In Zoom Out Reset image size

3.2.1. Cross-matching with NOMAD and SDSS

3.2.2. Cross-matching with X-Ray Catalog Stripe 82X

3.2.3. Insights from Machine Learning

We estimated the contamination in our high-z galaxy sample by simulating faint and low-S/N low-z interloper galaxies following a procedure used by Finkelstein et al. (2015). We did this by selecting a sample of bright low-z sources from our catalog, dimming and perturbing their fluxes, and assigning the appropriate uncertainties to their dimmed fluxes. This empirical test assumes that bright low-z galaxies have the same colors as the faint lower-z galaxies. The sources we dimmed all had m_r' = 14.3–18 mag, brighter than our brightest z = 3.5–4.5 candidate source, a best-fitting z_phot of 0.1 < z_phot <  0.6, and no photometric flags in any optical band (see Section 2.2 for definition of photometric flag). We reduced the source fluxes randomly to create a flat distribution of dimmed r'-band magnitude spanning the range of our candidate high-z galaxies (m_r' = 18–26). We assigned flux uncertainties to the dimmed fluxes by randomly drawing from the flux uncertainties of our candidate high-z galaxy sample. We then perturbed the dimmed fluxes, simulating photometric scatter by drawing flux perturbations from a Gaussian distribution with a standard deviation equal to the assigned flux uncertainties. Since the VISTA J and K photometry from the VICS82 Survey covered ∼85% of the total SHELA survey area, ∼15% of our high-z galaxy candidates do not have flux measurements in the VISTA bands. We incorporated this property of our catalog into the mock catalog by randomly omitting dimmed J and K fluxes at the same rate as the fraction of missing J and K fluxes in our final sample for a given m_r' bin.

We created 4.8 million artificially dimmed sources and ran them through EAZY and our selection criteria in an identical manner to our real catalog. We found that 30,594 artificially dimmed sources satisfied our z_phot = 4 selection criteria. We inspected the undimmed and unperturbed properties of these sources and found a range of r'− i' colors (0 < r' − i' < 1.4), and the objects with the reddest r' − i' colors contaminated at the highest rate. We defined the contamination rate in the ith r' − i' color bin as

$$\begin{eqnarray}&&{R}_{i}=\displaystyle \frac{{N}_{\mathrm{dimmed},\mathrm{selected},i}}{{N}_{\mathrm{dimmed},i}},\end{eqnarray} \tag{ 7 }$$

where N_{dimmed,selected,i} is the number of dimmed sources satisfying our z_phot = 4 criteria in the ith r' − i' color bin and N_dimmed,i is the total number of dimmed sources in the ith r' − i' color bin. We created seven r' − i' color bins spanning the range of r'− i' colors recovered (0 < r' − i' < 1.4), with each bin having width = 0.2 mag. The contamination fraction, F, in our high-z galaxy sample was defined as

$$\begin{eqnarray}&&F=\displaystyle \frac{{\displaystyle \sum }_{i}^{}\tfrac{{R}_{i}\times {N}_{\mathrm{total},i}}{(1-{R}_{i})}}{{N}_{z}},\end{eqnarray} \tag{ 8 }$$

where R_i is the contamination rate defined by Equation (7), N_{total, i} is the total number of sources in our object catalog in the ith r' − i' color bin with 0.1 < z_phot < 0.6, and N_z is the number of high-z candidates in our final sample in a given redshift bin. We calculated F as a function of m_r' and m_i' and found contamination fractions of between 0% and 20% generally increasing as a function of magnitude and not exceeding 25% in any magnitude bin above our 50%-completeness limit (m_i' > 23.5). We learned that we can improve our fidelity by implementing an r' − i' < 1.0 color cut, given that the objects with the reddest r' − i' colors contaminated at the highest rate. We then reran our contamination simulation with our final set of selection criteria and found improvement by a few percentage points in two of our brighter magnitude bins (m_i' = 18.75, 20.5) and, again, contamination fractions of between 0% and 20% generally increasing as a function of magnitude and not exceeding 25% in any magnitude bin brighter than our 50%-completeness limit as shown in Figure 5. The addition of the r' − i' < 1.0 color cut reduced our high-z candidate sample size from 3772 to the final size of 3740 with a median z_phot of 3.8. The measured i'-band magnitude and the redshift distribution of our final sample are shown in Figure 6. A summary of our sample selection criteria is listed in Table 3. Bright (m_i' < 22) candidates are plotted in Figure 4, showing the distinct color parameter space they occupy (along with SDSS spectroscopically classified AGNs) compared to the parameter space occupied by SDSS spectroscopically classified stars in the SHELA field.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Summary of z = 4 Sample Selection Criteria

Criterion	Section Reference
S/Nr' > 3.5	3.1
S/Ni' > 3.5	3.1
${\int }_{2.5}^{\infty }\mathrm{PDF}(z)\ {dz}\gt 0.8$	3.1
${\int }_{3.5}^{4.5}\mathrm{PDF}(z)\ {dz}\gt$ All other Δz = 1 bins	3.1
i' − [3.6] > −0.2	3.2.3
r'− i' < 1.0	3.3

Download table as:  ASCIITypeset image

A crucial component of calculating the UV luminosity function is measuring the effective volume of the survey. The effective volume measurement depends on quantifying the survey incompleteness due to image depth and selection effects. To quantify the survey incompleteness, we simulated a diverse population of high-z galaxies with assigned photometric properties and uncertainties consistent with our source catalog and measured the fraction of simulated sources that satisfied our high-z galaxy selection criteria as a function of the absolute magnitude and redshift.

The simulated mock galaxies were given properties drawn from distributions in redshift and dust attenuation (e.g., E[B − V]), while the ages and metallicities were fixed at 0.2 Gyr and solar (Z = Z_⊙), respectively. Because the fraction of recovered galaxies per redshift bin is independent of the number of simulated galaxies per redshift bin as long as low-number statistics are avoided, the redshift distribution was defined to be flat from 2 < z < 6, and the E(B − V) distribution was defined to be lognormal spanning 0< E(B − V) <1 and peaking at 0.2. Mock SEDs were then generated for each galaxy using pythonFSPS¹³ (Foreman-Mackey et al. 2014), a python package that calls the Flexible Stellar Population Synthesis Fortran library (Conroy et al. 2009; Conroy & Gunn 2010). We then integrated each galaxy SED through our nine filters (DECam u', g', r', i', and z'; VISTA J and K; and IRAC 3.6 and 4.5 μm). Each set of mock photometry was then scaled to have an r-band apparent magnitude within a lognormal distribution spanning 18 < m_r' < 27. This distribution ensured that we were simulating the most galaxies at the fainter magnitudes, where we expected to be incomplete, and fewer at bright magnitudes, where we expected to be very complete. Flux uncertainties and missing VISTA fluxes were assigned in the same way as during the contamination test in Section 3.3. Mock galaxies were not assigned photometric flags. This likely results in an incompleteness of only a couple percent, as the fraction of all sources affected by flags is small. All sources with an S/N of greater than 3.5 in any one band are flagged in another band less than 2% of the time on average (never exceeding 4%). In addition, all images have ∼0.5% of all pixels flagged on average (never exceeding 2%).

The photometric catalog of 600,000 mock galaxies was then run through EAZY to generate z_phot PDFs, and our high-z galaxy selection was applied. The completeness was defined as the number of mock galaxies recovered divided by the number of input mock galaxies, as a function of input absolute magnitude and redshift. Figure 7 shows the results of our simulation. We define the 50%-completeness limit as the absolute magnitude where the area under the curve falls to less than 50% of the areas under the average of the M_UV,i' = −25 to M_UV,i' = −28 curves, which we find to be at M_UV,i' = −22 (m = 24 for z = 4).

We utilize the effective volume method to correct for incompleteness in deriving our luminosity function. The effective volume (V_eff) can be estimated as

$$\begin{eqnarray}&&{V}_{\mathrm{eff}}({M}_{i^{\prime} })=\int \displaystyle \frac{{{dV}}_{C}}{{dz}}C({M}_{\mathrm{UV},i^{\prime} },z){dz},\end{eqnarray} \tag{ 9 }$$

where $\tfrac{{{dV}}_{C}}{{dz}}$ is the comoving volume element, which depends on the adopted cosmology, and C(M_UV,i', z) is the completeness as calculated in Section 3.4. The integral was evaluated over z = 3–5.

To calculate the luminosity function, we convert the apparent i-band AB magnitudes (m_i') to the absolute magnitude at rest-frame 1500 Å (M_UV,i') using the following formula:

$$\begin{eqnarray}&&{M}_{\mathrm{UV},i^{\prime} }={m}_{i}-5\mathrm{log}({d}_{L}/10\mathrm{pc})+2.5\mathrm{log}(1+z),\end{eqnarray} \tag{ 10 }$$

where d_L is the luminosity distance in pc. The second and third terms on the right-hand side are the distance modulus.

Our UV luminosity function is shown as red diamonds in Figure 8. We note that we do not include the luminosity function data points in bins M_UV,i' ≥ −22 in our analysis, as these bins are below our 50%-completeness limit as discussed in Section 3.4, where the completeness corrections are unreliable owing to the low S/N of our data. This is confirmed by comparing to the HST CANDELS results in these same magnitude bins, which are more reliable owing to their higher S/N. In Figure 8 we also include the UV luminosity function of z = 4 UV-selected galaxies from deeper HST imaging by Finkelstein et al. (2015) as green squares. In the bin where we overlap with this data set (M_UV,i' = −22.5), both results are consistent; however, if the luminosity function derived from HST imaging is extrapolated to brighter magnitudes, it would fall off more steeply than our luminosity function. While our luminosity function declines from M_UV,i' = −22 to M_UV,i' = −24, it flattens out at brighter magnitudes before turning over again.

Zoom In Zoom Out Reset image size

As our sample includes all galaxies that exhibit a Lyman break, we expect that this flattening is due to the increasing importance of AGNs at these magnitudes. The large volume surveyed by SDSS data has led to the selection and spectroscopic follow-up of AGNs at many redshifts. SDSS AGN studies have found the AGN UV luminosity function to exhibit a double power-law (DPL) shape (e.g., Richards et al. 2006). In Figure 8 we show as red crosses the z = 4 AGN UV number densities derived from the SDSS DR7 catalog (Schneider et al. 2010) by Akiyama et al. (2018), who select AGNs to M_UV,i' > −28.9. We can see that at the magnitudes where we overlap, −27 ≲ M_UV,i' ≲ −26, the agreement with our data is excellent. The only difference is at M = −28, where our survey detects two quasars while the AGN luminosity function by Akiyama et al. (2018) would predict less than one in our volume. We attribute this difference to cosmic variance. By combining our data with the star-forming galaxy number densities from HST imaging and the SDSS AGN number densities, our data can potentially provide a robust measurement of the bright-end slope of the star-forming galaxy luminosity function and the faint-end slope of the AGN luminosity function, which we explore in the following section.

With our luminosity function in agreement with the faint end of the AGN luminosity function from SDSS DR7 and the bright end of the star-forming galaxy luminosity function from CANDELS, we attempt to simultaneously fit empirically motivated functions to both components. For the AGN component we use a DPL function motivated by the AGN UV luminosity function work on large, homogeneous quasar samples (e.g., Boyle et al. 2000; Richards et al. 2006; Hopkins et al. 2007; Croom et al. 2009, and references therein). The function form of a DPL follows

$$\begin{eqnarray}&&{\rm{\Phi }}(M)=\displaystyle \frac{{{\rm{\Phi }}}^{* }}{{10}^{0.4(\alpha +1)(M-{M}^{* })}+{10}^{0.4(\beta +1)(M-{M}^{* })}},\end{eqnarray} \tag{ 11 }$$

where Φ* is the overall normalization, M* is the characteristic magnitude, α is the faint-end slope, and β is the bright-end slope.

For the star-forming galaxy UV luminosity function we consider separately both a Schechter function and a DPL, as well as including magnification via gravitational lensing with both functions. The Schechter (1976) function has been found to fit the star-forming galaxy UV luminosity function well across all redshifts (e.g., Steidel et al. 1999; Bouwens et al. 2007; Finkelstein et al. 2015). The Schechter function is described as

$$\begin{eqnarray}&&{\rm{\Phi }}(M)=\displaystyle \frac{0.4\ \mathrm{ln}(10)\ {{\rm{\Phi }}}^{* }}{{10}^{0.4(\alpha +1)(M-{M}^{* })}{e}^{{10}^{-0.4(M-{M}^{* })}}},\end{eqnarray} \tag{ 12 }$$

where Φ* is the overall normalization, M* is the characteristic magnitude, and α is the faint-end slope.

We consider the effects of gravitational lensing on the shape of the star-forming galaxy UV luminosity function. Gravitational lensing can distort the shape and magnification of distant sources as the paths of photons from the source get slightly perturbed into the line of sight of the observer. Lima et al. (2010) showed that this magnification can contribute to a bright-end excess where the slope of the intrinsic luminosity function is sufficiently steep. A magnification distribution for a given source redshift must be estimated by tracing rays through a series of lens planes derived from simulations such as the Millennium Simulation as done by Hilbert et al. (2007). Van der Burg et al. (2010) showed that a Schechter function corrected for magnification can fit the bright end of the luminosity function at z = 3 better than a Schechter fit alone. They inspect the sources that make up the excess and find nearby massive foreground galaxies or groups of galaxies that could act as lenses. We incorporate the effects of gravitational lensing in our fitting by creating a lensed Schechter function parameterization following the method of Ono et al. (2018), who adapt the method of Wyithe et al. (2011). We also produce a lensed DPL function. After performing our simultaneous fitting method, which we describe in the following paragraph, we found there to be no difference in the best-fitting parameters of the fits including and excluding the effects of lensing. This is consistent with the work of Ono et al. (2018), who found that taking into account the effects of lensing improves the galaxy luminosity function fit at z > 4 and not at z = 4, where a DPL fit is preferred. Therefore, we do not consider the lensed parameterizations further.

We employ a Markov chain Monte Carlo (MCMC) method to define the posteriors on our luminosity function parameterizations. We do this using an IDL implementation of the affine-invariant sampler (Goodman & Weare 2010) to sample the posterior, which is similar in production to the emcee package (Foreman-Mackey et al. 2013). Each of the 500 walkers was initialized by choosing a starting position with parameters determined by eye to exhibit a good fit, perturbed according to a normal distribution. We do not assume a prior for any of our free parameters.

We account for Eddington bias in our fitting routine. Rather than directly comparing the observations to a given model, we forward-model the effects of Eddington bias into the luminosity function model and compare this "convolved" model with our observations. We do this by realizing, for each set of luminosity function parameters, a mock sample of galaxies for that given function, where each galaxy has a magnitude according to the given luminosity function distribution. The magnitude of each simulated object is perturbed by an amount drawn from a normal distribution centered on zero with a width equal to the real sample median uncertainty in the corresponding magnitude bin. After perturbing, we then rebin the simulated luminosity distribution, and this binned luminosity function is used to calculate the χ² for that MCMC step. This is repeated for each step.

We ran our MCMC fitting algorithm twice. In both runs, we simultaneously fit a DPL function to the AGN component of the luminosity function while fitting the galaxy component. In the first run, we fit a Schechter function to the galaxy luminosity function, while in the second we fit a second DPL function to the galaxy luminosity function. We burn each chain for 2 × 10⁵ steps, which allows the chain to reach convergence for all free parameters, verified by examining the parameter distributions in independent groups of 10⁴ steps, which cease to evolve much past 10⁵ steps. We then measure the posterior for each parameter from the final 5 × 10⁴ steps. Our fiducial values for each parameter are then derived as the median and 68% central confidence range from the posterior distributions. The best-fitting parameters for the two fits and their corresponding χ² values are listed in Table 4. We overplot the fits to the luminosity function data in Figure 8. To calculate the χ² for our fits, we must compare the observed data with the luminosity function fits after forward-modeling the effects of Eddington bias into our luminosity function fits, just as we did during the fitting process. The absolute residuals of our "convolved" fits and the observed data are shown in the inset in Figure 8.

Table 4.  Fit Parameters for Luminosity Functions

	AGN Component					Galaxy Component
Fit Name	Log Φ*	M*	α	β		Log Φ*	M*	α	β	χ2
		(mag)	[Faint End]	[Bright End]			(mag)	[Faint End]	[Bright End]
DPL+DPL	−7.32${}_{-0.18}^{+0.21}$	−26.5${}_{-0.3}^{+0.4}$	−1.49${}_{-0.21}^{+0.30}$	−3.65${}_{-0.24}^{+0.21}$		−3.12${}_{-0.10}^{+0.09}$	−20.8${}_{-0.15}^{+0.16}$	−1.71 ± 0.08	−3.80 ± 0.10	42
DPL+Sch	−7.48${}_{-0.34}^{+0.58}$	−26.7${}_{-0.4}^{+1.1}$	−2.08${}_{-0.11}^{+0.18}$	−3.66${}_{-0.34}^{+0.68}$		−3.25 ± 0.06	−21.3 ± 0.06	−1.81 ± 0.05	⋯	71

Note. ${{\rm{\Phi }}}^{* }$ in units of Mpc⁻³ mag⁻¹. The parameters for the galaxy component of DPL+Sch correspond to a Schechter function (Equation (12)), and the remaining sets of parameters correspond to a double power-law function (Equation (11)).

Download table as:  ASCIITypeset image

These fitting results show that our data prefer the function that is a sum of two DPL functions, one for the AGN component and one for the galaxy component, henceforth the DPL+DPL fit. The DPL+DPL fit has a χ² = 42 over the entire range considered, −29 < M_UV,i' < −17. The fit that is a sum of a DPL AGN component and a Schechter galaxy component, henceforth the DPL+Schechter fit, has a χ² = 71. We investigate which fit is preferred by the data using the Bayesian information criterion for the two fits (Liddle 2007). The difference in the Bayesian information criterion value for our two fits can be defined as

$$\begin{eqnarray}&&{\rm{\Delta }}{\rm{BIC}}={\chi }_{2}^{2}-{\chi }_{1}^{2}+({k}_{2}-{k}_{1})\ {\rm{ln}}\ N,\end{eqnarray} \tag{ 13 }$$

where ${\chi }_{2}^{2}$ is the ${\chi }^{2}$ for the DPL+Schechter fit, ${\chi }_{2}^{2}$ is the χ² for the DPL+DPL fit, k₂ and k₁ are the number of fit parameters in the DPL+Schechter fit and the DPL+DPL fit, respectively, and N is the number of data points used during fitting. We find a ΔBIC = 26, which suggests that the DPL+DPL fit is very strongly preferred over the DPL+Schechter fit, as ΔBIC exceeds a significance value of 10 (Liddle 2007). Upon comparing the fits' residuals, we see that the data points in the magnitude bins at M_UV,i' = −24 mag and M_UV,i' = −23 mag are driving the preference for the DPL+DPL fit. If we exclude these two bins, the DPL+DPL fit χ² = 41, the DPL+Schechter fit χ² = 58, and the ΔBIC = 14, which indicates that the DPL+DPL fit is still strongly statistically preferred to the DPL+Schechter fit. However, given that the difference between the fits is driven by just a few data points, we do not believe we can firmly rule out a Schechter form for the star-forming galaxy component.

Given our method of fitting the luminosity function with a component for AGNs, we explored the SDSS spectral catalog for spectroscopically confirmed AGNs in SHELA and considered the effectiveness of our selection procedure at recovering AGNs. This is crucial, as our photometric redshift selection process did not include templates dominated by AGN emission, though the strong Lyman break inherent in these sources should still yield an accurate redshift. To confirm this, we queried the spectral catalog from SDSS (DR13; Albareti et al. 2017) using the SDSS CasJobs website¹⁴ and cross-matched the results with our entire DECam catalog. Each match required a separation of less than 04 and the SDSS spectra to be unflagged (i.e., ZWARNING = 0). We found zero spectra with an SDSS classification of "Galaxy" and 53 classified as "AGN" with spectroscopic redshifts (z_spec) greater than 3.2. The distribution of SDSS z_spec for this sample is shown in the left panel of Figure 9 in blue. Of the 32 AGNs with 3.4 < z_spec < 4.6, 23 were selected by our z_phot selection process, with seven of the nine missed AGNs having 3.4 < z_spec < 3.5. The z_spec distribution for the AGN subsample selected by our selection procedure is also shown in the left panel of Figure 9 in red. We used these samples to compute our differential completeness of AGNs and compared it to our simulated completeness in the right panel of Figure 9. We found that our completeness of spectroscopically confirmed AGNs is consistent with our simulated completeness except in the 3.4 < z < 3.6 bin, where we recovered only two of the nine spectroscopically confirmed AGNs when our simulations predicted we should recover seven to eight. This difference could be due to small-number statistics. We investigated why the seven spectroscopically confirmed AGNs in the 3.4 < z < 3.6 bin were not selected by our method and found that these sources had photometry, particularly the u' and g' bands, preferred by galaxies templates at lower redshift in our z_phot-fitting code EAZY. We attribute bright u' and g' fluxes to the larger far-UV continuum levels of AGNs or strong Lyα emission as compared to non-AGN galaxies, which would weaken the Lyα break in these sources. This could imply significant leaking Lyman-continuum radiation from these AGNs, which has implications on reionization (e.g., Smith et al. 2018). The reliability of our selection procedure to recover AGNs across the majority of our redshift range of interest and the fact that our luminosity function is consistent with the AGN luminosity function from SDSS suggest that our incompleteness to AGNs is not substantial.

Zoom In Zoom Out Reset image size

We compare our derived fits to the star-forming galaxy luminosity function with measurements from the literature in Figure 10. At magnitudes fainter than M_UV,i' > −22, the DPL galaxy component of the data-preferred DPL+DPL fit is very similar to the Schechter component of the DPL+Schechter fit. These results are in strong agreement with the luminosity function from the CFHT Deep Legacy Survey by van der Burg et al. (2010) at all magnitudes where they overlap. They are also in strong agreement with the luminosity function from HST legacy survey data by Finkelstein et al. (2015), which is included in our fitting. The luminosity function from the HST legacy survey data by Bouwens et al. (2015) and the luminosity function from Ono et al. (2018) using Subaru HSC data across 100 deg² in the HSC-SSP are generally consistent with the results presented here. The HSC-SSP luminosity function and that from Bouwens et al. (2015) have volume densities ∼2 times larger than the work by van der Burg et al. (2010) and Finkelstein et al. (2015) at magnitudes fainter than M_UV,i' > −21. This factor is larger than the ∼10% uncertainty expected as a result of cosmic variance in the HST fields, which cover  50 times less area than the HSC-SSP (Finkelstein et al. 2015). We do not know the cause of this difference, though we note that the HSC-SSP selection was done with optical imaging only, and we found in our study that without the addition of the IRAC data the contamination rate was significantly higher, although we acknowledge that there are multiple differences in the selection techniques between these studies.

Zoom In Zoom Out Reset image size

At the bright end, the luminosity function from Ono et al. (2018) extends to magnitudes brighter than M_UV,i' < −22.5 and falls between the bright end of our DPL galaxy component and our Schechter galaxy component. Ono et al. (2018) find that their luminosity function is best fit by a DPL with a steeper bright-end slope (β = −4.33) than our galaxy DPL component (β = −3.80). We also include the z = 4 galaxy luminosity function from the 4 deg² ALHAMBRA survey by Viironen et al. (2018), who used a z_phot PDF analysis to create a luminosity function marginalizing over both redshift and magnitude uncertainties. Viironen et al. (2018) find volume densities at the bright end larger than existing luminosity functions. In fact, their luminosity function follows a Schechter function with a normalization offset of ∼0.5 dex above the Schechter form from Finkelstein et al. (2015). The cause of this difference is unclear.

5.1.1. Comparison to Semianalytic Models (SAMs)

Figure 11 shows the predicted z = 4 UV luminosity functions from Yung et al. (2018). These predictions come from a set of SAMs, which contain the same physical processes as the models presented in Somerville (2015) but have been updated and recalibrated to the Planck 2016 cosmological parameters. We note that while these models include black hole accretion and the effects of AGN feedback, we do not examine the contribution to the UV luminosity function from black hole accretion here, focusing instead on the contribution from star formation. Their fiducial model with dust (solid black line) assumes that the molecular gas depletion time is shorter at higher gas densities (as motivated by observations), leading to an effective redshift dependence as high-redshift galaxies are more compact and have denser gas on average. On an SFR surface density versus gas surface density plot, this model would show a steeper dependence of SFR density on gas surface density than the classical Kennicutt–Schmidt relationship (e.g., Kennicutt & Evans 2012), above a critical gas surface density (for details see Somerville 2015). We also show their model with a fixed molecular gas depletion time, similar to that used in Somerville (2015), as seen in local spiral galaxies (e.g., Bigiel et al. 2008; Leroy et al. 2008).

Zoom In Zoom Out Reset image size

The Yung et al. (2018) dust models assume that the V-band dust optical depth is proportional to the cold gas metallicity times the cold gas surface density. The UV attenuation is then computed using a Calzetti attenuation curve and a "slab" model (for details see Somerville 2015). The normalization of the dust optical depth (physically equivalent to the dust-to-metal ratio) is allowed to vary as a function of redshift and was adjusted to fit the bright end of the luminosity from the previously published compilation of UV luminosity observations from Finkelstein (2016). (It was not adjusted to fit the new observations presented here.)

At the faint end, the model predictions are insensitive to the assumed star formation efficiency and mainly reflect the treatment of outflows driven by stars and supernovae. The models have a higher normalization than the observed luminosity function (which at these magnitudes comes from the CANDELS data set), although the faint-end slope is similar. This could be caused by stellar feedback in the models being too weak, as these models were tuned to match the z = 0 stellar mass function (Somerville 2015). This suggests that stellar feedback is stronger/more efficient at z = 4 (i.e., mass-loading rates are higher, or re-infall of ejected gas is slower) than at z = 0 (see also see White et al. 2015). This was also seen by Song et al. (2016) when comparing the z = 4 stellar mass function to a number of similar models—the observed stellar mass function also had a lower normalization at z = 4; as the stellar mass function steepened from z = 4 to 8, this discrepancy weakened, hence their conclusion of a weakening impact of feedback on the faint end with increasing redshift.

At the bright end, the Yung et al. (2018) model with dust is consistent with both of our fits to M_UV,i' > −22.5, lying closer to the Schechter fit at brighter magnitudes. Interestingly, at these bright magnitudes both the Schechter fit and DPL fit rule out at high significance the model with dust and fixed molecular gas depletion time, indicating that star formation must scale nonlinearly with molecular gas surface density (or some related quantity that evolves with redshift). This implies that star formation is more efficient at z = 4 than today. This is, of course, dependent on the dust model in these simulations—if bright galaxies had no dust, then these model predictions (which include dust; models without dust are shown for comparison) may not be accurate. However, there is a variety of evidence that bright/massive UV-selected galaxies at these redshifts do contain non-negligible amounts of dust (e.g., Finkelstein et al. 2012; Bouwens et al. 2014), and there is a not insignificant population of extremely massive dusty star-forming galaxies already in place (e.g., Casey et al. 2014). Finally, Finkelstein et al. (2015) compared a similar set of models to the CANDELS-only UV luminosity functions, finding that models with a diffuse dust component only (e.g., no birth cloud) provided the best match to the data (albeit at fainter magnitudes than considered here).

We compare our derived AGN luminosity function fits to measurements from the literature in Figure 12. At the bright end, our fits are consistent with previous studies where they overlap (−27.5 < M_UV,i' < −25.5). The previous studies include a study by Richards et al. (2006), who used the z = 4 AGN SDSS DR3 sample, and a study by Akiyama et al. (2018), who selected z = 4 AGNs from SDSS DR7. We convert the magnitudes M_i(z=2) at z = 3.75 from Richards et al. (2006) to M_UV,i' at z = 3.8 by adding an offset of 1.486 mag (Richards et al. 2006).

Zoom In Zoom Out Reset image size

At the faint end, we compare our fits to results from studies that derive AGN luminosity functions from spectroscopic observations of candidates selected via color and size criteria (Glikman et al. 2011; Ikeda et al. 2011; Niida et al. 2016). In addition, we include studies that rely on a z_phot selection using deep multiwavelength photometry in the COSMOS field (Masters et al. 2012) and the CANDELS GOODS-S field (Giallongo et al. 2015), where Giallongo et al. (2015) had the additional criteria of requiring an an X-ray detection of F_X > 1.5 × 10⁻¹⁷ erg cm⁻² s⁻¹ in deep 0.5–2 keV Chandra imaging. The average redshift of these samples is z = 4, slightly higher than our sample. For a consistent comparison with other works we scale the luminosity functions of Glikman et al. (2011), Ikeda et al. (2011), Niida et al. (2016), and Masters et al. (2012) up by a factor of 1.3 using the redshift evolution function (1 + z)^−6.9 from Richards et al. (2006). The sample used by Giallongo et al. (2015) had a redshift range of 4 < z < 4.5, so we the scale the luminosity function by a factor of 1.8. Finally, we consider the luminosity function of Akiyama et al. (2018) derived from the HSC Wide Survey optical photometry.

The faint end of the AGN DPL component in our DPL+DPL fit predicts volume densities at −26 < M_UV,i' < −23.5 about 0.3 dex lower than those found by Richards et al. (2006), Ikeda et al. (2011), Masters et al. (2012), Niida et al. (2016), and Akiyama et al. (2018). However, the luminosity function of Akiyama et al. (2018) flattens and falls toward our fit at M_UV,i' ∼ −22. Our faint-end slope ($\alpha =-{1.49}_{-0.21}^{+0.30}$) is in agreement with the values found by these studies. The cause of the 0.3 dex difference is unclear.

In the case of the AGN DPL component in our DPL+Schechter fit, the steeper faint-end slope ($\alpha =-{2.08}_{-0.11}^{+0.18}$) predicts volume densities in agreement with Glikman et al. (2011) and Giallongo et al. (2015), who predict a significantly higher volume density of faint AGNs than the other studies. We note that while these studies have small numbers of AGNs per magnitude bin, which may make the samples more sensitive to cosmic variance and false positives, Glikman et al. (2011) select candidates from a relatively large survey area of 3.76 deg² and observe broad emission lines, indicative of quasars, in every candidate spectrum. Furthermore, Giallongo et al. (2015) require candidates to have a significant X-ray detection in Chandra imaging.

In summary, our observed z = 4 UV luminosity function is best fit by the DPL+DPL fit, but both the DPL+DPL and the DPL+Schechter fits show agreement with existing AGN luminosity function studies at the bright end and around the knee. At the faint end the DPL+DPL fit predicts smaller volume densities of AGNs than other studies, while the DPL+Schechter fit predicts among the largest volume densities at the faintest magnitudes.

Our two fits differ in two ways: (1) the functional form used to fit the star-forming galaxy component is a DPL in the DPL+DPL fit and a Schechter in the DPL+Schechter fit, and (2) the component that accounts for the excess over an extrapolated Schechter at the bright end of existing star-forming galaxy luminosity functions. The DPL+Schechter fit accounts for the excess with a steeper faint end of the DPL AGN component, while the DPL+DPL fit accounts for the majority of the excess with the bright end of the DPL star-forming galaxy component. Thus, the fits predict dramatically different compositions of sources making up the bin (M_UV,i' = −23.5 mag) at the center of the excess. In this bin the DPL+Schechter fit predicts that AGNs outnumber galaxies by a 17:1 ratio (an AGN fraction of ∼94%), while the DPL+DPL fit predicts that galaxies outnumber AGNs by a 4.3:1 ratio (an AGN fraction of ∼18%). A simple experiment to test these predictions would be to use ground-based optical spectroscopy to follow up a fraction of the 298 high-z candidates we find in the M_UV,i' = −23.5 bin and count the fraction of spectra exhibiting broad emission lines and/or highly ionized lines (e.g., N v, He ii, C iv, Ne v), indicating accretion onto supermassive black holes. This experiment would also provide an independent measurement of our contamination fraction, which we estimate via simulations (Section 3.3). The measured fraction of AGNs can provide strong empirical proof in favor of the DPL+DPL fit or the DPL+Schechter fit.

One drawback of fitting the total UV luminosity function is the unknown contribution of composite objects. The total UV luminosity function likely contains a population of composite objects with both AGN activity and star formation contributing to their observed UV flux. This population may cause functions like a DPL or Schechter function, which have been used to fit AGN-only and star-forming-galaxy-only luminosity functions, to be poor fits to the total UV luminosity function. Spectroscopic follow-up as described in this section may aid in elucidating the frequency and impact of composite objects.

Here we calculate the rest-frame UV emissivities (also known as specific luminosity densities) and compare the output from AGNs and galaxies. We calculate the rest-frame UV emissivity of AGNs and galaxies by integrating each luminosity function from $-30\lt {M}_{\mathrm{UV},i^{\prime} }\lt -17$. Results are shown in Table 5. In the case of DPL+DPL fit, where the galaxy component of the UV luminosity function is represented by a DPL function and the AGN component is represented by a DPL, galaxies produce a UV luminosity density at 1500 Å (ρ₁₅₀₀) greater than that of AGNs by a factor of ∼90. In the case of the DPL+Schechter fit, where the galaxy component of the UV luminosity function is instead represented by a DPL function, the galaxy ρ₁₅₀₀ is greater than the AGN ρ₁₅₀₀ by a factor of ∼19. In either case, galaxies are the dominant nonionizing UV-emitting population, though if AGNs do have a steeper faint-end slope (as would need to be the case if galaxies follow a Schechter form), then AGNs are non-negligible.

Table 5.  UV Luminosity Densities

	AGN Component			Galaxy Component
Fit Name	ρ1500	ρ912		ρ1500	ρ912
DPL+DPL	2.0${}_{-0.4}^{+0.7}$	1.2${}_{-0.3}^{+0.4}$		184${}_{-7}^{+6}$	⋯
DPL+Sch	10${}_{-3}^{+4}$	5.8${}_{-1.8}^{+2.6}$		187${}_{-6}^{+6}$	⋯

Note. All values are in units of 10²⁴ ergs s⁻¹ Hz⁻¹ Mpc⁻³.

Download table as:  ASCIITypeset image

We convert the AGN ρ₁₅₀₀ to a UV luminosity density at 912 Å (ρ₉₁₂) and compare our results to other studies by reproducing Figure 1 from Madau & Haardt (2015) in Figure 13. To convert ρ₁₅₀₀ to ρ₉₁₂, we assume an AGN spectral shape of a DPL with a slope of α_ν = −1.41 (Shull et al. 2012) between 1000 and 1500 Å and a slope of α_ν = −0.83 (Stevans et al. 2014) between 912 and 1000 Å. We assume an AGN ionizing escape fraction of unity as is found for most bright AGNs (Giallongo et al. 2015). Results are shown in Table 5. If we directly follow the work of Giallongo et al. (2015) and assume a spectral shape with α_ν = −1.57 (Telfer et al. 2002) between 1200 and 1500 Å and a slope of α_ν = −0.44 (Vanden Berk et al. 2001) between 912 and 1200 Å, we find ρ₉₁₂ values that are 10% smaller. As shown in Figure 13, our two fits predict that values of ρ₉₁₂ straddle existing values from observations. The ρ₉₁₂ predicted by our preferred DPL+DPL fit (upward-pointing red triangle) falls below the line by Hopkins et al. (2007) and is too small to keep the IGM ionized at z ∼ 4. On the other hand, the DPL+Sch fit predicts a ρ₉₁₂ (downward-pointing red triangle) that falls near the points from studies that found numerous faint AGNs at z ∼ 4 (Glikman et al. 2011; Giallongo et al. 2015). This would imply that the AGN population could alone contribute enough hydrogen-ionizing emission to keep the universe ionized at z ∼ 4 and suggests that AGNs may have played a significant role at early times. This scenario can be further constrained by the spectroscopic experiment proposed in the preceding subsection.

Zoom In Zoom Out Reset image size