Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

Bridging Star-Forming Galaxy and AGN Ultraviolet Luminosity Functions at $z=4$ with the SHELA Wide-Field Survey

Bridging Star-Forming Galaxy and AGN Ultraviolet Luminosity Functions at $z=4$ with the SHELA... Draft version June 15, 2018 Typeset using L T X twocolumn style in AASTeX62 Bridging Star-Forming Galaxy and AGN Ultraviolet Luminosity Functions at z = 4 with the SHELA Wide-Field Survey 1 1 1 2, 3, 4 2, 3 Matthew L. Stevans, Steven L. Finkelstein, Isak Wold, Lalitwadee Kawinwanichakij, Casey Papovich, 1 5, 6 1 5, 6 1 Sydney Sherman, Robin Ciardullo, Jonathan Florez, Caryl Gronwall, Shardha Jogee, 7, 8 7 Rachel S. Somerville, and L. Y. Aaron Yung Department of Astronomy, University of Texas at Austin, Austin, TX 78705 Department of Physics and Astronomy, Texas A&M University, College Station, TX, 77843-4242 USA George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M University, College Station, TX, 77843-4242 USA LSSTC Data Science Fellow Department of Astronomy & Astrophysics, The Pennsylvania State University, University Park, PA 16802 Institute for Gravitation and the Cosmos, The Pennsylvania State University, University Park, PA 16802 Department of Physics & Astronomy, Rutgers University, Piscataway, New Jersey 08854 Center for Computational Astrophysics, Flatiron Institute, New York, NY, 10010 (Received March 6, 2018; Revised June 8, 2018; Accepted June 9, 2018) Submitted to ApJ ABSTRACT We present a joint analysis of the rest-frame ultraviolet (UV) luminosity functions of continuum- selected star-forming galaxies and galaxies dominated by active galactic nuclei (AGNs) at z  4. These 3,740 z  4 galaxies are selected from broad-band imaging in nine photometric bands over 18 deg in the Spitzer/HETDEX Exploratory Large Area Survey (SHELA) eld. The large area and moderate depth of our survey provide a unique view of the intersection between the bright end of the galaxy UV luminosity function (M < 22) and the faint end of the AGN UV luminosity function. We AB do not separate AGN-dominated galaxies from star-formation-dominated galaxies, but rather t both luminosity functions simultaneously. These functions are best t with a double power-law (DPL) for both the galaxy and AGN components, where the galaxy bright-end slope has a power-law index of +0:30 3:80  0:10, and the corresponding AGN faint-end slope is = 1:49 . We cannot rule AGN 0:21 out a Schechter-like exponential decline for the galaxy UV luminosity function, and in this scenario +0:18 the AGN luminosity function has a steeper faint-end slope of 2:08 . Comparison of our galaxy 0:11 luminosity function results with a representative cosmological model of galaxy formation suggests that the molecular gas depletion time must be shorter, implying that star formation is more ecient in bright galaxies at z = 4 than at the present day. If the galaxy luminosity function does indeed have a power-law shape at the bright end, the implied ionizing emissivity from AGNs is not inconsistent with previous observations. However, if the underlying galaxy distribution is Schechter, it implies a signi cantly higher ionizing emissivity from AGNs at this epoch. Keywords: galaxies: high-redshift, galaxies: active{quasars, galaxies: luminosity function 1. INTRODUCTION omy. With the number of massive galaxies increasing from z  4 2 (Marchesini et al. 2009; Muzzin et al. Explaining how galaxies grow and evolve over cosmic 2013) and the positive relation between stellar mass and time is one of the main goals of extragalactic astron- star formation rate, by studying the properties of galax- ies with the highest star formation rates at z  4 we Corresponding author: Matthew L. Stevans can glean how the most massive galaxies built up their [email protected] stellar mass. The use of multi-wavelength photometry arXiv:1806.05187v1 [astro-ph.GA] 13 Jun 2018 2 Stevans et al. and the Lyman break technique has revolutionized the point sources. However, this method is less reliable study of galaxies in the z > 2 universe (e.g., Steidel near the photometric limit especially in ground-based et al. 1996). These tools are currently the most ecient imaging with poor seeing. For example, recent work by for selecting large samples of high-redshift star-forming Akiyama et al. (2018) has shown such a morphological galaxies for extensive study. A power tool for under- selection can select a sample of point sources with only standing the distribution of star-formation at high red- 40% completeness and 30% contamination at i = 24 mag shifts is the rest-frame UV luminosity function. This in photometry with median seeing conditions of 0: 7 and probes recent unobscured star-formation directly over 5- depths of i = 26:4 mags. the last 100 Myr and is, therefore, a fundamental tracer The shape of the AGN luminosity function is of in- of galaxy evolution. terest as well, as a steep faint end can result in a non- The shape of the star-forming galaxy UV luminos- negligible contribution of ionizing photons from AGNs ity function at z = 4 has been dicult to pin down to the total ionizing budget. Current uncertainties in the at the bright end. The characteristic luminosity of the literature at z  4 are large (Glikman et al. 2011, Mas- Schechter function, which is often used to describe the ters et al. 2012, Giallongo et al. 2015), thus AGNs have luminosity function in eld environments, ranges over a received renewed interest in the literature with regards few orders of magnitude (e.g., Steidel et al. 1999; Finkel- to reionization (e.g., Madau & Haardt 2015; McGreer stein et al. 2015; Viironen et al. 2017), and there is grow- et al. 2017). ing evidence of an excess of galaxies over the exponen- Studying both AGN-dominated and star-formation- tially declining bright end of the Schechter function (e.g., dominated UV-emitting galaxies simultaneously is pos- van der Burg et al. 2010; Ono et al. 2018). The uncer- sible given a large enough volume. The Hyper Suprime- tainty at the bright end is due in part to cosmic variance Cam (HSC) Subaru strategic program (SSP) has de- and the small area of past surveys which miss the bright- tected large numbers of both types of objects at z = 4 est galaxies with the lowest surface density. The largest using their optical-only data. However, Akiyama et al. z  4 spectroscopically observed sample used in a pub- (2018) opt to use their excellent ground-based resolution lished luminosity function is from the VIMOS VLT Deep (0.6 ; 4.27 kpc at z = 4) to remove extended sources Survey (VVDS, Le F evre et al. 2013) consisting of 129 and focus on the AGN population separately. Ono et al. 2 0 spectra from  1 deg (Cucciati et al. 2012). This small (2018) selects z  4 galaxies (and AGNs) as g -band sample size limits the analysis of how galaxy growth dropouts in the HSC SSP using strict color cuts includ- properties (e.g., star-formation rate (SFR)) depend on ing the requirement that sources are not signi cantly 0 0 properties like stellar mass and environment, especially detected ( < 2) in the g band. This g band could re- at the bright-end. The large cost of spectroscopically move UV-bright AGNs and could explain why Ono et al. surveying faint sources leaves the most ecient method (2018) nd less sources at M < 24 mag than Akiyama of using multi-wavelength photometry as the best way et al. (2018) (see Figure 7 in Ono et al. 2018). to collect larger samples of star-forming galaxies. For Here, we make use of the 24 deg Spitzer -HETDEX example, a few thousand z = 4 galaxies were detected Exploratory Large Area (SHELA) survey dataset in the four 1 deg elds of the CFHT Survey (van der to probe both AGN-dominated and star-formation- Burg et al. 2010). dominated UV-emitting galaxies over a large area. The Another challenge in measuring the bright end of SHELA dataset includes deep (22.6 AB mag, 50% com- the UV luminosity function is the existence of AGNs plete) 3.6 m and 4.5 m imaging from Spitzer /IRAC 0 0 0 0 0 and their photometric similarities with UV-bright galax- (Papovich et al. 2016) and u g r i z imaging from the ies. The spectral energy distributions (SEDs) of AGN- Dark Energy Camera over 18 deg (DECam; Wold et dominated galaxies are characterized by a power-law al. , in prep). Because SHELA falls within SDSS Stripe- continuum and highly ionized emission lines in the rest- 82 there exists a large library of ancillary data, which frame UV (e.g., Stevans et al. 2014), and like high-z we take advantage of by including in our analysis the UV-bright galaxies the observed SEDs of high-z AGNs VISTA J and Ks photometry from the VICS82 survey exhibit a Lyman break feature due to absorption from (Geach et al. 2017) to help rule out low-z interloping intervening neutral hydrogen in the IGM. Thus, any UV- galaxies. In addition, there is deep X-ray imaging in bright galaxy selection technique relying on the Lyman this eld from the Stripe-82X survey (LaMassa et al. break will also select AGNs. Some have attempted to 2016), which could be used to identify bright AGNs. use a morphological cut to break the color degeneracy We select objects at z > 4 based on photometric crite- of UV-bright galaxies and AGNs by assuming the for- ria. Our sample includes, therefore, both galaxies whose mer will appear extended and the latter will be strictly light is powered by star-formation and AGN activity. Galaxy and AGN Luminosity Functions at z = 4 with SHELA 3 As the bulk of AGNs at z  4 are too faint to detect K = 20:9 mag. Figure 1 shows the lter transmission in the existing X-ray data, we include all z = 4 can- curves for the nine photometric bands used overplot- didate galaxies, regardless of powering source, in our ted with model high-z galaxy spectra illustrating the sample and use our large dynamic range in luminosity { wavelength coverage of our dataset. combined with very bright AGNs from SDSS (Richards 2.2. DECam Data Reduction and Photometric et al. 2006, Akiyama et al. 2018) and very faint galax- Calibration ies from the deeper, narrower Hubble Space Telescope surveys (Finkelstein et al. 2015) { to t the luminosity The DECam images were processed by the NOAO functions of both populations simultaneously. Impor- Community Pipeline (CP). A detailed description of the tantly, our sample is selected using both optical and Community Pipeline reduction procedure can be found Spitzer mid-infrared data, which results in an improved in the DECam Data Handbook on the NOAO website , contamination rate over optical data alone. however, we provide a brief summary of the procedure This paper is organized as follows. The SHELA eld here. First, the DECam images were calibrated using dataset used in this paper is summarized in 2.1. The calibration exposures from the observing run. The main DECam reduction are discussed in Sections 2.2{2.4 and calibration steps included an electronic bias calibration, the IRAC data reduction and photometry in 2.5. Sam- saturation masking, bad pixel masking and interpola- ple selection and contamination are discussed in Section tion, dark count calibration, linearity correction, and 3. Our UV luminosity function is presented in Section 4. at- eld calibration. Next, the images were astrometri- The implications of our results are discussed in section 5. cally calibrated with 2MASS reference images. Finally, We summarize our work and discuss future work in Sec- the images were remapped to a grid where each pixel is tion 6. Throughout this paper we assume a Planck 2013 a square with a side length of 0.27 . Observations taken 1 1 cosmology, with H = 67:8 km s Mpc , = 0:307 on the same night were then co-added. 0 M and = 0:693 (Planck Collaboration et al. 2014). All The CP data products for the SHELA eld were down- magnitudes given are in the AB system (Oke & Gunn loaded from the NOAO Science Archive . The data 1983). products include the co-added images, remapped im- ages, data quality maps (DQMs), exposure time maps 2. DATA REDUCTION AND PHOTOMETRY (ETMs), and weight maps (WMs). The co-added im- In this section, we describe our dataset, image reduc- ages from the CP were not intended for scienti c use, tion, and source extraction procedures. The procedures so we opted to co-add the remapped images. We fol- applied to DECam imaging are largely similar to those lowed the co-adding procedure of Wold et al. (in prep.), used with these data in Wold et al. (in prep), thus we which we summarize here. Using the software package direct the reader there for more detailed information. SWARP (Bertin et al. 2002) the sub-images stored in the FITS les of the remapped images were stitched to- 2.1. Overview of Dataset gether and background subtracted. The remapped im- In this study, we use imaging in nine photometric ages were combined using a weighted mean procedure bands spanning the optical to mid-IR in the SHELA optimized for point-sources. The weighting of each im- Field. The SHELA Field is centered at R.A. = age is a function of the seeing, transparency, and sky h m s  0 00 1 22 00 , declination = +0 00 00 (J2000) and extends brightness and is de ned by Equation A3 in Gawiser approximately 6: 5 in R.A. and 1: 25 in declination. et al. (2006) as 0 0 0 0 0 The optical bands consist of u , g , r , i , and z from the Dark Energy Camera (DECam) and covers 17.8 factor PS 2 w = ; (1) deg of the SHELA footprint (Wold et al. in prep). i scale  rms i i The mid-IR bands include the 3.6 m and 4.5 m from the Infrared Array Camera (IRAC) aboard the Spitzer where scale is the image transparency (de ned as the Space Telescope and covers 24 deg (Papovich et al. median brightness of the bright unsaturated stars after 2016). In addition, we include near-IR photometry in J normalizing the brightness measurement of each star by and K from the February 2017 version of the VISTA- s its median brightness across all exposures), rms is the CFHT Stripe 82 Near-infrared Survey (VICS82; Geach root mean square of the uctuations in background pix- et al. 2017), which covers 85% of the optical imag- els, and factor is de ned as ing footprint and has 5- depths of J = 21:3 mag and http://ast.noao.edu/data/docs 1 3 http://stri-cluster.herts.ac.uk/vics82/ http://archive.noao.edu/ 4 Stevans et al. 1.0 J Ks 0.8 ′ ′ i z [4.5] 0.6 [3.6] 0.4 z=3.5 0.2 z=4.5 0.0 0.40 0.75 1.25 2.10 3.60 4.50 Wavelength (μm) Figure 1. The lter transmission curves for the nine photometric bands used in this study (curves are labelled in the gure) and two model star-forming galaxy spectra with redshifts z = 3:5 and z = 4:5, respectively (black). The model spectra have units of Jy and are arbitrarily scaled. The z = 3:5 galaxy spectrum falls completely red-ward of the u band transmission curve and the z = 4:5 galaxy spectra has almost zero ux falling in the g bandpass. Table 1. DECam Imaging Summary - Seeing and Limiting Magnitude 0 0 0 0 0 u g r i z Sub-Field R.A. Dec. FWHM depth FWHM depth FWHM depth FWHM depth FWHM depth 00 00 00 00 00 ID No. (J2000) (J2000) ( ) (mag) ( ) (mag) ( ) (mag) ( ) (mag) ( ) (mag) h m s  0 00 SHELA-1 1 00 52:8 -0 00 36 1.12 25.2 1.06 24.8 1.0 24.8 0.87 24.5 0.86 23.8 h m s  0 00 SHELA-2 1 07 02:4 -0 00 36 1.2 25.2 1.26 24.8 1.28 24.6 1.39 23.9 0.96 23.6 h m s  0 00 SHELA-3 1 13 12:0 -0 00 36 1.22 25.4 1.36 25.0 1.14 24.7 1.04 24.5 1.21 23.5 h m s  0 00 SHELA-4 1 19 21:6 -0 00 36 1.15 25.3 1.4 24.4 1.05 24.3 1.0 22.1 1.13 23.7 h m s  0 00 SHELA-5 1 25 31:2 -0 00 36 1.21 25.1 1.07 24.9 1.02 24.3 0.93 23.9 0.85 23.6 h m s  0 00 SHELA-6 1 31 40:8 -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 1.89"-diameter apertures on the PSF-matched images (see Section 2.3). The initial rms per pixel was de ned as the inverse of the square-root of the exposure time. The median of FWHM stack factor = 1 exp 1:3 ; (2) the rms map is scaled to the global pixel-to-pixel rms FWHM which is de ned as the standard deviation of the uxes in good-quality, blank sky pixels. Good-quality, blank where FWHM is the median FWMH of bright un- stack sky pixels are pixels not included in a source according saturated stars in an unweighted stacked image and to our initial Source Extractor catalog (see Section 2.3 FWHM is the median FWHM of bright unsaturated for discussion of our source extraction procedure), and stars in each individual exposure. have an exposure time greater than 0.9 times the median The seeing and transparency measurements were de- value. termined using a preliminary source catalog generated The DECam imaging data were photometrically cal- for each resampled image using the Source Extractor ibrated with photometry from the Sloan Digital Sky software package (Bertin & Arnouts 1996).The seeing Survey (SDSS) data release 11 (DR11; Eisenstein et al. in the nal stacked images are listed in Table 1. 2011) using only F0 stars. F0 stars were used because After discovering the original WMs from the Com- their spectral energy distribution span all ve optical munity Pipeline had values inconsistent with the ETM, lters while appearing in the sky at a suciently high we created custom rms maps for the co-added images. Transmission Galaxy and AGN Luminosity Functions at z = 4 with SHELA 5 surface density to provide statistically signi cant num- deconvolution routine with an increasing number of it- bers in each DECam image. We began by creating a erations until the PSF of the convolved image (again preliminary source catalog for the stacked DECam im- measured from stacking stars) had a ux within a 7- ages using Source Extractor and position matching to pixel (1.89 ) diameter aperture matched to that of the the SDSS source catalog. Then we selected F0 stars us- PSF of the target image to within 5%. The total uxes ing SDSS colors by integrating an F0-star model spec- were measured in 30-pixel (8.1 ) diameter circular aper- trum from the 1993 Kurucz Stellar Atmospheres Atlas tures. In Figure 2 we show the results of PSF-matching (Kurucz 1979) with each of the ve optical SDSS lter in each sub- eld by displaying a comparison of the en- curves. For sources in the catalog to be identi ed as larged PSFs to the largest PSF as a percent di erence. an F0 star, the total color di erences, using colors for Due to variations in intrinsic galaxy colors and vari- all adjacent bands, added in quadrature must have been able image depth and sky coverage, some galaxies will less than 0.35. We then calculate the expected mag- not appear in all bands. To get photometric measure- nitude o set between SDSS and DECam lters for F0 ments of all sources in every DECam image we com- 0 0 0 stars, which are as follows: u : 0.33, g : 0.02, r : -0.001, bined the information in the ve optical band images 0 0 i : -0.02, and z : -0.01. The zero-point for each lter into a single detection image. We followed the proce- was then calculated as dure of Szalay et al. (1999) and summed the square of the signal-to-noise ratio in each band pixel-by-pixel as AB ZPT = median(m m m ); (3) DECam o set SDSS follows: band;i where m is the expected magnitude o set between o set D = ; (4) band;i SDSS and DECam lters for F0 stars. After the zero- points were applied to each stacked image, the image where D is the detection image ith pixel value, F i band;i pixel values were converted to units of nJy. is the ith pixel ux in the band image, and  is band;i the rms at that ith pixel pulled from the the band rms 2.3. DECam Photometry image. A weight image associated with this detection Studying galaxy properties relies on accurate mea- image was created where pixels associated with detec- surements of galaxy colors. One way to obtain ac- tion map pixels with data in at least one band have a curate colors is to perform xed-aperture photometry value of unity and pixels associated with detection map where you measure source uxes in every band using pixels without data have a value of zero. the same sized aperture. However, since the DECam im- Photometry was measured on the PSF-matched im- ages where taken in varying seeing conditions they have ages using the Source Extractor software (v2.19.5, point spread functions (PSFs) with a range of full width Bertin & Arnouts 1996). Catalogs were created for at half maximum (FWHM). To perform xed-aperture each of the six SHELA sub- elds with Source Extractor photometry on these images, the PSFs of all the imaging in two image mode using the detection image described covering a single patch of sky must be adjusted to have above, and cycling through the ve DECam bands as a similar PSF size. We divided the DECam imaging the measurement image. In the nal source catalog, we into six sub- elds (each de ned as roughly one DECam maximized the detection of faint sources while minimiz- pointing). For each sub- eld, we enlarged the PSFs of ing false detections by optimizing the combination of the the stacked images to match the PSF of the stacked im- SExtractor parameters DETECT THRESH and DE- age with the largest PSF in that sub- eld. For example, TECT MINAREA. We did this by running SExtractor 0 0 in sub- eld SHELA-1 we matched the PSFs of the g , r , with an array of combinations of DETECT THRESH 0 0 0 i , and z stacked images to the PSF of the u band. To and DETECT MINAREA and chose the combination enlarge the PSFs we adopted the procedure of Finkel- of 1.6 and 3, respectively, which detected all sources stein et al. (2015) who used the IDL deconv tool Lucy- that appeared real by visual inspection and included Richardson deconvolution routine. This routine takes as the fewest false positive detections from random noise inputs two PSFs (the desired larger PSF and the starting uctuations. smaller PSF) and the number of iterations to run and We measure source colors in 1.89 diameter circu- outputs a convolution kernel. The input image PSFs lar apertures (which corresponds to an enclosed ux were produced by median combining the 100 brightest fraction of 59-75% for unresolved sources in our sub- stars (sources with stellar classi cations in SDSS DR11) elds). To obtain the total ux, we derived an aper- in each image. Before combining, the stars were over- ture correction de ned as the ux in a 1.89 diame- sampled by a factor of 10, re-centered, and then binned ter aperture divided by the ux in a Kron aperture by ten to ensure the star centroids aligned. We ran the (i.e., FLUX AUTO), using the default Kron aperture 6 Stevans et al. 20 20 20 g' u' u' Field: 1 Field: 2 Field: 3 15 15 15 r' g' r' Matched to: u' Matched to: i' Matched to: g' i' r' i' 10 10 10 z' z' z' 5 5 5 0 0 0 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 Aperture Diameter (arcsec) Aperture Diameter (arcsec) Aperture Diameter (arcsec) 20 20 20 u' g' u' Field: 4 Field: 5 Field: 6 15 15 15 r' r' r' Matched to: g' Matched to: u' Matched to: g' i' i' i' 10 10 10 z' z' z' 5 5 5 0 0 0 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 Aperture Diameter (arcsec) Aperture Diameter (arcsec) Aperture Diameter (arcsec) Figure 2. The results of PSF-matching the DECam images. Each panel shows the percent di erence between the enlarged PSFs and the largest PSF per SHELA sub- eld. The colored lines correspond to the four bands listed in each panel's legend. The vertical dashed line denotes 1.89" which is the aperture diameter at which we compared PSFs during the PSF-matching procedure (see Section 2.3 for details). The horizontal dashed line was placed to zero to guide the eye. This gure illustrates that all bands in all sub- elds have PSFs that collect the same fraction of light as their respective largest PSF to within 5% except for the z band in sub- eld 3 which matches to about 6%. parameters of PHOT AUTOPARAM= 2.5, 3.5, which compared the extinction-corrected DECam photometry has been shown to calculate the total ux to within 5% to SDSS DR14 per sub- eld and band and found agree- (Finkelstein et al. 2015). This correction was derived in ment for point sources to better than m < 0:05 0:2 the r -band on a per-object basis to account for di erent mag in terms of scatter. source sizes and ellipticities and was applied to the uxes in the other DECam bands per sub- eld. In areas where 2.4. DECam Photometric Errors the sub- elds overlapped, sources with positions that We estimated photometric uncertainties in the DE- matched to within 1.2 in neighboring sub- eld catalogs Cam images by estimating the image noise in aper- had their uxes mean-combined after being weighted by tures as a function of pixels per aperture, N, following the inverse square of their uncertainties (see Section 2.4). the procedure described in Section 2 of Papovich et al. The DECam source uxes were corrected for Galactic (2016). There are two limiting cases for the uncertainty extinction using the color excess E(B-V) measurements in apertures with N pixels,  . If pixel errors are com- by Schla y & Finkbeiner (2011). We obtained E(B-V) pletely uncorrelated, the aperture uncertainty scales as values using the Galactic Dust Reddening and Extinc- the square root of the number of pixels,  =   N , N 1 tion application on the NASA/IPAC Infrared Science where  is the standard deviation of sky background 4 1 Archive (IRAS) website . We queried IRAS for E(B-V) pixels. If pixel errors are completely correlated then values (the mean value within a 5 radius) for a grid =   N (Quadri et al. 2006). Thus, the aperture 0 N 1 of points across the SHELA eld with 4 spacing and uncertainty will scale as N with 0:5 < < 1. assigned each source the E(B-V) value from the closest To estimate the aperture noise as a function of N grid point. The Cardelli et al. (1989) Milky Way redden- pixels, we measured the sky counts in 80,000 randomly ing curve parameterized by R = 3:1 was used to derive 00 00 placed apertures ranging in diameter from 0: 27 to 8: 1 the corrections at each band's central wavelength. We across each stacked DECam image. We required aper- tures to fall in regions of the background sky, which we http://irsa.ipac.caltech.edu/applications/DUST/ de ne as the region where the exposure time map has the value of the at least the median exposure time (en- % Diff. in Total Fractional Flux % Diff. in Total Fractional Flux % Diff. in Total Fractional Flux % Diff. in Total Fractional Flux % Diff. in Total Fractional Flux % Diff. in Total Fractional Flux Galaxy and AGN Luminosity Functions at z = 4 with SHELA 7 Table 2. Fit Parameters for Background Fluctuations as Function β = 1 of Aperture Size Using Eq.(5) y 800 i Sub-Field ID Band   rms 1 med (nJy) (nJy) SHELA-1 u 5.02 0.09 0.89 1.33 0.35 5.58 β = 0.5 g 4.07 0.21 0.83 1.71 0.43 8.85 σ r 4.14 0.18 0.91 2.37 0.36 10.2 i 4.1 0.13 1.0 2.88 0.42 14.5 200 z 7.16 0.19 0.93 2.55 0.45 25.0 SHELA-2 u 1.21 0.24 0.91 2.41 0.55 5.88 0 0 5 10 15 20 25 g 1.8 0.35 0.87 2.3 0.5 8.26 r 3.02 0.21 0.94 2.68 0.41 10.2 i 14.3 0.05 1.0 1.37 0.36 15.8 0 Figure 3. Background noise uctuations,  , in an aperture z 5.36 0.34 0.9 2.48 0.51 27.7 of N pixels plotted as a function of the square root of the SHELA-3 u 1.68 0.23 0.89 2.2 0.42 4.81 number of pixels for the ve DECam band images in sub- g 5.36 0.12 0.88 1.33 0.33 6.25 eld SHELA-1. The colored dotted lines are the measured r 2.77 0.19 0.94 2.63 0.45 9.8 aperture uxes and the solid lines are ts to the data. See legend insert for color coding information. The dot-dashed i 3.06 0.32 0.92 3.07 0.38 12.6 0 line shows the relation assuming uncorrelated pixels, z 11.0 0.15 0.91 1.95 0.41 24.1 p N . The dashed line shows the relation assuming perfectly SHELA-4 u 1.08 0.26 0.95 2.54 0.49 4.64 correlated pixels (  N ; Quadri et al. 2006). g 7.36 0.13 0.85 1.24 0.34 8.31 r 2.86 0.25 0.91 2.44 0.51 12.1 measured ux uncertainty in a given aperture,  , as i 24.9 1.42 0.6 1.27 0.6 63.5 a function of the square root of the number of pixels z 5.48 0.26 0.93 2.79 0.45 22.6 in each aperture, N , for the ve DECam bands in the SHELA-5 u 5.25 0.11 0.87 1.25 0.33 5.8 sub- eld SHELA-1. g 3.11 0.23 0.87 1.92 0.43 7.86 Following Papovich et al. 2016 we t a parameterized r 4.78 0.38 0.81 2.02 0.46 15.8 function to the noise in an aperture of N pixels,  , as, i 5.53 0.33 0.85 2.32 0.46 21.9 z 6.86 0.23 0.92 2.53 0.49 28.0 =  ( N + N ); (5) N 1 SHELA-6 u 1.66 0.15 0.95 2.38 0.42 4.6 g 6.15 0.12 0.86 1.27 0.32 6.8 where  is the pixel-to-pixel standard deviation in the r 4.81 0.21 0.87 2.09 0.4 12.1 sky background, and , , , and  are free parameters. i 6.87 0.14 0.93 2.14 0.37 14.7 The best- tting parameters in Equation 5 for the com- bined DECam images are listed in Table 2. While the z 6.09 0.51 0.82 2.13 0.53 32.8 second term was intended to aid in tting the data at Note| Typical values of  0.65-0.70 when using a two parameter large N values, in actuality the second term contributed t (i.e. with only and ) suggest slightly correlated noise between signi cantly at all N values resulting in 0.9-1. Nev- pixels. ertheless our functional ts reproduce the data well as can be seen, for example, in Figure 3. To estimate how correlated the pixel-to-pixel noise is, we t  with only the rst term in Equation 5 and found typical values suring >50% of each image was considered), excluding of  0.65-0.70 suggesting slightly correlated pixel-to- detected sources and pixels agged in the DQM. We also pixel noise. required the apertures do not overlap with each other. The photometric errors estimated by Equation 5 were We then estimated  for each aperture with N pixels scaled to apply to ux measurements outside the region by computing the standard deviation of the distribution with the median exposure time. The ux uncertainty for of aperture uxes from the normalized median absolute the i-th source in band b in the sub- eld f is calculated deviation,  (Beers et al. 1990). We calculated nmad 1 as by computing  for all pixels in the background sky ! nmad rms as de ned above. Figure 3 shows an example of the i;f;b 2 2 =  ; (6) i;f;b N;f;b rms med;f;b (nJy) N 8 Stevans et al. where  is given by Equation 5 for each band and uxes is closest to the actual image pixels, with respect N;f;b sub- eld, rms is the value of the rms map at the to the noise model. We describe how we use the Tractor i;f;b central pixel location of the i-th source in each band to perform forced photometry on IRAC images in detail and sub- eld, and rms is the median value of the below. med;f;b rms map in each band and sub- eld. The photometric We begin our source modeling procedure by selecting error estimates exclude Poisson photon errors, which we the ducial band high-resolution model of each source. estimate to contribute <5% uncertainty to the optical We use the uxes and surface brightness pro le shape uxes of our high-z candidates. parameters measured in our DECam detection image be- As described in Section 2.3, all source uxes are mea- cause the image combines the information of all sources sured in circular apertures of 1.89 diameter and scaled in the ve optical band images (as described in Section to total on a per-object basis. Likewise, the ux uncer- 2.3). Second, we use one-component circular Gaussian tainties in the apertures were scaled by the same amount to model the PSF. During the modeling of each source, to determine the total ux uncertainties, so that the S/N we allow Tractor to optimize the Gaussian  value, in for the total ux is the same as for the aperture ux. addition to optimize a source ux. We nd that the me- dian of the optimized Gaussian  is 0: 80 (equivalent to a full-width at half maximum of 1: 88) for both IRAC 2.5. IRAC data reduction and photometry 3.6 and 4.5 m images, consistent with those measured As discussed below, we wish to enhance the validity of from an empirical point response function for the 3.6m our z = 4 galaxy sample by including IRAC photometry and 4.5m image (FWHM of 1: 97, see Papovich 2016, in our galaxy sample selection. While we could allow Section 3.4). this by position-matching the published Spitzer/IRAC In practice, we extracted an IRAC image cutout of catalog from Papovich et al. (2016) to our DECam cat- each source in the input DECam catalog. We selected alog, this is not optimal for two reasons. First, the Pa- 00 00 the cutout size of 16  16 . This cutout size represents povich et al. (2016) catalog is IRAC-detected, and so a trade-o between minimizing computational costs re- only includes sources with signi cant IRAC ux, while lated to larger cutout sizes and ensuring that the sources for our purposes, even a non-detection in IRAC can be lie well within the cutout extent. The sources of interest useful for calculating a photometric redshift. Second, within cutout are modeled as either unresolved (i.e., a this catalog uses apertures de ned on the positions and point source) or resolved based on the DECam detec- shapes of the IRAC sources, while the larger PSF of tion image. We considered a source to be resolved if the IRAC data results in signi cant blending, which is an estimated radius r > 0: 1. We de ne the radius as, a larger issue at fainter magnitudes, where we expect r = a  b=a, where a is a semi-major axis and source to nd the bulk of our sources of interest. For these b=a is an axis ratio. We perform the photometry for reasons, we applied the Tractor image modeling code resolved sources using a deVaucouleurs pro le (equiva- (Lang et al. 2016a,b) to perform \forced photometry", lent to S ersic pro le with n=4) with shape parameters which employs prior measurements of source positions (semi-major axis, position angle, and axis ratio) mea- and surface brightness pro les from a high-resolution sured using our DECam detection image. We have also band to model and t the uxes of the source in the re- performed the photometry using an exponential pro le maining bands, splitting the ux in overlapping objects (equivalent to S ersic pro le with n=1), but we do not into their respective sources. We speci cally used the nd any signi cant di erence between the IRAC ux Tractor to optimize the likelihood for the photometric measurements for the two galaxy pro les. Therefore, properties of DECam sources in each of IRAC 3.6 and we adopt a deVaucouleurs pro le to model all resolved 4.5 m bands given initial information on the source sources. The Tractor simultaneously modeled and op- and image parameters. The input image parameters of timized the sources of interest and neighboring sources IRAC 3.6 and 4.5 m images included a noise mode, within the cutout. Finally, the Tractor provided the a point spread function (PSF) model, image astromet- measurement IRAC ux of each DECam source with ric information (WCS), and calibration information (the the lowest reduced chi-squared value. We validated the \sky noise" or rms of the image background). The in- Tractor-based IRAC uxes by comparing the uxes of put source parameters included the DECam source posi- isolated sources (no neighbors within 3 ) to the pub- tions, brightness, and surface brightness pro le shapes. lished Spitzer/IRAC catalog from Papovich et al. (2016). The Tractor proceeds by rendering a model of a galaxy For both bands, we found good agreement with a bias or a point source convolved with the image PSF model o set of m < 0:05 mag and a scatter of <0.11 mag at each IRAC band and then performs a linear least- square t for source uxes such that the sum of source Galaxy and AGN Luminosity Functions at z = 4 with SHELA 9 down to m = 20:5 mag and a bias o set of m < 0:13 overcome this we required a source to pass the z se- phot mag and a scatter of <0.26 down to m = 22 mag. lection process with and without including the VISTA J and K photometry. This double z selection pro- phot 3. SAMPLE SELECTION cess resulted in an initial sample of 4,364 high-z objects. Next, we performed an investigation into possible con- 3.1. Photometric Redshifts and Selection Criteria tamination, which resulted in additional selection crite- We selected our sample of high-redshift galaxies using ria and a re ned sample. a selection procedure that relies on photometric redshift (z ) tting, which leverages the combined information phot 3.2. Identifying Contamination in all photometric bands used. We obtained z 's and phot Photometric studies of high-z objects can be contam- z probability distribution functions (PDFs) from the phot inated by galactic stars and low-z galaxies whose 4000 EAZY software package (Brammer et al. 2008). For A break can mimic the Lyman- break of high-z galax- this analysis, we use the \z a" redshift column from 2 ies. The inclusion of IRAC photometry is crucially im- EAZY which is produced by minimizing the  in the portant for removing galactic stars from our sample as all-template linear combination mode. We also tried the galactic stars have optical colors very similar to z = 4 \z peak" column and our resulting luminosity function objects (Figure 5). While inspecting the photometry did not change signi cantly. EAZY assumes the inter- and best- tting SED templates to con rm our ts were galactic medium (IGM) prescription of Madau (1995). robust and our high-z galaxy candidates were convincing We did not use any magnitude priors based on galaxy we found evidence of contamination in our preliminary luminosity functions when running EAZY as the exist- photometrically selected sample. We explored ways of ing uncertainties in the bright-end would bias our re- identifying and removing the contamination including sults (e ectively assuming a at prior). We then ap- cross-matching our sample with proper motion catalogs, plied a number of selection criteria using the z PDFs phot SDSS spectroscopy, and X-ray catalogs. Additionally, from EAZY following the procedure of Finkelstein et al. we implemented machine learning methods. (2015), which we summarize here, to construct a z 4 galaxy sample. First, we required a source to have a 3.2.1. Cross-matching with NOMAD and SDSS signal-to-noise ratio of greater than 3.5 in the two pho- 0 0 tometric bands (r and i bands) that probe the UV While inspecting the photometry of our preliminary continuum of our galaxy sample, which has been shown sample derived before the inclusion of the VISTA data 0 0 to limit contamination by noise to negligible amounts by we found a fraction of candidates had red r i colors 0 0 Finkelstein et al. (2015). Next we required the integral and excesses in the i and z bands with respect to the of each source's z PDF for z > 2:5 to be greater than best tting z  4 template. We investigated whether phot 0.8 and the integral of the z PDF from z = 3.5-4.5 these objects had low redshift origins by cross-matching phot to be greater than the integral of the z PDF in all our sample with the Naval Observatory Merged As- phot other z = 1 bins centered on integer-valued redshifts. trometric Dataset (NOMAD) proper motion catalog Finally, we required a source to have no photometric (Zacharias et al. 2004) and the SDSS spectral catalog 0 0 0 0 ags on its u , g , r , and i ux measurements in our (DR13; Albareti et al. 2017). The cross-matching with photometric catalog. These four bands are crucial for NOMAD resulted in the identi cation of 16 objects with probing both sides of the Lyman break of galaxies in proper motion measurements, 6 of which were large in the redshift range of interest. Photometric ags indi- magnitude suggesting that these sources were stellar cate saturated pixels, transient sources, or bad pixels as contaminants. Cross-matching with SDSS spectra re- de ned by the CP (see section 2.2). sulted in the identi cation of 23 z  4 QSOs, two low Initially, we required sources to pass this z selec- mass stars, and one low-z galaxy in our sample. All of phot 0 0 tion procedure when only the DECam and IRAC bands the stellar objects and the low-z galaxy had red r i were used. After inspecting candidates from this selec- colors con rming our suspicion that a fraction of our tion and nding low-z contaminants (see Section 3.2) we sample had low-z origins. We removed from our candi- moved to include the VISTA J and K photometry when date sample the 6 objects with proper motions above 50 available. The candidate sample of high-z galaxies after mas/yr in NOMAD. We chose the threshold 50 mas/yr including the VISTA J and K photometry signi cantly because some SDSS QSOs are reported to have small reduced the amount of obvious low-z contamination as nonzero proper motions in the NOMAD catalog. After described in the following section, however, a new class inclusion of the VISTA J and K photometry, 5 of the of contaminant appeared in the candidate sample due 6 rejected objects became best- t by a low-z solution to erroneous VISTA photometry for some sources. To and therefore were rejected by our selection process au- 10 Stevans et al. tomatically. Results like this suggested the inclusion of measurements. The obviously high-z objects and the the VISTA data improved the delity of our sample. SDSS QSOs had roughly at rest-UV spectral slopes 0 0 0 (i.e., the i , r , and z uxes were comparable) while the 3.2.2. Cross-matching with X-ray catalog Stripe 82X obviously low-z objects had a clear spectral peak be- 0 0 In principle, AGNs can be distinguished from UV- tween the optical and mid-IR, speci cally the r i and 0 0 0 bright star-forming galaxies at high z by their X-ray i z colors were quite red while the z [3:6] was blue. emission as AGNs dominate X-ray number counts down We then selected three data "features" or measurements 17 2 1 to F = 1  10 erg cm s in 0.5-2 keV Chan- for the decision tree to choose from to predict the clas- dra imaging (Lehmer et al. 2012). In fact, Giallongo si cations of the training set: 1. The  of the best 0 0 et al. (2015) found 22 AGN candidates at z > 4 by tting high-z galaxy template, 2. The r i color, and measuring uxes in deep 0.5-2 keV Chandra imaging 3. The signal-to-noise ratio in the u band. We include at the positions of their optically-selected candidates. this signal-to-noise ratio because most obviously high- They required AGN candidates to have an X-ray detec- z sources had a signal-to-noise ratio of less than 3 (as 17 2 1 tion of F > 1:5 10 erg cm s . We attempted they should, as this lter is completely blue-ward of the to identify AGNs in our candidate sample by position Lyman break for the target redshift range) while the ob- matching with the 31 deg Stripe 82X X-ray catalog viously low-z source did not. We restricted the training (LaMassa et al. 2016). Unfortunately, the ux limit set to include only the obviously high-z objects (includ- 16 2 1 at 0.5-2 keV was 8:7  10 erg cm s , shallower ing SDSS classi ed QSOs) and obviously low-z objects. than the detection threshold used by Giallongo et al. To evaluate the decision tree performance we assigned (2015). Cross-matching the Stripe 82X X-ray Catalog 34% of the training set to a test set and left the test set with our candidate sample resulted in 8 matches within out of the rst round of training. The rst round of 7 or the matching radius used by LaMassa et al. (2016) training was validated using three-fold cross-validation to match the Stripe 82X X-ray Catalog with ancillary with a score of 94+/-2%. Three-fold cross validation data. Two of the X-ray-matched sources (m 19) are veri es the model is not over- tting or overly dependent SDSS AGNs, although one has two extended objects on a small randomly selected training (or validation) set. within the 7 matching radius, drawing into question The process of three-fold cross-validation involves split- the likelihood of the match. Another X-ray-matched ting the training sample into three sub-samples, train- 0 0 source (m 22) is a very red object (r i = 0.9 and ing the model independently on each combination of i [3:6] =3.4) without SDSS spectroscopy. The remain- two sub-samples while validating on the remaining sub- ing ve X-ray-matched sources are fainter (m 24) and sample, and then averaging the validation scores of all without SDSS spectroscopy, and two have large separa- the combinations. A validation score of 100% indicates tions (> 6 ) with their X-ray counterpart. All 8 X-ray- each sub-sample combination trained a model that suc- matched sources satis ed each of our selection criteria cessfully predicted the classi cations of every object in and made it into our nal sample. the remaining sub-sample. After validation, we applied the model to the test set and achieved a test score of 3.2.3. Insights from machine learning 91%. We then incorporated the test set into the train- To further understand and minimize the contamina- ing set and retrained the model. This second round tion in our preliminary sample we utilized two machine of training had a 3-fold cross validation score of 94 learning algorithms: a decision tree algorithm and a ran- 5%. The algorithm determined that classi cation of the dom forest algorithm. First, we tried the decision tree training set could be predicted at 92% accuracy using algorithm from the sci-kit learn Python package (Pe- 0 0 2 the r i color and the  only, where obviously high-z dregosa et al. 2011). Executing the decision tree algo- 0 0 2 objects have r i < 0.555 and  < 70.6. rithm involved creating a training set by classifying (via After performing the decision tree machine learning visual inspection) the 300 brightest objects as one of we tried the random forest algorithm by using the scikit- ve types: 1. obviously high-z galaxy or high-z AGN, learn routine RandomForestClassi er (Pedregosa et al. 2. obviously low-z galaxy, 3. indistinguishable between 2011). The random forest algorithm uses the combined high-z or low-z object, 4. SDSS spectroscopically iden- information of all the features of our dataset instead ti ed QSO, 5. spurious source. The classi cation was of using only the three well-motivated features that we driven by a combination of the shape of the optical SED, provided the decision tree algorithm. The random for- the  of the best tting high-z (z  4) galaxy tem- est algorithm ts a number of decision tree classi ers on plate, the di erence in  between the best tting high-z a randomly drawn and bootstrapped subset of the data galaxy template and the best tting low-z galaxy tem- using a randomly drawn subset of the dataset features. plate, SDSS spectral classi cation, and proper motion Galaxy and AGN Luminosity Functions at z = 4 with SHELA 11 The decision tree classi ers are averaged to maximize ac- duced the source uxes randomly to create a at dis- curacy and control over- tting. To provide the best clas- tribution of dimmed r -band mag spanning the range si cations to the algorithm we re-inspected the photo- of our candidate high-z galaxies (m = 18-26). We metry and the besting- tting EAZY template SEDs of assigned ux uncertainties to the dimmed uxes by ran- the brightest 311 objects (m < 22:9) and classi ed by domly drawing from the ux uncertainties of our can- eye each with a probability of being at high z, at low z, didate high-z galaxy sample. We then perturbed the and a spurious source. We also inspected the best- tting dimmed uxes simulating photometric scatter by draw- EAZY template at the redshift of the second highest ing ux perturbations from a Gaussian distribution with peak in the redshift PDF, which was usually between a standard deviation equal to the assigned ux uncer- z = 0:1 0:5. The features we provided the algorithm tainties. Since the VISTA J and K photometry from included all colors using the DECam and IRAC pho- the VICS82 Survey covered 85% of the total SHELA tometric bands, the  value of the best- tting EAZY survey area, 15% of our high-z galaxy candidates do template, and the u-band S/N ratio. We ran the random not have ux measurements in the VISTA bands. We forest algorithm using 1000 estimators (or decision trees) incorporated this property of our catalog into the mock with a max depth of two and balanced class weights. catalog by randomly omitting dimmed J and K uxes While the classi cation accuracy of the random forest at the same rate as the fraction of missing J and K was only marginally better than the decision tree, we uxes in our nal sample for a given m 0 bin. were able to determine the relative importance of the We created 4.8 million arti cially dimmed sources and features within the random forest. The most important ran them through EAZY and our selection criteria in an feature was the i [3:6] color, which was consistent identical manner as our real catalog. We found 30,594 with our prior observations of the low-z contaminants arti cially dimmed sources satis ed our z = 4 selec- phot having blue z [3:6] while the SDSS classi ed AGNs tion criteria. We inspected the undimmed and unper- had at or red z [3:6]. We compared the predictive turbed properties of these sources and found a range of 0 0 0 0 power of the random forest algorithm to that of a simple r i colors (0 < r i < 1:4), and the objects with the 0 0 0 i [3:6] color cut{where high-z candidates were required reddest r i colors contaminated at the highest rate. 0 0 0 to have an i [3:6] > 0:2{and found the simple color We de ned the contamination rate in the i-th r i cut to be just as predictive. We, therefore, elected to color bin as adopt the simple i [3:6] color cut as an additional dimmed,selected,i step in our selection process. After we applied this cut R = ; (7) our high-z object sample consisted of 3,833 objects. We N dimmed,i then inspected the brightest 3,200 brightest candidates where N is the number of dimmed dimmed,selected,i (m 0 < 24:0) and removed 61 spurious sources cutting 0 0 sources satisfying our z = 4 criteria in the i-th r i phot our candidate sample to 3,772 objects. The nal step in 0 0 color bin and N is the total number of dimmed our candidate selection process was an r i color cut, dimmed,i 0 0 sources in the i-th r i color bin. We created seven which was a result of the contamination test described 0 0 0 0 r i color bins spanning the range of r i colors in the following subsection. 0 0 recovered (0 < r i < 1.4) with each bin having width=0.2 mag. The contamination fraction, F , in our 3.3. Estimating Contamination Using Dimmed Real high-z galaxy sample was de ned as Sources R N P i We estimated the contamination in our high-z galaxy total;i (1R ) sample by simulating faint and low-S/N low-z interloper F = ; (8) galaxies following a procedure used by Finkelstein et al. (2015). We did this by selecting a sample of bright where R is the contamination rate de ned by Equation low-z sources from our catalog, dimming and perturbing 7, N is the total number of sources in our object total;i 0 0 their uxes, and assigned the appropriate uncertainties catalog in the i-th r i color bin with 0:1 < z < 0:6 phot to their dimmed uxes. This empirical test assumes that and N is the number of high-z candidates in our nal bright low-z galaxies have the same colors as the faint sample in a given redshift bin. We calculated F as a 0 0 0 lower-z galaxies. The sources we dimmed all had m = function of m and m and found contamination frac- r r i 14.3 - 18 mag, brighter than our brightest z =3.5-4.5 tions of between 0-20% generally increasing as a function candidate source, a best- tting z of 0.1 < z < of magnitude and not exceeding 25% in any magnitude phot phot 0.6, and no photometric ags in any optical band (see bin above our 50%-completeness limit (m 0 > 23:5). We section 2.2 for de nition of photometric ag). We re- learned we can improve our delity by implementing a 12 Stevans et al. Table 3. Summary of z = 4 Sample Selection Criteria 3.4. Completeness Simulations A crucial component of calculating the UV luminosity Criterion Section Reference function is measuring the e ective volume of the survey. S/N > 3:5 (3.1) The e ective volume measurement depends on quanti- S/N > 3:5 (3.1) fying the survey incompleteness due to image depth and PDF(z) dz > 0:8 (3.1) selection e ects. To quantify the survey incompleteness 2:5 4:5 PDF(z) dz > All other z = 1 bins (3.1) we simulated a diverse population of high-z galaxies with 3:5 i [3:6] > 0:2 (3.2.3) assigned photometric properties and uncertainties con- 0 0 r i < 1:0 (3.3) sistent with our source catalog and measured the frac- tion of simulated sources that satis ed our high-z galaxy selection criteria as a function of the absolute magnitude and redshift. The simulated mock galaxies were given properties 0 0 r i < 1:0 color cut, given that the objects with the 0 0 drawn from distributions in redshift and dust atten- reddest r i colors contaminated at the highest rate. uation (e.g., E[B-V]) while the ages and metallicities We then re-ran our contamination simulation with our were xed at 0.2 Gyr and solar (Z = Z ), respectively. nal set of selection criteria and found improvement by Because the fraction of recovered galaxies per redshift a few percentage points in two of our brighter magnitude bin is independent of the number of simulated galax- bins (m = 18.75, 20.5) and, again, contamination frac- ies per redshift bin as long as low-number statistics are tions of between 0-20% generally increasing as a function avoided, the redshift distribution was de ned to be at of magnitude and not exceeding 25% in any magnitude from 2 < z < 6, and the E(B-V) distribution was de ned bin brighter than our 50%-completeness limit as shown 0 0 to be log-normal spanning 0 < E(B-V) < 1 and peaking in Figure 6. The addition of the r i < 1:0 color cut re- at 0.2. Mock SEDs were then generated for each galaxy duced our high-z candidate sample size from 3,772 to the using pythonFSPS (Foreman-Mackey et al. 2014), a nal size of 3,740 with a median z of 3.8. The mea- phot 0 python package that calls the Flexible Stellar Popula- sured i -band magnitude and the redshift distribution tion Synthesis Fortran library (Conroy et al. 2009; Con- of our nal sample is shown in Figure 4. A summary of roy & Gunn 2010). We then integrated each galaxy our sample selection criteria is listed in Table 3. Bright 0 0 0 0 SED through our nine lters (DECam u , g , r , i , and (m 0 < 22) candidates are plotted in Figure 5 showing z ; VISTA J and K ; and IRAC 3.6 m and 4.5 m). the distinct color parameter space they occupy (along Each set of mock photometry was then scaled to have with SDSS spectroscopically classi ed AGNs) compared an r-band apparent magnitude within a log-normal dis- to the parameter space occupied by SDSS spectroscopi- tribution spanning 18 < m 0 < 27. This distribution cally classi ed stars in the SHELA eld. ensured we were simulating the most galaxies at the fainter magnitudes where we expected to be incomplete 4 3 and fewer at bright magnitudes where we expected to be 10 10 very complete. Flux uncertainties and missing VISTA 2 uxes were assigned in the same way as during the con- 2 tamination test in Section 3.3. Mock galaxies were not assigned photometric ags. This likely results in an in- completeness of only a couple percent as the fraction of all sources a ected by ags is small. All sources with a signal-to-noise ratio of greater than 3.5 in any one band, 20 22 24 26 3.0 3.5 4.0 4.5 m are agged in another band less than 2% of the time on i Redshift average (never exceeding 4%). In addition all images have 0.5% of all pixels agged on average (never ex- Figure 4. The i -band magnitude distribution (left) and the ceeding 2%). photometric redshift distribution (right) of our nal sample of z  4 candidates. Photometric redshifts are the redshifts The photometric catalog of 600,000 mock galaxies was where  is minimized for the all-template linear combina- then run through EAZY to generate z PDFs and our phot tion mode from the EAZY software (z a). We have found high-z galaxy selection was applied. The completeness high-z sources over a wide range of brightnesses and across the entire z =3.5-4.5 range. http://dfm.io/python-fsps/current/ Number Number Galaxy and AGN Luminosity Functions at z = 4 with SHELA 13 3.0 2 AGN (3. 2 ≤ z < 3. 5) AGN (3. 2 ≤ z < 3. 5) spec spec AGN (4. 5 ≤ z < 4. 8) AGN (4. 5 ≤ z < 4. 8) spec spec 2.5 2.0 1.5 1.0 −1 Final Candidates Final Candidates 0.5 (m ′ < 22) (m ′ < 22) i i Candidates omitted Candidates omitted ′ ′ b colors (m < 22) b colors (m < 22) i i −2 0.0 ′ ′ Stars (S/N > 100) Stars (S/N > 100) i i AGN (3. 5 ≤ z < 4. 5) AGN (3. 5 ≤ z < 4. 5) spec spec −0.5 −0.5 0.0 0.5 1.0 1.5 2.0 2.5 −0.5 0.0 0.5 1.0 1.5 2.0 2.5 ′ ′ ′ ′ r − i [mag] r − i [mag] 0 0 0 0 Figure 5. Left: Color-color plot showing DECam g r vs r i . Blue points correspond to bright (S=N > 100) sources classi ed as stars within SDSS DR14. Sources spectroscopic identi ed as QSO within SDSS DR13 at 3:2  z < 3:5, 3:5  z < 4:5, and 4:5  z < 4:8 are orange circles, black \x"s, and green circles, respectively. Bright candidate objects from our study are 0 0 0 0 shown as squares (see legend insert for color coding). Right: Color-color plot showing i [3:6] vs r i . The r i < 1:0 and the i [3:6] > 0:2 color selection criteria are denoted as the vertical dashed line and the horizontal dashed line, respectively. No candidates are hidden by the legend inserts. These plots illustrate how the inclusion of IRAC photometry breaks the optical color degeneracy of z = 4 sources and galactic stars. In the left panel where only optical colors are used the z  4 AGNs share the parameter space of galactic stars, while in the right panel where one color includes the [3:6] IRAC band there is a larger separation. was de ned as the number of mock galaxies recovered magnitude at rest-frame 1500 A (M ) using the fol- UV;i divided by the number of input mock galaxies, as a func- lowing formula tion of input absolute magnitude and redshift. Figure 7 shows the results of our simulation. We de ne the 50%- M 0 = m 5 log(d =10pc) + 2:5 log(1 + z); (10) UV;i i L completeness limit as the absolute magnitude where the area under the curve falls to less than 50% the areas un- where d is the luminosity distance in pc. The second der the average of the M 0 = 25 to M 0 = 28 and third terms of the right side are the distance mod- UV;i UV;i curves, which we nd to be at M = 22 (m = 24 UV;i ulus. for z = 4). Our UV luminosity function is shown as red diamonds in Figure 8. We note that we do not include the lumi- 4. RESULTS nosity function data points in bins M 0 -22 in our UV;i analysis as these bins are below our 50%-completeness 4.1. The Rest-Frame UV z=4 Luminosity Function limit as discussed in Section 3.4, where the complete- We utilize the e ective volume method to correct for ness corrections are unreliable due to the low S/N of incompleteness in deriving our luminosity function. The our data. This is con rmed by comparing to the HST e ective volume (V ) can be estimated as eff CANDELS results in these same magnitude bins, which are more reliable due to their higher S/N. In Figure 8 dV V (M 0 ) = C (M 0; z)dz; (9) eff i UV;i we also include the UV luminosity function of z = 4 dz UV-selected galaxies from deeper Hubble imagining by dV where is the co-moving volume element, which de- Finkelstein et al. (2015) as green squares. In the bin dz 0 0 pends on the adopted cosmology, and C (M ; z) is the where we overlap with this dataset (M = -22.5) UV;i UV;i completeness as calculated in Section 3.4. The integral both results are consistent, however, if the luminosity was evaluated over z = 3 5. function derived from Hubble imaging is extrapolated to To calculate the luminosity function, we convert the brighter magnitudes, it would fall o more steeply than apparent i-band AB magnitudes (m 0 ) to the absolute our luminosity function. While our luminosity function ′ ′ g r [mag] i [3. 6] [mag] 14 Stevans et al. P(z) only UV, i 1.0 1.0 P(z) + i − [3.6] -28.0 0 0 0 P(z) + i − [3.6] + r − i -27.0 50%-Completeness Limit 0.8 -26.0 0.8 -25.0 -24.0 0.6 0.6 -23.0 -22.5 -22.0 0.4 0.4 -21.5 0.2 0.2 0.0 0.0 3.0 3.5 4.0 4.5 5.0 18 19 20 21 22 23 24 25 Input z Apparent Magnitude (mag) Figure 6. The results of our contamination simulations Figure 7. The results of our completeness simulations, estimating the contamination fraction using dimmed real showing the fraction of simulated sources recovered as a func- sources, showing the contamination fraction, F , as a func- tion of input redshift. Each colored line is for the correspond- 0 0 tion of m (solid and dashed colored lines) and m (dotted ing 0.5 magnitude bin according to the legend. We de ne i r colored lines). The color of the lines represent the selection the 50%-completeness limit as the absolute magnitude curve criteria applied before F was calculated: blue for F after with an integral value of less than 50% the area under the only the z PDF selection cuts, purple for F after the average of the M = 25 to M = 28 curves. We nd phot 1450 1450 z PDF selection and the i [3:6] color cuts, and orange the 50%-completeness to be M = 22:0. phot 1450 for F after the z PDF selection, the i [3:6] color, and the phot 0 0 r i color cuts. The error bars indicate the standard devia- Hubble imaging and the SDSS AGN number densities, tion of the mean from subdividing our 4.8 million simulated our data can potentially provide a robust measurement sources into 20 subsamples. The nal contamination frac- of the bright-end slope of the star-forming galaxy lu- tion was found to be between 0-20% generally increasing as minosity function and the faint-end slope of the AGN a function of magnitude and not exceeding 25% in any mag- nitude bin brighter than our 50%-completeness limit m 0 luminosity function, which we explore in the following 23.5 (black solid line). section. 4.2. Fitting the Luminosity Function 0 0 declines from M = 22 to M = 24, it attens UV;i UV;i out at brighter magnitudes before turning over again. With our luminosity function in agreement with the As our sample includes all galaxies which exhibit a faint end of the AGN luminosity function from SDSS Lyman break, we expect that this attening is due to DR7 and the bright end of the star-forming galaxy lu- the increasing importance of AGN at these magnitudes. minosity function from CANDELS, we attempt to si- The large volume surveyed by SDSS data has led to multaneously t empirically motivated functions to both components. For the AGN component we use a DPL the selection and spectroscopic follow-up of AGNs at many redshifts. SDSS AGN studies have found the AGN function motivated by the AGN UV luminosity func- tion work on large, homogeneous quasar samples (e.g., UV luminosity function to exhibit a double power law (DPL) shape (e.g., Richards et al. 2006). In Figure 8 we Boyle et al. 2000; Richards et al. 2006; Croom et al. show as red "x"s the z = 4 AGN UV number densities 2009; Hopkins et al. 2007, and references therein). The derived from the SDSS DR7 catalog (Schneider et al. function form of a DPL follows 2010) by Akiyama et al. (2018), who select AGNs to M > 28:9. We can see that at the magnitudes UV;i (M ) = ; (11) where we overlap 27 . M 0 . 26, the agreement UV;i 0:4( +1)(MM ) 0:4( +1)(MM ) 10 + 10 is excellent with our data. The only di erence is at where  is the overall normalization, M is the char- M = 28, where our survey detects two quasars when acteristic magnitude, is the faint-end slope, and is the AGN luminosity function by Akiyama et al. (2018) the bright-end slope. would predict less than one in our volume. We attribute this di erence to cosmic variance. By combining our For the star-forming galaxy UV luminosity function we consider sepraratley both a Schechter function and a data with the star-forming galaxy number densities from DPL, as well as including magni cation via gravitational Contamination Fraction, F Completeness Galaxy and AGN Luminosity Functions at z = 4 with SHELA 15 -2 3.0 -3 2.0 1.0 0.0 -4 27 26 25 24 23 UV, i -5 -6 DPL+DPL Fit Galaxy DPL Component -7 AGN DPL Component DPL+Sch Fit Galaxy Schechter Component -8 AGN DPL Component This Work This Work -9 (<50% complete) Finkelstein+15 SDSS DR7 (Akiyama+18) -10 Ono+18 28 26 24 22 20 18 UV, i Figure 8. The rest-frame UV z = 4 luminosity function of star-forming galaxies and AGNs from the SHELA Field shown as red diamonds with Poisson statistic error bars. The open red diamonds are the luminosity function points in bins below our 50%-completeness limit as discussed in Section 3.4, where the completeness corrections are unreliable due to the low S/N of our data. We constrain the form of the luminosity function by including fainter galaxies from Hubble elds (Finkelstein et al. 2015; open black squares) and brighter AGNs from SDSS DR7 (Akiyama et al. 2018; black \x"s). For comparison, we overplot as gray circles the g -band dropout luminosity function from the > 100 sq. deg. HSC SSP by Ono et al. (2018), which shows lower number densities and larger error bars in the regime (M < 23:5 mag) where AGN likely dominate (see Section 1). UV;i Our measured luminosity function is consistent with these works where they overlap. Our two best- tting functional forms are shown, as discussed in Section 4.2. The t with the smallest  is the sum of two DPL functions (DPL+DPL Fit; red solid line), one for the AGN component (red dash-dotted line) and one for the galaxy component (red dashed line). The second best t (DPL+Sch Fit; blue densely dotted line) is comprised of a DPL function for the AGN component (blue dotted line) and a Schechter function for the galaxy component (blue dash-dot-dotted line). The absolute value of the residuals of the two ts are shown in the inset plot for a subset of the data in units of the uncertainty in each bin. The data favors the DPL+DPL Fit suggesting the M 0 = 23:5 mag bin is dominated by star-forming galaxies, though this is dominated by the observed UV;i number densities in just a few bins, thus it is dicult to rule out a Schechter form for the galaxy component. −1 −3 Φ [Number mag Mpc ] Absolute Residual (σ) 16 Stevans et al. lensing with both functions. The Schechter (1976) func- We account for Eddington bias in our tting routine. tion has been found to t the star-forming galaxy UV Rather than directly comparing the observations to a luminosity function well across all redshifts (e.g., Steidel given model, we forward model the e ects of Eddington et al. 1999; Bouwens et al. 2007; Finkelstein et al. 2015). bias into the luminosity function model, and compare The Schechter function is described as this \convolved" model to our observations. We do this by, for each set of luminosity function parameters, re- 0:4 ln(10) (M ) = ; (12) 0:4(MM ) alizing a mock sample of galaxies for that given func- 0:4( +1)(MM ) 10 10 e tion, where each galaxy has a magnitude according to where  is the overall normalization, M is the char- the given luminosity function distribution. The magni- acteristic magnitude, and is the faint-end slope. tude of each simulated object is perturbed by an amount We consider the e ects of gravitational lensing on the drawn from a normal distribution centered on zero with shape of the star-forming galaxy UV luminosity func- a width equal to the real sample median uncertainty in tion. Gravitational lensing can distort the shape and the corresponding magnitude bin. After perturbing, we magni cation of distant sources as the paths of photons then re-bin the simulated luminosity distribution and from the source get slightly perturbed into the line of this binned luminosity function is used to calculate the sight of the observer. Lima et al. (2010) showed that for that MCMC step. This is repeated for each step. this magni cation can contribute to a bright-end ex- We ran our MCMC tting algorithm twice. In both cess where the slope of the intrinsic luminosity function runs, we simultaneously t a DPL function to the AGN is suciently steep. A magni cation distribution for a component of the luminosity function while tting the given source redshift must be estimated by tracing rays galaxy component. In the rst run, we t a Schechter through a series of lens planes derived from simulations function to the galaxy luminosity function, while in the such as the Millennium Simulation as done by Hilbert second we t a second DPL function to the galaxy lu- et al. (2007). van der Burg et al. (2010) showed that minosity function. We burn each chain for 210 steps, a Schechter function corrected for magni cation can t which allows the chain to reach convergence for all free the bright-end of the luminosity function at z = 3 better parameters, veri ed by examining the parameter distri- than a Schechter t alone. They inspect the sources that butions in independent groups of 10 steps, which cease make up the excess and nd nearby massive foreground to evolve much past 10 steps. We then measure the galaxies or groups of galaxies that could act as lenses. posterior for each parameter from the nal 510 steps. We incorporate the e ects of gravitational lensing in our Our ducial values for each parameter are then derived tting by creating a lensed Schechter function parame- as the median and 68% central con dence range from terization following the method of Ono et al. (2018) who the posterior distributions. The best tting parameters adapts the method of Wyithe et al. (2011). We also for the two ts and their corresponding  values are produce a lensed DPL function. After performing our listed in Table 4. We over plot the ts to the luminosity simultaneous tting method, which we describe in the function data in Figure 8. To calculate the  for our following paragraph, we found there to be no di erence ts, we must compare the observed data to the lumi- in the best tting parameters of the ts including and nosity function ts after forward modeling the e ects of excluding the e ects of lensing. This is consistent with Eddington bias into our luminosity function ts just as the work of Ono et al. (2018) who found that taking we did during the tting process. The absolute residuals into account the e ects of lensing improves the galaxy of our \convolved" ts and the observed data are shown luminosity function t at z > 4 and not at z = 4 where in the insert in Figure 8. a DPL t is preferred. Therefore we do not consider the These tting results show that our data prefers the lensed parameterizations further. function that is a sum of two DPL functions, one for We employ a Markov Chain Monte Carlo (MCMC) the AGN component and one for the galaxy compo- method to de ne the posteriors on our luminosity func- nent, henceforth the DPL+DPL Fit. The DPL+DPL tion parameterizations. We do this using an IDL imple- Fit has a  = 42 over the entire range considered, mentation of the ane-invariant sampler (Goodman & 29 < M 0 < 17. The t that is a sum of a UV;i Weare 2010) to sample the posterior, which is similar DPL AGN component and a Schechter galaxy compo- in production to the emcee package (Foreman-Mackey nent, henceforth the DPL+Schetcher Fit, has a  = 71. et al. 2013). Each of the 500 walkers was initialized We investigate which t is preferred by the data using by choosing a starting position with parameters deter- the Bayesian information criterion for the two ts (Lid- mined by-eye to exhibit a good t, perturbed according dle 2007). The di erence in the Bayesian information to a normal distribution. We do not assume a prior for criterion value for our two ts can be de ned as: any of our free parameters. Galaxy and AGN Luminosity Functions at z = 4 with SHELA 17 Table 4. Fit Parameters for Luminosity Functions AGN Component Galaxy Component Fit Name Log  M Log  M (mag) [Faint End] [Bright End] (mag) [Faint End] [Bright End] +0:21 +0:4 +0:30 +0:21 +0:09 +0:16 DPL+DPL -7.32 -26.5 -1.49 -3.65 -3.12 -20.8 -1.710.08 -3.800.10 42 0:18 0:3 0:21 0:24 0:10 0:15 +0:58 +1:1 +0:18 +0:68 DPL+Sch -7.48 -26.7 -2.08 -3.66 -3.250.06 -21.30.06 -1.810.05  71 0:34 0:4 0:11 0:34 3 1 Note| 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). 1.0 2 2 SDSS Sample BIC =   + (k k ) ln N; (13) 2 1 2 1 Selected 0.8 2 2 2 where  is the  for the DPL+Schechter Fit,  is 2 2 0.6 the  for the DPL+DPL Fit, k and k are the num- 10 2 1 0.4 ber of t parameters in the DPL+Schechter Fit and the DPL+DPL Fit, respectively, and N is the number of 0.2 data points used during tting. We nd a BIC = 26 0 0.0 which suggests the DPL+DPL Fit is very strongly pre- 3.5 4.0 4.5 5.0 3.0 3.5 4.0 4.5 5.0 z z spec spec ferred over the DPL+Schechter Fit as BIC exceeds a signi cance value of 10 (Liddle 2007). Upon compar- ing the the ts' residuals we see the data points in the Figure 9. Left: The SDSS spectroscopic redshift distribu- tion of AGNs in SHELA (blue) overplotted with the same magnitude bins at M 0 = 24 mag and M 0 = 23 UV;i UV;i distribution for the sub-sample selected by our selection pro- mag are driving the preference for the DPL+DPL Fit. If cedure (red). Right: The di erential completeness fraction we exclude these two bins, the DPL+DPL Fit  = 41, of AGNs with SDSS spectroscopic redshifts (blue) and the the DPL+Schechter Fit  = 58, and the BIC = 14 expected di erential completeness fraction for objects with which indicates the DPL+DPL Fit is still strongly sta- a comparable M (green). Our completeness to AGNs is UV;i tistically preferred to the DPL+Schechter Fit. However, similar to what we expect from our simulations of galaxies given that the di erence between the ts is driven by except at z = 3:5 where we are less complete to AGNs due to 0 0 signi cant u and g ux driving the redshift probability dis- just a few data points, we do not believe we can rmly tribution functions to peak at redshifts below our selection rule out a Schechter form for the star-forming galaxy window (see Section 4.3 for details). component. 4.3. A Sample of Spectroscopically Con rmed AGNs Each match required a separation of less than 0.4" and the SDSS spectra to be un agged (i.e., ZWARNING = Given our method of tting the luminosity function 0). We found zero spectra with an SDSS classi cation of with a component for AGNs, we explored the SDSS Galaxy and 53 classi ed as AGN with spectroscopic red- spectral catalog for spectroscopically con rmed AGNs shifts (z ) greater than 3.2. The distribution of SDSS spec in SHELA and considered the e ectiveness of our se- z for this sample is shown in the left panel of Figure 9 spec lection procedure at recovering AGNs. This is crucial in blue. Of the 32 AGNs with 3:4 < z < 4:6, 23 were spec as our photometric redshift selection process did not in- selected by our z selection process, with seven of the phot clude templates dominated by AGN emission, though nine missed AGNs having 3:4 < z < 3:5. The z spec spec the strong Lyman break inherent in these sources should distribution for the AGN subsample selected by our se- still yield an accurate redshift. To con rm this, we lection procedure is also shown in the left panel of Figure queried the spectral catalog from SDSS (DR13; Albareti 9 in red. We used these samples to compute our di eren- et al. 2017) using the SDSS CasJobs website and cross- tial completeness of AGNs and compared it to our sim- matched the results with our entire DECam catalog. ulated completeness in the right panel of Figure 9. We found our completeness of spectroscopically con rmed http://skyserver.sdss.org/casjobs/ AGNs is consistent with our simulated completeness ex- Number Completeness 18 Stevans et al. cept in the 3:4 < z < 3:6 bin where we recovered only data, the contamination rate was signi cantly higher, two of the nine spectroscopically con rmed AGNs when although we acknowledge there are multiple di erences our simulations predicted we should recover 78. This in the selection techniques between these studies. di erence could be due to small number statistics. We At the bright end, the luminosity function from investigated why the seven spectroscopically con rmed Ono et al. (2018) extends to magnitudes brighter than AGNs in the 3:4 < z < 3:6 bin were not selected by our M 0 < 22:5 and falls between the bright end of UV;i method and found that these sources had photometry, our DPL galaxy component and our Schechter galaxy 0 0 particularly the u and g bands, preferred by galax- component. Ono et al. (2018) nd their luminosity ies templates at lower redshift in our z - tting code function is best t by a DPL with a steeper bright-end phot 0 0 EAZY. We attribute bright u and g uxes to the larger slope ( = 4:33) than our galaxy DPL component far-UV continuum levels of AGNs or strong Lyman- ( = 3:80). We also include the z = 4 galaxy lumi- emission as compared to non-AGN galaxies, which nosity function from the 4 deg ALHAMBRA survey by would weaken the Lyman- break in these sources. This Viironen et al. (2017) who used a z PDF analysis phot could imply signi cant leaking Lyman-continuum radia- to create a luminosity function marginalizing over both tion from these AGNs, which has implications on reion- redshift and magnitude uncertainties. Viironen et al. ization (e.g., Smith et al. 2016). The reliability of our (2017) nd volume densities at the bright-end larger selection procedure to recover AGNs across the major- than existing luminosity functions. In fact, their lu- ity of our redshift range of interest and the fact that our minosity function follows a Schechter function with a luminosity function is consistent with the AGN luminos- normalization o set of 0.5 dex above the Schechter ity function from SDSS suggests our incompleteness to form from Finkelstein et al. (2015). The cause of this AGNs is not substantial. di erence is unclear. 5. DISCUSSION 5.1.1. Comparison to Semi-Analytic Models 5.1. Comparison to z=4 Galaxy Studies Figure 11 shows the predicted z = 4 UV luminosity We compare our derived ts to the star-forming galaxy functions from Yung et al. (2018). These predictions luminosity function to measurements from the litera- come from a set of semi-analytic models (SAMs), which ture in Figure 10. At magnitudes fainter than M 0 > contain the same physical processes as the models pre- UV;i 22, the DPL galaxy component of the data-preferred sented in Somerville et al. (2015), but have been up- DPL+DPL Fit is very similar to the Schechter compo- dated and recalibrated to the Planck 2016 Cosmological nent of the DPL+Schechter Fit. These results are in parameters. We note that while these models include strong agreement with the luminosity function from the black hole accretion and the e ects of AGN feedback, CFHT Deep Legacy Survey by van der Burg et al. (2010) we do not examine the contribution to the UV luminos- at all magnitudes where they overlap. They are also ity function from black hole accretion here, and focus in strong agreement with the luminosity function from instead on the contribution from star formation. Their Hubble legacy survey data by Finkelstein et al. (2015) ducial model with dust (solid black line) assumes that which is included in our tting. The luminosity func- the molecular gas depletion time is shorter at higher gas tion from the Hubble legacy survey data by Bouwens densities (as motivated by observations), leading to an et al. (2015) and the luminosity function from Ono et al. e ective redshift dependence as high redshift galaxies (2018) using Subaru Hyper Suprime-Cam (HSC) data are more compact and have denser gas on average. On across 100 deg in the HSC Subaru strategic program a SFR surface density versus gas surface density plot, (SSP) are generally consistent with the results presented this model would show a steeper dependence of star for- here. The HSC SSP luminosity function and that from mation rate density on gas surface density than the clas- Bouwens et al. (2015) have volume densities of 2x sical Kennicutt-Schmidt relationship (e.g., Kennicutt & larger than the work by van der Burg et al. (2010) Evans 2012), above a critical gas surface density (for and Finkelstein et al. (2015) at magnitudes fainter than details see Somerville et al. 2015). We also show their M > 21. This factor is larger than the 10% un- model with a xed molecular gas depletion time, similar UV;i certainty expected due to cosmic variance in the Hubble to that used in Somerville et al. (2015), as seen in local elds, which cover 50x less area than the HSC SSP spiral galaxies (e.g., Bigiel et al. 2008; Leroy et al. 2008). (Finkelstein et al. 2015). We do not know the cause of The Yung et al. (2018) dust models assume that the this di erence, though we note that the HSC SSP selec- V-band dust optical depth is proportional to the cold gas tion was done with optical imaging only, and we found metallicity times the cold gas surface density. The UV in our study that without the addition of the IRAC attenuation is then computed using a Calzetti attenua- Galaxy and AGN Luminosity Functions at z = 4 with SHELA 19 -2 This Work -2 Finkelstein+15 -3 -3 -4 -4 -5 -5 Galaxy DPL Component Galaxy DPL Component -6 Galaxy Schechter Component -6 Galaxy Schechter Component Finkelstein+15 SAMs w/dust Bouwens+15 -7 -7 SAMs w/dust & fixed H t 10 2 dep 10 van der Burg+10 SAMs w/no dust Ono+18 -8 SAMs w/no dust & fixed H t -8 2 dep Viironen+17 10 10 24 23 22 21 20 19 18 17 25 24 23 22 21 20 19 18 M 0 M 0 UV, i UV, i Figure 10. The star-forming galaxies components of our Figure 11. The star-forming galaxy component of our ts ts to the total z = 4 rest-frame UV luminosity function to the total z = 4 rest-frame UV luminosity function com- compared to the data from other star-forming galaxy stud- pared to predictions from SAMs (Yung et al. 2018). Our ies See the legend insert for the list of compared works. Our measured total luminosity function (including star-forming DPL galaxy luminosity function is in agreement with lumi- galaxies and AGNs) is shown as red diamonds with Poisson nosity functions from the literature around the knee and to statistic error bars. The luminosity function from Finkel- fainter magnitudes, but has the shallowest bright-end slope stein et al. (2015) is shown as open black squares. At the ( = 3:80). faint end the luminosity functions from the SAMs with an evolving H gas depletion time (solid black line) and a xed H gas depletion time (dashed black line) have a higher nor- tion curve and a \slab" model (for details see Somerville malization than the observed luminosity function suggesting et al. 2015). The normalization of the dust optical depth stellar feedback in the models is too weak. At the bright (physically equivalent to the dust-to-metal ratio) is al- end the model with a xed molecular depletion timescale is lowed to vary as a function of redshift, and was adjusted robustly ruled out by either of our parameterized ts, sug- gesting that star formation scales with molecular gas surface to t the bright end of the luminosity from the previ- density, and thus is more ecient at z = 4 than today. ously published compilation of UV luminosity observa- tions from Finkelstein (2016). (It was not adjusted to t the new observations presented here). At the bright end, the Yung et al. (2018) model with At the faint end, the model predictions are insensitive dust is consistent with both of our ts to M > UV;i to the assumed star formation eciency, and mainly re- 22.5, lying closer to the Schechter t at brighter mag- ect the treatment of out ows driven by stars and su- nitudes. Interestingly, at these bright magnitudes both pernovae. The models have a higher normalization than the Schechter and DPL t rule out at high signi cance the observed luminosity function (which at these mag- the model with dust and xed molecular gas depletion nitudes comes from the CANDELS dataset), although time, indicating that star formation must scale non- the faint-end slope is similar. This could be caused by linearly with molecular gas surface density (or some re- stellar feedback in the models being too weak, as these lated quantity which evolves with redshift). This implies models were tuned to match the z = 0 stellar mass func- that star formation is more ecient at z = 4 than to- tion (Somerville et al. 2015). This suggests that stellar day. This is of course dependant on the dust model feedback is stronger/more ecient at z = 4 (i.e., mass in these simulations { if bright galaxies had no dust, loading rates are higher, or re-infall of ejected gas is then these model predictions (which include dust; mod- slower) than at z = 0 (also see White et al. 2015). This els without dust are shown for comparison) may not be was also seen by Song et al. (2016) when comparing the accurate. However, there is a variety of evidence that z = 4 stellar mass function to a number of similar mod- bright/massive UV-selected galaxies at these redshifts els { the observed stellar mass function also had a lower do contain non-negligible amounts of dust (e.g., Finkel- normalization at z = 4; as the stellar mass function stein et al. 2012; Bouwens et al. 2014), and there is a steepened from z = 4 to 8, this discrepancy weakened, not-insigni cant population of extremely massive dusty hence their conclusion of a weakening impact of feedback star-forming galaxies already in place (e.g., Casey et al. on the faint end with increasing redshift. 2014). Finally, Finkelstein et al. (2015) compared a sim- −1 −3 Φ [Number mag Mpc ] −1 −3 Φ [Number mag Mpc ] 20 Stevans et al. +0:18 ilar set of models to the CANDELS-only UV luminos- 2:08 ) predicts volume densities in agreement with 0:11 ity functions, nding that models with a di use dust Glikman et al. (2011) and Giallongo et al. (2015), who component only (e.g., no birth cloud) provided the best predict a signi cantly higher volume density of faint match to the data (albeit at fainter magnitudes than AGNs than the other studies. We note that while these considered here). studies have small numbers of AGN per magnitude bin, which may make the samples more sensitive to cosmic variance and false positives, Glikman et al. (2011) se- 5.2. Comparison to z = 4 AGN Studies lects candidates from a relatively large survey area of We compare our derived AGN luminosity function ts 3.76 deg and observes broad emission lines, indicative to measurements from the literature in Figure 12. At the of quasars, in every candidate spectra. Furthermore, bright end, our ts are consistent with previous studies Giallongo et al. (2015) requires candidates to have a where they overlap (27:5 < M < 25:5). The pre- UV;i signi cant X-ray detection in Chandra imaging. vious studies include a study by Richards et al. (2006) In summary, our observed z = 4 UV luminosity func- who used the z = 4 AGN SDSS DR3 sample and a tion is best t by the DPL+DPL Fit, but both the study by Akiyama et al. (2018) who selected z = 4 AGNs DPL+DPL and the DPL+Schechter ts show agreement from SDSS DR7. We convert the magnitudes M at i(z=2) with existing AGN luminosity function studies at the z = 3:75 from Richards et al. (2006) to M 0 at z = 3:8 UV;i bright end and around the knee. At the faint end the by adding an o set of 1.486 mag (Richards et al. 2006). DPL+DPL predicts smaller volume densities of AGNs At the faint end, we compare our ts to results from than other studies while the DPL+Schechter predicts studies that derive AGN luminosity functions from spec- among the largest volume densities at the faintest mag- troscopic observations of candidates selected via color nitudes. and size criteria (Glikman et al. 2011, Ikeda et al. 2011, and Niida et al. 2016). In addition, we include studies that rely on a z selection using deep multi wave- phot 5.3. Comparing Predictions of Our Two Fits - Is the length photometry in the COSMOS eld (Masters et al. UV luminosity function a double power law? 2012) and the CANDELS GOODS-S eld (Giallongo et al. 2015) where Giallongo et al. (2015) had the ad- Our two ts di er in two ways: 1.) The functional ditional criteria of requiring an an X-ray detection of form used to t the star-forming galaxy component is 17 2 1 F > 1:5 10 erg cm s in deep 0.5-2 keV Chan- a DPL in the DPL+DPL Fit and a Schechter in the dra imaging. The average redshift of these samples is DPL+Schechter Fit, and 2.) The component that ac- z = 4, slightly higher than our sample. For a consistent counts for the excess over an extrapolated Schechter at comparison with other works we scale the luminosity the bright end of existing star-forming galaxy luminos- functions of Glikman et al. (2011), Ikeda et al. (2011), ity functions. The DPL+Schechter accounts for the ex- Niida et al. (2016), and Masters et al. (2012) up by a fac- cess with a steeper faint end of the DPL AGN compo- 6:9 tor of 1.3 using the redshift evolution function (1+z) nent while the DPL+DPL Fit accounts for the major- from Richards et al. (2006). The sample used by (Gial- ity of the excess with the bright end of the DPL star- longo et al. 2015) had a redshift range of 4 < z < 4:5, so forming galaxy component. Thus the ts predict dra- we the scale the luminosity function by a factor of 1.8. matically di erent compositions of sources making up Finally, we consider the luminosity function of Akiyama the bin (M 0 = 23:5 mag) at the center of the ex- UV;i et al. (2018) derived from the Hyper Suprime-Cam Wide cess. In this bin the DPL+Schechter Fit predicts AGNs Survey optical photometry. outnumber galaxies by a 17:1 ratio (an AGN fraction The faint end of the AGN DPL component in our of 94%), while the DPL+DPL Fit predicts galaxies DPL+DPL Fit predicts volume densities at 26 < to outnumber AGNs by a 4.3:1 ratio (an AGN fraction M 0 < 23:5 about 0.3 dex lower than those found by of 18%). A simple experiment to test these predic- UV;i Richards et al. (2006), Ikeda et al. (2011), Masters et al. tions would be to use ground-based optical spectroscopy (2012), Niida et al. (2016), Akiyama et al. (2018). How- to follow-up a fraction of the 298 high-z candidates we ever, the luminosity function of Akiyama et al. (2018) nd in the M = 23:5 bin and count the fraction UV;i attens and falls towards our t at M 0  22. Our of spectra exhibiting broad emission lines and/or highly UV;i +0:30 faint-end slope ( = 1:49 ) is in agreement with ionized lines (e.g., N V, He II, C IV, Ne V) indicating ac- 0:21 the values found by these studies. The cause of the 0.3 cretion onto supermassive black holes. This experiment dex di erence is unclear. would also provide an independent measurement of our In the case of the AGN DPL component in our contamination fraction which we estimate via simula- DPL+Schechter Fit, the steeper faint-end slope ( = tions (Section 3.3). The measured fraction of AGNs can Galaxy and AGN Luminosity Functions at z = 4 with SHELA 21 -4 Richards+06 SDSS DR7 (Akiyama+18) HSC (Akiyama+18) -5 Glikman+11 -6 -7 -8 AGN DPL with galaxy DPL AGN DPL with galaxy Schechter Ikeda+11 -9 Niida+16 Masters+12 -10 Giallongo+15 28 26 24 22 20 UV, i Figure 12. The AGN components of our ts to the total z = 4 rest-frame UV luminosity function compared to the data from other AGN studies. See the legend insert for the list of compared works. The AGN luminosity function from the DPL+DPL Fit has number densities similar to existing luminosity function measurements at M < 25 and predicts relatively low number UV;i densities at fainter magnitudes. The DPL+Schechter Fit has an AGN faint-end slope that predicts number densities at the larger end of the range previously published. Table 5. UV Luminosity Densities provide strong empirical proof in favor of the DPL+DPL Fit or the DPL+Schechter Fit. AGN Component Galaxy Component One drawback of tting the total UV luminosity func- tion is the unknown contribution of composite objects. The total UV luminosity function likely contains a pop- Fit Name 1500 912 1500 912 ulation of composite objects with both AGN activity +0:7 +0:4 +6 DPL+DPL 2.0 1.2 184 0:4 0:3 7 and star-formation contributing to their observed UV +4 +2:6 +6 DPL+Sch 10 5.8 187 3 1:8 6 ux. This population may cause functions like a DPL 24 1 1 3 Note|All values are in units of 10 ergs s Hz Mpc . or Schechter function, which have been used to t AGN- only and star-forming-galaxy-only luminosity functions, to be poor ts to the total UV luminosity function. Spectroscopic follow-up as described in this section may aid in elucidating the frequency and impact of composite that of AGNs by a factor of 90. In the case of the objects. DPL+Schechter Fit where the galaxy component of the UV luminosity function is instead represented by a DPL 5.4. Rest-Frame UV Emissivities function, the galaxy  is greater than the AGN 1500 1500 by a factor of 19. In either case, galaxies are the dom- Here we calculate the rest-frame UV emissivities (also inant non-ionizing UV emitting population, though if known as speci c luminosity densities) and compare the AGNs do have a steeper faint-end slope (as would need output from AGNs and galaxies. We calculate the rest- to be the case if galaxies follow a Schechter form), then frame UV emissitvity of AGNs and galaxies by integrat- AGNs are non-negligible. ing each luminosity function from 30 < M 0 < 17. UV;i We convert the AGN  to a UV luminosity den- Results are shown in Table 5. In the case of DPL+DPL sity at 912 A ( ) and compare our results to other where the galaxy component of the UV luminosity func- studies by reproducing Figure 1 from Madau & Haardt tion is represented by a DPL function and the AGN (2015) in Figure 13. To convert  to  we as- 1500 912 component is represented by a DPL, galaxies produce a sume an AGN spectral shape of a DPL with a slope UV luminosity density at 1500 A ( ) greater than −1 −3 Φ [Number mag Mpc ] 22 Stevans et al. faint AGNs at z  4 Glikman et al. 2011; Giallongo et al. 2015. This would imply the AGN population could alone contribute enough hydrogen ionizing emission required to keep the universe ionized at z  4 and suggests AGNs may have played a signi cant role at early times. This scenario can be further constrained by the spectroscopic experiment proposed in the preceding sub-section. DPL+Sch Fit 6. CONCLUSIONS DPL+DPL Fit In this study, we measure the bright end of the rest- Madau & Haardt 15 frame UV luminosity function of z = 4 star-forming Hopkins+07 galaxies and the faint end of the rest-frame UV lumi- 0 1 2 3 4 5 6 7 nosity function of z = 4 AGNs. We use nine photo- 0 0 0 0 0 metric bands (u g r i z from DECam, J and K from VISTA, and 3:6 and 4:5 m from Spitzer /IRAC) cov- ering the wide area (18 deg ) SHELA Field to select Figure 13. The AGN hydrogen ionizing emissivity from 3,740 candidate z  4 galaxies via a photometric red- this work and others. The AGN emissivity predicted by the shift selection procedure. From simulations, we nd a DPL+DPL Fit and the DPL+Schechter Fit are represented relatively low contamination rate of interloping low-z as red triangles. The range they span is marked by the thick red line. The values from other works are inferred from Fig galaxies and galactic stars of 20% near our complete- 1 of Madau & Haardt (2015). The original sources of each ness limit (m  23) due in large part to the inclusion dataset are as follows: Schulze et al. (2009) (cyan pentagon), of IRAC photometry. Palanque-Delabrouille et al. (2013) (orange triangles), Bon- Our conclusions are as follows: giorno et al. (2007) (magenta circles), Masters et al. (2012) We combine our candidate sample with a sample of (black pentagons), Glikman et al. (2011) (blue square), and Giallongo et al. (2015) (green squares). The solid blue line bright AGNs from SDSS and fainter galaxies from is the functional form derived by Madau & Haardt (2015) deep Hubble imaging (including the HUDF and to coincide with the plotted observation from the literature. CANDELS) to produce a rest-frame UV luminos- The dotted green line shows the LyC AGN emissivity from ity function that spans the range -29< M 0 <- UV;i Hopkins et al. (2007). The  from our DPL+DPL Fit 17. This range contains both AGNs and star- is below the line by Hopkins et al. (2007) indicating AGNs forming galaxies several magnitudes above and contribute only a small fraction of the total ionizing back- below their respective characteristic luminosities, ground at z = 4, while the  from the DPL+Schechter t suggests AGNs would contribute signi cantly to the total thus we implement a tting procedure that simul- ionizing background. taneously ts the AGN luminosity function and the star-forming galaxy luminosity function with of = 1:41 (Shull et al. 2012) between 1000 A and independent functions. This simpli es the source selection process by not requiring a step for classi- 1500 A and a slope of = 0:83 (Stevans et al. 2014) between 912 A and 1000 A . We assume an AGN ioniz- fying objects as either an AGN or galaxy, which is commonly done with morphological criteria. We ing escape fraction of unity as is found for most bright nd the data is best t by our DPL+DPL Fit AGNs (Giallongo et al. 2015). Results are shown in Ta- which is a combination of a DPL function for ble 5. If we directly follow the work of Giallongo et al. the AGN component and a DPL function for the (2015) and assume a spectral shape with = 1:57 galaxy component. The DPL+DPL t is preferred (Telfer et al. 2002) between 1200 A and 1500 A and over the DPL+Shechter Fit, which is a combina- a slope of = 0:44 (Vanden Berk et al. 2001) be- tween 912 A and 1200 A , we nd  values that are tion of a DPL function for the AGN component and a Schechter function for the galaxy compo- 10% smaller. As shown in Figure 13, our two ts pre- dict values of  straddle existing values from observa- nent, and this excess over Schechter cannot be explained by the e ects of gravitational lensing. tions. The  predicted by our preferred DPL+DPL Fit (upward-pointed red triangle) falls below the line by We note that we cannot signi cantly rule out a Schechter form. Hopkins et al. (2007) and is too small to keep the IGM ionized at z  4. On the other hand, the DPL+Sch We compare our measured luminosity functions Fit predicts a  (downward-pointed red triangle) that to the literature and nd our DPL galaxy lumi- falls near the points from studies that found numerous nosity function is in agreement with luminosity −1 −1 −3 ρ [erg s Hz Mpc ] 912 Galaxy and AGN Luminosity Functions at z = 4 with SHELA 23 functions from the literature around the knee and Future work is needed to con rm the shape of the to fainter magnitudes while having the shallow- star-forming galaxy and AGN luminosity functions, es- est bright-end slope. The AGN luminosity func- pecially where they intersect. We discuss a simple ex- tion from the DPL+DPL Fit has number densities periment to measure the relative number densities of similar to existing luminosity functions at magni- AGNs and galaxy at luminosities where the respective tudes up to M = 25:5 mag while under pre- luminosity functions intersect. Spectroscopic follow-up UV;i dicting number densities by  0:3 dex at fainter of a sample of our z = 4 candidates in the M -23.5 UV;i magnitudes. The DPL+Schechter Fit has an AGN bin where our two ts predict di erent AGNs to galaxies faint-end slope that is among the steepest values ratios is underway. Imaging from space-based telescopes published. The shape of the galaxy bright end is such as HST or JWST would facilitate a robust mor- consistent with model predictions where star for- phological classi cation of our bright candidates. Other mation is more ecient at higher redshift due to possibilities for distinguishing AGNs including taking increased gas densities. deeper X-ray imaging in the eld and using JWST to measure mid-IR SEDs (Kirkpatrick et al. 2012). The authors acknowledge Raquel Martinez, Karl Geb- We measure  by integrating the rest-frame hardt, and Eric J. Gawiser for the interesting discussions UV luminosity function ts and nd that galax- we had and their suggestions which improved the quality ies dominate the production of non-ionizing ux of this research. We thank Neal J. Evans and William P. at z = 4 for both possible ts. Speci cally, galax- Bowman for useful comments on the draft of this paper. ies produce a factor of 90 more non-ionizing UV M. L. S. and S. L. F. acknowledge support from the Uni- output than AGNs according to the DPL+DPL versity of Texas at Austin, the NASA Astrophysics and Fit, while the DPL+Schechter Fit predicts galax- Data Analysis Program through grant NNX16AN46G, ies produce a factor of 19 more non-ionizing UV and the National Science Foundation AAG Award AST- output than AGNs. 1614798. The work of C. P. And L. K. is supported by the National Science Foundation through grants AST 413317, and 1614668. The Institute for Gravitation and We convert the AGN  to  and nd AGNs 1500 912 the Cosmos is supported by the Eberly College of Sci- do not produce enough ionizing radiation to keep ence and the Oce of the Senior Vice President for Re- the universe ionized at z = 4 by themselves if the search at the Pennsylvania State University. S. J., S. AGN is truly represent by the DPL component S. and J. F. acknowledge support from the University of in our DPL+DPL Fit. This suggests AGNs are Texas at Austin and NSF grants NSF AST-1614798 and not the dominant contributor to cosmic reioniza- NSF AST-1413652. R. S. S. and A. Y. thank the Downs- tion at earlier times. On the other hand, if the brough family for their generous support, and gratefully DPL+Schechter Fit is true, AGNs could alone pro- acknowledge funding from the Simons Foundation. duce the  needed to maintain the ionized state of the universe at z = 4 and perhaps at earlier times. Facilities: CTIO,SST(IRAC) REFERENCES Akiyama, M., He, W., Ikeda, H., et al. 2018, PASJ, 70, S34, Bertin, E., Mellier, Y., Radovich, M., et al. 2002, in Astronomical Society of the Paci c Conference Series, doi: 10.1093/pasj/psx091 Vol. 281, Astronomical Data Analysis Software and Albareti, F. D., Allende Prieto, C., Almeida, A., et al. 2017, Systems XI, ed. D. A. Bohlender, D. Durand, & T. H. Handley, 228 ApJS, 233, 25, doi: 10.3847/1538-4365/aa8992 Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846, Beers, T. C., Flynn, K., & Gebhardt, K. 1990, doi: 10.1088/0004-6256/136/6/2846 Astronomical Journal (ISSN 0004-6256), 100, 32 Bouwens, R. J., Illingworth, G. D., Franx, M., & Ford, H. 2007, The Astrophysical Journal, 670, 928 Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. doi: 10.1051/aas:1996164 2014, ApJ, 793, 115, doi: 10.1088/0004-637X/793/2/115 24 Stevans et al. |. 2015, ApJ, 803, 34, doi: 10.1088/0004-637X/803/1/34 Kirkpatrick, A., Pope, A., Alexander, D. M., et al. 2012, ApJ, 759, 139, doi: 10.1088/0004-637X/759/2/139 Boyle, B. J., Shanks, T., Croom, S. M., et al. 2000, Monthly Notices of the Royal Astronomical Society, 317, 1014 Kurucz, R. L. 1979, ApJS, 40, 1, doi: 10.1086/190589 Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, LaMassa, S. M., Urry, C. M., Cappelluti, N., et al. 2016, The Astrophysical Journal, 686, 1503 The Astrophysical Journal, 817, 172 Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, Lang, D., Hogg, D. W., & Mykytyn, D. 2016a, The Tractor: Astrophysical Journal, 345, 245 Probabilistic astronomical source detection and measurement, Astrophysics Source Code Library. Casey, C. M., Scoville, N. Z., Sanders, D. B., et al. 2014, ApJ, 796, 95, doi: 10.1088/0004-637X/796/2/95 http://ascl.net/1604.008 Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833, Lang, D., Hogg, D. W., & Schlegel, D. J. 2016b, AJ, 151, doi: 10.1088/0004-637X/712/2/833 36, doi: 10.3847/0004-6256/151/2/36 Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486, Le F evre, O., Cassata, P., Cucciati, O., et al. 2013, doi: 10.1088/0004-637X/699/1/486 Astronomy & Astrophysics, 559, A14 Lehmer, B. D., Xue, Y. Q., Brandt, W. N., et al. 2012, The Croom, S. M., Richards, G. T., Shanks, T., et al. 2009, Monthly Notices of the Royal Astronomical Society, 399, Astrophysical Journal, 752, 46 1755 Leroy, A. K., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2782, doi: 10.1088/0004-6256/136/6/2782 Cucciati, O., Tresse, L., Ilbert, O., et al. 2012, Astronomy & Astrophysics, 539, A31 Liddle, A. R. 2007, MNRAS, 377, L74, Eisenstein, D. J., Weinberg, D. H., Agol, E., et al. 2011, doi: 10.1111/j.1745-3933.2007.00306.x AJ, 142, 72, doi: 10.1088/0004-6256/142/3/72 Lima, M., Jain, B., & Devlin, M. 2010, MNRAS, 406, 2352, Finkelstein, S. L. 2016, PASA, 33, e037, doi: 10.1111/j.1365-2966.2010.16884.x doi: 10.1017/pasa.2016.26 Madau, P. 1995, ApJ, 441, 18, doi: 10.1086/175332 Finkelstein, S. L., Papovich, C., Salmon, B., et al. 2012, Madau, P., & Haardt, F. 2015, The Astrophysical Journal ApJ, 756, 164, doi: 10.1088/0004-637X/756/2/164 Letters, 813, L8 Finkelstein, S. L., Song, M., Behroozi, P., et al. 2015, Marchesini, D., van Dokkum, P. G., F orster Schreiber, ArXiv e-prints. https://arxiv.org/abs/1504.00005 N. M., et al. 2009, ApJ, 701, 1765, Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, doi: 10.1088/0004-637X/701/2/1765 J. 2013, PASP, 125, 306, doi: 10.1086/670067 Masters, D., Capak, P., Salvato, M., et al. 2012, The Astrophysical Journal, 755, 169 Gawiser, E., van Dokkum, P. G., Herrera, D., et al. 2006, ApJS, 162, 1, doi: 10.1086/497644 McGreer, I. D., Fan, X., Jiang, L., & Cai, Z. 2017, arXiv.org Geach, J. E., Lin, Y. T., Makler, M., et al. 2017, The Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, Astrophysical Journal Supplement Series, 231, 7 777, 18, doi: 10.1088/0004-637X/777/1/18 Giallongo, E., Grazian, A., Fiore, F., et al. 2015, Niida, M., Nagao, T., Ikeda, H., et al. 2016, The Astronomy & Astrophysics, 578, A83 Astrophysical Journal, 832, 208 Glikman, E., Djorgovski, S. G., Stern, D., et al. 2011, The Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713, Astrophysical Journal Letters, 728, L26 doi: 10.1086/160817 Goodman, J., & Weare, J. 2010, Communications in Ono, Y., Ouchi, M., Harikane, Y., et al. 2018, PASJ, 70, Applied Mathematics and Computational Science, Vol. 5, S10, doi: 10.1093/pasj/psx103 No. 1, p. 65-80, 2010, 5, 65, Papovich, C., Shipley, H. V., Mehrtens, N., et al. 2016, doi: 10.2140/camcos.2010.5.65 ApJS, 224, 28, doi: 10.3847/0067-0049/224/2/28 Hilbert, S., White, S. D. M., Hartlap, J., & Schneider, P. Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, 2007, MNRAS, 382, 121, Journal of Machine Learning Research, 12, 2825 doi: 10.1111/j.1365-2966.2007.12391.x Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. Hopkins, P. F., Richards, G. T., & Hernquist, L. 2007, The 2014, A&A, 571, A16, doi: 10.1051/0004-6361/201321591 Astrophysical Journal, 654, 731 Quadri, R., Marchesini, D., van Dokkum, P., et al. 2006, Ikeda, H., Nagao, T., Matsuoka, K., et al. 2011, The arXiv.org, 1103 Astrophysical Journal Letters, 728, L25 Richards, G. T., Strauss, M. A., Fan, X., et al. 2006, Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531, arXiv.org, 2766 doi: 10.1146/annurev-astro-081811-125610 Schechter, P. 1976, Astrophysical Journal, 203, 297 Galaxy and AGN Luminosity Functions at z = 4 with SHELA 25 Schla y, E. F., & Finkbeiner, D. P. 2011, The Szalay, A. S., Connolly, A. J., & Szokoly, G. P. 1999, The Astronomical Journal, 117, 68 Astrophysical Journal, 737, 103 Telfer, R. C., Zheng, W., Kriss, G. A., & Davidsen, A. F. Schneider, D. P., Richards, G. T., Hall, P. B., et al. 2010, 2002, The Astrophysical Journal, 565, 773 The Astronomical Journal, 139, 2360 van der Burg, R. F. J., Hildebrandt, H., & Erben, T. 2010, Shull, J. M., Stevans, M., & Danforth, C. W. 2012, The A&A, 523, A74, doi: 10.1051/0004-6361/200913812 Astrophysical Journal, 752, 162 Vanden Berk, D. E., Richards, G. T., Bauer, A., et al. 2001, Smith, B. M., Windhorst, R. A., Jansen, R. A., et al. 2016, The Astronomical Journal, 122, 549 ArXiv e-prints. https://arxiv.org/abs/1602.01555 Viironen, K., L opez-Sanjuan, C., Hern andez-Monteagudo, C., et al. 2017, arXiv.org, arXiv:1712.01028 Somerville, R. S., Popping, G., & Trager, S. C. 2015, White, C. E., Somerville, R. S., & Ferguson, H. C. 2015, MNRAS, 453, 4337, doi: 10.1093/mnras/stv1877 ApJ, 799, 201, doi: 10.1088/0004-637X/799/2/201 Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, Wyithe, J. S. B., Yan, H., Windhorst, R. A., & Mao, S. ApJ, 825, 5, doi: 10.3847/0004-637X/825/1/5 2011, Nature, 469, 181, doi: 10.1038/nature09619 Steidel, C. C., Adelberger, K. L., Giavalisco, M., Dickinson, Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., M., & Pettini, M. 1999, ApJ, 519, 1, doi: 10.1086/307363 Popping, G., & Dav e, R. 2018, ArXiv e-prints. Steidel, C. C., Giavalisco, M., Dickinson, M., & Adelberger, https://arxiv.org/abs/1803.09761 Zacharias, N., Monet, D. G., Levine, S. E., et al. 2004, in K. L. 1996, AJ, 112, 352, doi: 10.1086/118019 Bulletin of the American Astronomical Society, Vol. 36, Stevans, M. L., Shull, J. M., Danforth, C. W., & Tilton, American Astronomical Society Meeting Abstracts, 1418 E. M. 2014, The Astrophysical Journal, 794, 75 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Astrophysics arXiv (Cornell University)

Bridging Star-Forming Galaxy and AGN Ultraviolet Luminosity Functions at $z=4$ with the SHELA Wide-Field Survey

Loading next page...
 
/lp/arxiv-cornell-university/bridging-star-forming-galaxy-and-agn-ultraviolet-luminosity-functions-no3uFT1Std

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

ISSN
1538-4357
eISSN
ARCH-3330
DOI
10.3847/1538-4357/aacbd7
Publisher site
See Article on Publisher Site

Abstract

Draft version June 15, 2018 Typeset using L T X twocolumn style in AASTeX62 Bridging Star-Forming Galaxy and AGN Ultraviolet Luminosity Functions at z = 4 with the SHELA Wide-Field Survey 1 1 1 2, 3, 4 2, 3 Matthew L. Stevans, Steven L. Finkelstein, Isak Wold, Lalitwadee Kawinwanichakij, Casey Papovich, 1 5, 6 1 5, 6 1 Sydney Sherman, Robin Ciardullo, Jonathan Florez, Caryl Gronwall, Shardha Jogee, 7, 8 7 Rachel S. Somerville, and L. Y. Aaron Yung Department of Astronomy, University of Texas at Austin, Austin, TX 78705 Department of Physics and Astronomy, Texas A&M University, College Station, TX, 77843-4242 USA George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M University, College Station, TX, 77843-4242 USA LSSTC Data Science Fellow Department of Astronomy & Astrophysics, The Pennsylvania State University, University Park, PA 16802 Institute for Gravitation and the Cosmos, The Pennsylvania State University, University Park, PA 16802 Department of Physics & Astronomy, Rutgers University, Piscataway, New Jersey 08854 Center for Computational Astrophysics, Flatiron Institute, New York, NY, 10010 (Received March 6, 2018; Revised June 8, 2018; Accepted June 9, 2018) Submitted to ApJ ABSTRACT We present a joint analysis of the rest-frame ultraviolet (UV) luminosity functions of continuum- selected star-forming galaxies and galaxies dominated by active galactic nuclei (AGNs) at z  4. These 3,740 z  4 galaxies are selected from broad-band imaging in nine photometric bands over 18 deg in the Spitzer/HETDEX Exploratory Large Area Survey (SHELA) eld. The large area and moderate depth of our survey provide a unique view of the intersection between the bright end of the galaxy UV luminosity function (M < 22) and the faint end of the AGN UV luminosity function. We AB do not separate AGN-dominated galaxies from star-formation-dominated galaxies, but rather t both luminosity functions simultaneously. These functions are best t with a double power-law (DPL) for both the galaxy and AGN components, where the galaxy bright-end slope has a power-law index of +0:30 3:80  0:10, and the corresponding AGN faint-end slope is = 1:49 . We cannot rule AGN 0:21 out a Schechter-like exponential decline for the galaxy UV luminosity function, and in this scenario +0:18 the AGN luminosity function has a steeper faint-end slope of 2:08 . Comparison of our galaxy 0:11 luminosity function results with a representative cosmological model of galaxy formation suggests that the molecular gas depletion time must be shorter, implying that star formation is more ecient in bright galaxies at z = 4 than at the present day. If the galaxy luminosity function does indeed have a power-law shape at the bright end, the implied ionizing emissivity from AGNs is not inconsistent with previous observations. However, if the underlying galaxy distribution is Schechter, it implies a signi cantly higher ionizing emissivity from AGNs at this epoch. Keywords: galaxies: high-redshift, galaxies: active{quasars, galaxies: luminosity function 1. INTRODUCTION omy. With the number of massive galaxies increasing from z  4 2 (Marchesini et al. 2009; Muzzin et al. Explaining how galaxies grow and evolve over cosmic 2013) and the positive relation between stellar mass and time is one of the main goals of extragalactic astron- star formation rate, by studying the properties of galax- ies with the highest star formation rates at z  4 we Corresponding author: Matthew L. Stevans can glean how the most massive galaxies built up their [email protected] stellar mass. The use of multi-wavelength photometry arXiv:1806.05187v1 [astro-ph.GA] 13 Jun 2018 2 Stevans et al. and the Lyman break technique has revolutionized the point sources. However, this method is less reliable study of galaxies in the z > 2 universe (e.g., Steidel near the photometric limit especially in ground-based et al. 1996). These tools are currently the most ecient imaging with poor seeing. For example, recent work by for selecting large samples of high-redshift star-forming Akiyama et al. (2018) has shown such a morphological galaxies for extensive study. A power tool for under- selection can select a sample of point sources with only standing the distribution of star-formation at high red- 40% completeness and 30% contamination at i = 24 mag shifts is the rest-frame UV luminosity function. This in photometry with median seeing conditions of 0: 7 and probes recent unobscured star-formation directly over 5- depths of i = 26:4 mags. the last 100 Myr and is, therefore, a fundamental tracer The shape of the AGN luminosity function is of in- of galaxy evolution. terest as well, as a steep faint end can result in a non- The shape of the star-forming galaxy UV luminos- negligible contribution of ionizing photons from AGNs ity function at z = 4 has been dicult to pin down to the total ionizing budget. Current uncertainties in the at the bright end. The characteristic luminosity of the literature at z  4 are large (Glikman et al. 2011, Mas- Schechter function, which is often used to describe the ters et al. 2012, Giallongo et al. 2015), thus AGNs have luminosity function in eld environments, ranges over a received renewed interest in the literature with regards few orders of magnitude (e.g., Steidel et al. 1999; Finkel- to reionization (e.g., Madau & Haardt 2015; McGreer stein et al. 2015; Viironen et al. 2017), and there is grow- et al. 2017). ing evidence of an excess of galaxies over the exponen- Studying both AGN-dominated and star-formation- tially declining bright end of the Schechter function (e.g., dominated UV-emitting galaxies simultaneously is pos- van der Burg et al. 2010; Ono et al. 2018). The uncer- sible given a large enough volume. The Hyper Suprime- tainty at the bright end is due in part to cosmic variance Cam (HSC) Subaru strategic program (SSP) has de- and the small area of past surveys which miss the bright- tected large numbers of both types of objects at z = 4 est galaxies with the lowest surface density. The largest using their optical-only data. However, Akiyama et al. z  4 spectroscopically observed sample used in a pub- (2018) opt to use their excellent ground-based resolution lished luminosity function is from the VIMOS VLT Deep (0.6 ; 4.27 kpc at z = 4) to remove extended sources Survey (VVDS, Le F evre et al. 2013) consisting of 129 and focus on the AGN population separately. Ono et al. 2 0 spectra from  1 deg (Cucciati et al. 2012). This small (2018) selects z  4 galaxies (and AGNs) as g -band sample size limits the analysis of how galaxy growth dropouts in the HSC SSP using strict color cuts includ- properties (e.g., star-formation rate (SFR)) depend on ing the requirement that sources are not signi cantly 0 0 properties like stellar mass and environment, especially detected ( < 2) in the g band. This g band could re- at the bright-end. The large cost of spectroscopically move UV-bright AGNs and could explain why Ono et al. surveying faint sources leaves the most ecient method (2018) nd less sources at M < 24 mag than Akiyama of using multi-wavelength photometry as the best way et al. (2018) (see Figure 7 in Ono et al. 2018). to collect larger samples of star-forming galaxies. For Here, we make use of the 24 deg Spitzer -HETDEX example, a few thousand z = 4 galaxies were detected Exploratory Large Area (SHELA) survey dataset in the four 1 deg elds of the CFHT Survey (van der to probe both AGN-dominated and star-formation- Burg et al. 2010). dominated UV-emitting galaxies over a large area. The Another challenge in measuring the bright end of SHELA dataset includes deep (22.6 AB mag, 50% com- the UV luminosity function is the existence of AGNs plete) 3.6 m and 4.5 m imaging from Spitzer /IRAC 0 0 0 0 0 and their photometric similarities with UV-bright galax- (Papovich et al. 2016) and u g r i z imaging from the ies. The spectral energy distributions (SEDs) of AGN- Dark Energy Camera over 18 deg (DECam; Wold et dominated galaxies are characterized by a power-law al. , in prep). Because SHELA falls within SDSS Stripe- continuum and highly ionized emission lines in the rest- 82 there exists a large library of ancillary data, which frame UV (e.g., Stevans et al. 2014), and like high-z we take advantage of by including in our analysis the UV-bright galaxies the observed SEDs of high-z AGNs VISTA J and Ks photometry from the VICS82 survey exhibit a Lyman break feature due to absorption from (Geach et al. 2017) to help rule out low-z interloping intervening neutral hydrogen in the IGM. Thus, any UV- galaxies. In addition, there is deep X-ray imaging in bright galaxy selection technique relying on the Lyman this eld from the Stripe-82X survey (LaMassa et al. break will also select AGNs. Some have attempted to 2016), which could be used to identify bright AGNs. use a morphological cut to break the color degeneracy We select objects at z > 4 based on photometric crite- of UV-bright galaxies and AGNs by assuming the for- ria. Our sample includes, therefore, both galaxies whose mer will appear extended and the latter will be strictly light is powered by star-formation and AGN activity. Galaxy and AGN Luminosity Functions at z = 4 with SHELA 3 As the bulk of AGNs at z  4 are too faint to detect K = 20:9 mag. Figure 1 shows the lter transmission in the existing X-ray data, we include all z = 4 can- curves for the nine photometric bands used overplot- didate galaxies, regardless of powering source, in our ted with model high-z galaxy spectra illustrating the sample and use our large dynamic range in luminosity { wavelength coverage of our dataset. combined with very bright AGNs from SDSS (Richards 2.2. DECam Data Reduction and Photometric et al. 2006, Akiyama et al. 2018) and very faint galax- Calibration ies from the deeper, narrower Hubble Space Telescope surveys (Finkelstein et al. 2015) { to t the luminosity The DECam images were processed by the NOAO functions of both populations simultaneously. Impor- Community Pipeline (CP). A detailed description of the tantly, our sample is selected using both optical and Community Pipeline reduction procedure can be found Spitzer mid-infrared data, which results in an improved in the DECam Data Handbook on the NOAO website , contamination rate over optical data alone. however, we provide a brief summary of the procedure This paper is organized as follows. The SHELA eld here. First, the DECam images were calibrated using dataset used in this paper is summarized in 2.1. The calibration exposures from the observing run. The main DECam reduction are discussed in Sections 2.2{2.4 and calibration steps included an electronic bias calibration, the IRAC data reduction and photometry in 2.5. Sam- saturation masking, bad pixel masking and interpola- ple selection and contamination are discussed in Section tion, dark count calibration, linearity correction, and 3. Our UV luminosity function is presented in Section 4. at- eld calibration. Next, the images were astrometri- The implications of our results are discussed in section 5. cally calibrated with 2MASS reference images. Finally, We summarize our work and discuss future work in Sec- the images were remapped to a grid where each pixel is tion 6. Throughout this paper we assume a Planck 2013 a square with a side length of 0.27 . Observations taken 1 1 cosmology, with H = 67:8 km s Mpc , = 0:307 on the same night were then co-added. 0 M and = 0:693 (Planck Collaboration et al. 2014). All The CP data products for the SHELA eld were down- magnitudes given are in the AB system (Oke & Gunn loaded from the NOAO Science Archive . The data 1983). products include the co-added images, remapped im- ages, data quality maps (DQMs), exposure time maps 2. DATA REDUCTION AND PHOTOMETRY (ETMs), and weight maps (WMs). The co-added im- In this section, we describe our dataset, image reduc- ages from the CP were not intended for scienti c use, tion, and source extraction procedures. The procedures so we opted to co-add the remapped images. We fol- applied to DECam imaging are largely similar to those lowed the co-adding procedure of Wold et al. (in prep.), used with these data in Wold et al. (in prep), thus we which we summarize here. Using the software package direct the reader there for more detailed information. SWARP (Bertin et al. 2002) the sub-images stored in the FITS les of the remapped images were stitched to- 2.1. Overview of Dataset gether and background subtracted. The remapped im- In this study, we use imaging in nine photometric ages were combined using a weighted mean procedure bands spanning the optical to mid-IR in the SHELA optimized for point-sources. The weighting of each im- Field. The SHELA Field is centered at R.A. = age is a function of the seeing, transparency, and sky h m s  0 00 1 22 00 , declination = +0 00 00 (J2000) and extends brightness and is de ned by Equation A3 in Gawiser approximately 6: 5 in R.A. and 1: 25 in declination. et al. (2006) as 0 0 0 0 0 The optical bands consist of u , g , r , i , and z from the Dark Energy Camera (DECam) and covers 17.8 factor PS 2 w = ; (1) deg of the SHELA footprint (Wold et al. in prep). i scale  rms i i The mid-IR bands include the 3.6 m and 4.5 m from the Infrared Array Camera (IRAC) aboard the Spitzer where scale is the image transparency (de ned as the Space Telescope and covers 24 deg (Papovich et al. median brightness of the bright unsaturated stars after 2016). In addition, we include near-IR photometry in J normalizing the brightness measurement of each star by and K from the February 2017 version of the VISTA- s its median brightness across all exposures), rms is the CFHT Stripe 82 Near-infrared Survey (VICS82; Geach root mean square of the uctuations in background pix- et al. 2017), which covers 85% of the optical imag- els, and factor is de ned as ing footprint and has 5- depths of J = 21:3 mag and http://ast.noao.edu/data/docs 1 3 http://stri-cluster.herts.ac.uk/vics82/ http://archive.noao.edu/ 4 Stevans et al. 1.0 J Ks 0.8 ′ ′ i z [4.5] 0.6 [3.6] 0.4 z=3.5 0.2 z=4.5 0.0 0.40 0.75 1.25 2.10 3.60 4.50 Wavelength (μm) Figure 1. The lter transmission curves for the nine photometric bands used in this study (curves are labelled in the gure) and two model star-forming galaxy spectra with redshifts z = 3:5 and z = 4:5, respectively (black). The model spectra have units of Jy and are arbitrarily scaled. The z = 3:5 galaxy spectrum falls completely red-ward of the u band transmission curve and the z = 4:5 galaxy spectra has almost zero ux falling in the g bandpass. Table 1. DECam Imaging Summary - Seeing and Limiting Magnitude 0 0 0 0 0 u g r i z Sub-Field R.A. Dec. FWHM depth FWHM depth FWHM depth FWHM depth FWHM depth 00 00 00 00 00 ID No. (J2000) (J2000) ( ) (mag) ( ) (mag) ( ) (mag) ( ) (mag) ( ) (mag) h m s  0 00 SHELA-1 1 00 52:8 -0 00 36 1.12 25.2 1.06 24.8 1.0 24.8 0.87 24.5 0.86 23.8 h m s  0 00 SHELA-2 1 07 02:4 -0 00 36 1.2 25.2 1.26 24.8 1.28 24.6 1.39 23.9 0.96 23.6 h m s  0 00 SHELA-3 1 13 12:0 -0 00 36 1.22 25.4 1.36 25.0 1.14 24.7 1.04 24.5 1.21 23.5 h m s  0 00 SHELA-4 1 19 21:6 -0 00 36 1.15 25.3 1.4 24.4 1.05 24.3 1.0 22.1 1.13 23.7 h m s  0 00 SHELA-5 1 25 31:2 -0 00 36 1.21 25.1 1.07 24.9 1.02 24.3 0.93 23.9 0.85 23.6 h m s  0 00 SHELA-6 1 31 40:8 -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 1.89"-diameter apertures on the PSF-matched images (see Section 2.3). The initial rms per pixel was de ned as the inverse of the square-root of the exposure time. The median of FWHM stack factor = 1 exp 1:3 ; (2) the rms map is scaled to the global pixel-to-pixel rms FWHM which is de ned as the standard deviation of the uxes in good-quality, blank sky pixels. Good-quality, blank where FWHM is the median FWMH of bright un- stack sky pixels are pixels not included in a source according saturated stars in an unweighted stacked image and to our initial Source Extractor catalog (see Section 2.3 FWHM is the median FWHM of bright unsaturated for discussion of our source extraction procedure), and stars in each individual exposure. have an exposure time greater than 0.9 times the median The seeing and transparency measurements were de- value. termined using a preliminary source catalog generated The DECam imaging data were photometrically cal- for each resampled image using the Source Extractor ibrated with photometry from the Sloan Digital Sky software package (Bertin & Arnouts 1996).The seeing Survey (SDSS) data release 11 (DR11; Eisenstein et al. in the nal stacked images are listed in Table 1. 2011) using only F0 stars. F0 stars were used because After discovering the original WMs from the Com- their spectral energy distribution span all ve optical munity Pipeline had values inconsistent with the ETM, lters while appearing in the sky at a suciently high we created custom rms maps for the co-added images. Transmission Galaxy and AGN Luminosity Functions at z = 4 with SHELA 5 surface density to provide statistically signi cant num- deconvolution routine with an increasing number of it- bers in each DECam image. We began by creating a erations until the PSF of the convolved image (again preliminary source catalog for the stacked DECam im- measured from stacking stars) had a ux within a 7- ages using Source Extractor and position matching to pixel (1.89 ) diameter aperture matched to that of the the SDSS source catalog. Then we selected F0 stars us- PSF of the target image to within 5%. The total uxes ing SDSS colors by integrating an F0-star model spec- were measured in 30-pixel (8.1 ) diameter circular aper- trum from the 1993 Kurucz Stellar Atmospheres Atlas tures. In Figure 2 we show the results of PSF-matching (Kurucz 1979) with each of the ve optical SDSS lter in each sub- eld by displaying a comparison of the en- curves. For sources in the catalog to be identi ed as larged PSFs to the largest PSF as a percent di erence. an F0 star, the total color di erences, using colors for Due to variations in intrinsic galaxy colors and vari- all adjacent bands, added in quadrature must have been able image depth and sky coverage, some galaxies will less than 0.35. We then calculate the expected mag- not appear in all bands. To get photometric measure- nitude o set between SDSS and DECam lters for F0 ments of all sources in every DECam image we com- 0 0 0 stars, which are as follows: u : 0.33, g : 0.02, r : -0.001, bined the information in the ve optical band images 0 0 i : -0.02, and z : -0.01. The zero-point for each lter into a single detection image. We followed the proce- was then calculated as dure of Szalay et al. (1999) and summed the square of the signal-to-noise ratio in each band pixel-by-pixel as AB ZPT = median(m m m ); (3) DECam o set SDSS follows: band;i where m is the expected magnitude o set between o set D = ; (4) band;i SDSS and DECam lters for F0 stars. After the zero- points were applied to each stacked image, the image where D is the detection image ith pixel value, F i band;i pixel values were converted to units of nJy. is the ith pixel ux in the band image, and  is band;i the rms at that ith pixel pulled from the the band rms 2.3. DECam Photometry image. A weight image associated with this detection Studying galaxy properties relies on accurate mea- image was created where pixels associated with detec- surements of galaxy colors. One way to obtain ac- tion map pixels with data in at least one band have a curate colors is to perform xed-aperture photometry value of unity and pixels associated with detection map where you measure source uxes in every band using pixels without data have a value of zero. the same sized aperture. However, since the DECam im- Photometry was measured on the PSF-matched im- ages where taken in varying seeing conditions they have ages using the Source Extractor software (v2.19.5, point spread functions (PSFs) with a range of full width Bertin & Arnouts 1996). Catalogs were created for at half maximum (FWHM). To perform xed-aperture each of the six SHELA sub- elds with Source Extractor photometry on these images, the PSFs of all the imaging in two image mode using the detection image described covering a single patch of sky must be adjusted to have above, and cycling through the ve DECam bands as a similar PSF size. We divided the DECam imaging the measurement image. In the nal source catalog, we into six sub- elds (each de ned as roughly one DECam maximized the detection of faint sources while minimiz- pointing). For each sub- eld, we enlarged the PSFs of ing false detections by optimizing the combination of the the stacked images to match the PSF of the stacked im- SExtractor parameters DETECT THRESH and DE- age with the largest PSF in that sub- eld. For example, TECT MINAREA. We did this by running SExtractor 0 0 in sub- eld SHELA-1 we matched the PSFs of the g , r , with an array of combinations of DETECT THRESH 0 0 0 i , and z stacked images to the PSF of the u band. To and DETECT MINAREA and chose the combination enlarge the PSFs we adopted the procedure of Finkel- of 1.6 and 3, respectively, which detected all sources stein et al. (2015) who used the IDL deconv tool Lucy- that appeared real by visual inspection and included Richardson deconvolution routine. This routine takes as the fewest false positive detections from random noise inputs two PSFs (the desired larger PSF and the starting uctuations. smaller PSF) and the number of iterations to run and We measure source colors in 1.89 diameter circu- outputs a convolution kernel. The input image PSFs lar apertures (which corresponds to an enclosed ux were produced by median combining the 100 brightest fraction of 59-75% for unresolved sources in our sub- stars (sources with stellar classi cations in SDSS DR11) elds). To obtain the total ux, we derived an aper- in each image. Before combining, the stars were over- ture correction de ned as the ux in a 1.89 diame- sampled by a factor of 10, re-centered, and then binned ter aperture divided by the ux in a Kron aperture by ten to ensure the star centroids aligned. We ran the (i.e., FLUX AUTO), using the default Kron aperture 6 Stevans et al. 20 20 20 g' u' u' Field: 1 Field: 2 Field: 3 15 15 15 r' g' r' Matched to: u' Matched to: i' Matched to: g' i' r' i' 10 10 10 z' z' z' 5 5 5 0 0 0 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 Aperture Diameter (arcsec) Aperture Diameter (arcsec) Aperture Diameter (arcsec) 20 20 20 u' g' u' Field: 4 Field: 5 Field: 6 15 15 15 r' r' r' Matched to: g' Matched to: u' Matched to: g' i' i' i' 10 10 10 z' z' z' 5 5 5 0 0 0 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 Aperture Diameter (arcsec) Aperture Diameter (arcsec) Aperture Diameter (arcsec) Figure 2. The results of PSF-matching the DECam images. Each panel shows the percent di erence between the enlarged PSFs and the largest PSF per SHELA sub- eld. The colored lines correspond to the four bands listed in each panel's legend. The vertical dashed line denotes 1.89" which is the aperture diameter at which we compared PSFs during the PSF-matching procedure (see Section 2.3 for details). The horizontal dashed line was placed to zero to guide the eye. This gure illustrates that all bands in all sub- elds have PSFs that collect the same fraction of light as their respective largest PSF to within 5% except for the z band in sub- eld 3 which matches to about 6%. parameters of PHOT AUTOPARAM= 2.5, 3.5, which compared the extinction-corrected DECam photometry has been shown to calculate the total ux to within 5% to SDSS DR14 per sub- eld and band and found agree- (Finkelstein et al. 2015). This correction was derived in ment for point sources to better than m < 0:05 0:2 the r -band on a per-object basis to account for di erent mag in terms of scatter. source sizes and ellipticities and was applied to the uxes in the other DECam bands per sub- eld. In areas where 2.4. DECam Photometric Errors the sub- elds overlapped, sources with positions that We estimated photometric uncertainties in the DE- matched to within 1.2 in neighboring sub- eld catalogs Cam images by estimating the image noise in aper- had their uxes mean-combined after being weighted by tures as a function of pixels per aperture, N, following the inverse square of their uncertainties (see Section 2.4). the procedure described in Section 2 of Papovich et al. The DECam source uxes were corrected for Galactic (2016). There are two limiting cases for the uncertainty extinction using the color excess E(B-V) measurements in apertures with N pixels,  . If pixel errors are com- by Schla y & Finkbeiner (2011). We obtained E(B-V) pletely uncorrelated, the aperture uncertainty scales as values using the Galactic Dust Reddening and Extinc- the square root of the number of pixels,  =   N , N 1 tion application on the NASA/IPAC Infrared Science where  is the standard deviation of sky background 4 1 Archive (IRAS) website . We queried IRAS for E(B-V) pixels. If pixel errors are completely correlated then values (the mean value within a 5 radius) for a grid =   N (Quadri et al. 2006). Thus, the aperture 0 N 1 of points across the SHELA eld with 4 spacing and uncertainty will scale as N with 0:5 < < 1. assigned each source the E(B-V) value from the closest To estimate the aperture noise as a function of N grid point. The Cardelli et al. (1989) Milky Way redden- pixels, we measured the sky counts in 80,000 randomly ing curve parameterized by R = 3:1 was used to derive 00 00 placed apertures ranging in diameter from 0: 27 to 8: 1 the corrections at each band's central wavelength. We across each stacked DECam image. We required aper- tures to fall in regions of the background sky, which we http://irsa.ipac.caltech.edu/applications/DUST/ de ne as the region where the exposure time map has the value of the at least the median exposure time (en- % Diff. in Total Fractional Flux % Diff. in Total Fractional Flux % Diff. in Total Fractional Flux % Diff. in Total Fractional Flux % Diff. in Total Fractional Flux % Diff. in Total Fractional Flux Galaxy and AGN Luminosity Functions at z = 4 with SHELA 7 Table 2. Fit Parameters for Background Fluctuations as Function β = 1 of Aperture Size Using Eq.(5) y 800 i Sub-Field ID Band   rms 1 med (nJy) (nJy) SHELA-1 u 5.02 0.09 0.89 1.33 0.35 5.58 β = 0.5 g 4.07 0.21 0.83 1.71 0.43 8.85 σ r 4.14 0.18 0.91 2.37 0.36 10.2 i 4.1 0.13 1.0 2.88 0.42 14.5 200 z 7.16 0.19 0.93 2.55 0.45 25.0 SHELA-2 u 1.21 0.24 0.91 2.41 0.55 5.88 0 0 5 10 15 20 25 g 1.8 0.35 0.87 2.3 0.5 8.26 r 3.02 0.21 0.94 2.68 0.41 10.2 i 14.3 0.05 1.0 1.37 0.36 15.8 0 Figure 3. Background noise uctuations,  , in an aperture z 5.36 0.34 0.9 2.48 0.51 27.7 of N pixels plotted as a function of the square root of the SHELA-3 u 1.68 0.23 0.89 2.2 0.42 4.81 number of pixels for the ve DECam band images in sub- g 5.36 0.12 0.88 1.33 0.33 6.25 eld SHELA-1. The colored dotted lines are the measured r 2.77 0.19 0.94 2.63 0.45 9.8 aperture uxes and the solid lines are ts to the data. See legend insert for color coding information. The dot-dashed i 3.06 0.32 0.92 3.07 0.38 12.6 0 line shows the relation assuming uncorrelated pixels, z 11.0 0.15 0.91 1.95 0.41 24.1 p N . The dashed line shows the relation assuming perfectly SHELA-4 u 1.08 0.26 0.95 2.54 0.49 4.64 correlated pixels (  N ; Quadri et al. 2006). g 7.36 0.13 0.85 1.24 0.34 8.31 r 2.86 0.25 0.91 2.44 0.51 12.1 measured ux uncertainty in a given aperture,  , as i 24.9 1.42 0.6 1.27 0.6 63.5 a function of the square root of the number of pixels z 5.48 0.26 0.93 2.79 0.45 22.6 in each aperture, N , for the ve DECam bands in the SHELA-5 u 5.25 0.11 0.87 1.25 0.33 5.8 sub- eld SHELA-1. g 3.11 0.23 0.87 1.92 0.43 7.86 Following Papovich et al. 2016 we t a parameterized r 4.78 0.38 0.81 2.02 0.46 15.8 function to the noise in an aperture of N pixels,  , as, i 5.53 0.33 0.85 2.32 0.46 21.9 z 6.86 0.23 0.92 2.53 0.49 28.0 =  ( N + N ); (5) N 1 SHELA-6 u 1.66 0.15 0.95 2.38 0.42 4.6 g 6.15 0.12 0.86 1.27 0.32 6.8 where  is the pixel-to-pixel standard deviation in the r 4.81 0.21 0.87 2.09 0.4 12.1 sky background, and , , , and  are free parameters. i 6.87 0.14 0.93 2.14 0.37 14.7 The best- tting parameters in Equation 5 for the com- bined DECam images are listed in Table 2. While the z 6.09 0.51 0.82 2.13 0.53 32.8 second term was intended to aid in tting the data at Note| Typical values of  0.65-0.70 when using a two parameter large N values, in actuality the second term contributed t (i.e. with only and ) suggest slightly correlated noise between signi cantly at all N values resulting in 0.9-1. Nev- pixels. ertheless our functional ts reproduce the data well as can be seen, for example, in Figure 3. To estimate how correlated the pixel-to-pixel noise is, we t  with only the rst term in Equation 5 and found typical values suring >50% of each image was considered), excluding of  0.65-0.70 suggesting slightly correlated pixel-to- detected sources and pixels agged in the DQM. We also pixel noise. required the apertures do not overlap with each other. The photometric errors estimated by Equation 5 were We then estimated  for each aperture with N pixels scaled to apply to ux measurements outside the region by computing the standard deviation of the distribution with the median exposure time. The ux uncertainty for of aperture uxes from the normalized median absolute the i-th source in band b in the sub- eld f is calculated deviation,  (Beers et al. 1990). We calculated nmad 1 as by computing  for all pixels in the background sky ! nmad rms as de ned above. Figure 3 shows an example of the i;f;b 2 2 =  ; (6) i;f;b N;f;b rms med;f;b (nJy) N 8 Stevans et al. where  is given by Equation 5 for each band and uxes is closest to the actual image pixels, with respect N;f;b sub- eld, rms is the value of the rms map at the to the noise model. We describe how we use the Tractor i;f;b central pixel location of the i-th source in each band to perform forced photometry on IRAC images in detail and sub- eld, and rms is the median value of the below. med;f;b rms map in each band and sub- eld. The photometric We begin our source modeling procedure by selecting error estimates exclude Poisson photon errors, which we the ducial band high-resolution model of each source. estimate to contribute <5% uncertainty to the optical We use the uxes and surface brightness pro le shape uxes of our high-z candidates. parameters measured in our DECam detection image be- As described in Section 2.3, all source uxes are mea- cause the image combines the information of all sources sured in circular apertures of 1.89 diameter and scaled in the ve optical band images (as described in Section to total on a per-object basis. Likewise, the ux uncer- 2.3). Second, we use one-component circular Gaussian tainties in the apertures were scaled by the same amount to model the PSF. During the modeling of each source, to determine the total ux uncertainties, so that the S/N we allow Tractor to optimize the Gaussian  value, in for the total ux is the same as for the aperture ux. addition to optimize a source ux. We nd that the me- dian of the optimized Gaussian  is 0: 80 (equivalent to a full-width at half maximum of 1: 88) for both IRAC 2.5. IRAC data reduction and photometry 3.6 and 4.5 m images, consistent with those measured As discussed below, we wish to enhance the validity of from an empirical point response function for the 3.6m our z = 4 galaxy sample by including IRAC photometry and 4.5m image (FWHM of 1: 97, see Papovich 2016, in our galaxy sample selection. While we could allow Section 3.4). this by position-matching the published Spitzer/IRAC In practice, we extracted an IRAC image cutout of catalog from Papovich et al. (2016) to our DECam cat- each source in the input DECam catalog. We selected alog, this is not optimal for two reasons. First, the Pa- 00 00 the cutout size of 16  16 . This cutout size represents povich et al. (2016) catalog is IRAC-detected, and so a trade-o between minimizing computational costs re- only includes sources with signi cant IRAC ux, while lated to larger cutout sizes and ensuring that the sources for our purposes, even a non-detection in IRAC can be lie well within the cutout extent. The sources of interest useful for calculating a photometric redshift. Second, within cutout are modeled as either unresolved (i.e., a this catalog uses apertures de ned on the positions and point source) or resolved based on the DECam detec- shapes of the IRAC sources, while the larger PSF of tion image. We considered a source to be resolved if the IRAC data results in signi cant blending, which is an estimated radius r > 0: 1. We de ne the radius as, a larger issue at fainter magnitudes, where we expect r = a  b=a, where a is a semi-major axis and source to nd the bulk of our sources of interest. For these b=a is an axis ratio. We perform the photometry for reasons, we applied the Tractor image modeling code resolved sources using a deVaucouleurs pro le (equiva- (Lang et al. 2016a,b) to perform \forced photometry", lent to S ersic pro le with n=4) with shape parameters which employs prior measurements of source positions (semi-major axis, position angle, and axis ratio) mea- and surface brightness pro les from a high-resolution sured using our DECam detection image. We have also band to model and t the uxes of the source in the re- performed the photometry using an exponential pro le maining bands, splitting the ux in overlapping objects (equivalent to S ersic pro le with n=1), but we do not into their respective sources. We speci cally used the nd any signi cant di erence between the IRAC ux Tractor to optimize the likelihood for the photometric measurements for the two galaxy pro les. Therefore, properties of DECam sources in each of IRAC 3.6 and we adopt a deVaucouleurs pro le to model all resolved 4.5 m bands given initial information on the source sources. The Tractor simultaneously modeled and op- and image parameters. The input image parameters of timized the sources of interest and neighboring sources IRAC 3.6 and 4.5 m images included a noise mode, within the cutout. Finally, the Tractor provided the a point spread function (PSF) model, image astromet- measurement IRAC ux of each DECam source with ric information (WCS), and calibration information (the the lowest reduced chi-squared value. We validated the \sky noise" or rms of the image background). The in- Tractor-based IRAC uxes by comparing the uxes of put source parameters included the DECam source posi- isolated sources (no neighbors within 3 ) to the pub- tions, brightness, and surface brightness pro le shapes. lished Spitzer/IRAC catalog from Papovich et al. (2016). The Tractor proceeds by rendering a model of a galaxy For both bands, we found good agreement with a bias or a point source convolved with the image PSF model o set of m < 0:05 mag and a scatter of <0.11 mag at each IRAC band and then performs a linear least- square t for source uxes such that the sum of source Galaxy and AGN Luminosity Functions at z = 4 with SHELA 9 down to m = 20:5 mag and a bias o set of m < 0:13 overcome this we required a source to pass the z se- phot mag and a scatter of <0.26 down to m = 22 mag. lection process with and without including the VISTA J and K photometry. This double z selection pro- phot 3. SAMPLE SELECTION cess resulted in an initial sample of 4,364 high-z objects. Next, we performed an investigation into possible con- 3.1. Photometric Redshifts and Selection Criteria tamination, which resulted in additional selection crite- We selected our sample of high-redshift galaxies using ria and a re ned sample. a selection procedure that relies on photometric redshift (z ) tting, which leverages the combined information phot 3.2. Identifying Contamination in all photometric bands used. We obtained z 's and phot Photometric studies of high-z objects can be contam- z probability distribution functions (PDFs) from the phot inated by galactic stars and low-z galaxies whose 4000 EAZY software package (Brammer et al. 2008). For A break can mimic the Lyman- break of high-z galax- this analysis, we use the \z a" redshift column from 2 ies. The inclusion of IRAC photometry is crucially im- EAZY which is produced by minimizing the  in the portant for removing galactic stars from our sample as all-template linear combination mode. We also tried the galactic stars have optical colors very similar to z = 4 \z peak" column and our resulting luminosity function objects (Figure 5). While inspecting the photometry did not change signi cantly. EAZY assumes the inter- and best- tting SED templates to con rm our ts were galactic medium (IGM) prescription of Madau (1995). robust and our high-z galaxy candidates were convincing We did not use any magnitude priors based on galaxy we found evidence of contamination in our preliminary luminosity functions when running EAZY as the exist- photometrically selected sample. We explored ways of ing uncertainties in the bright-end would bias our re- identifying and removing the contamination including sults (e ectively assuming a at prior). We then ap- cross-matching our sample with proper motion catalogs, plied a number of selection criteria using the z PDFs phot SDSS spectroscopy, and X-ray catalogs. Additionally, from EAZY following the procedure of Finkelstein et al. we implemented machine learning methods. (2015), which we summarize here, to construct a z 4 galaxy sample. First, we required a source to have a 3.2.1. Cross-matching with NOMAD and SDSS signal-to-noise ratio of greater than 3.5 in the two pho- 0 0 tometric bands (r and i bands) that probe the UV While inspecting the photometry of our preliminary continuum of our galaxy sample, which has been shown sample derived before the inclusion of the VISTA data 0 0 to limit contamination by noise to negligible amounts by we found a fraction of candidates had red r i colors 0 0 Finkelstein et al. (2015). Next we required the integral and excesses in the i and z bands with respect to the of each source's z PDF for z > 2:5 to be greater than best tting z  4 template. We investigated whether phot 0.8 and the integral of the z PDF from z = 3.5-4.5 these objects had low redshift origins by cross-matching phot to be greater than the integral of the z PDF in all our sample with the Naval Observatory Merged As- phot other z = 1 bins centered on integer-valued redshifts. trometric Dataset (NOMAD) proper motion catalog Finally, we required a source to have no photometric (Zacharias et al. 2004) and the SDSS spectral catalog 0 0 0 0 ags on its u , g , r , and i ux measurements in our (DR13; Albareti et al. 2017). The cross-matching with photometric catalog. These four bands are crucial for NOMAD resulted in the identi cation of 16 objects with probing both sides of the Lyman break of galaxies in proper motion measurements, 6 of which were large in the redshift range of interest. Photometric ags indi- magnitude suggesting that these sources were stellar cate saturated pixels, transient sources, or bad pixels as contaminants. Cross-matching with SDSS spectra re- de ned by the CP (see section 2.2). sulted in the identi cation of 23 z  4 QSOs, two low Initially, we required sources to pass this z selec- mass stars, and one low-z galaxy in our sample. All of phot 0 0 tion procedure when only the DECam and IRAC bands the stellar objects and the low-z galaxy had red r i were used. After inspecting candidates from this selec- colors con rming our suspicion that a fraction of our tion and nding low-z contaminants (see Section 3.2) we sample had low-z origins. We removed from our candi- moved to include the VISTA J and K photometry when date sample the 6 objects with proper motions above 50 available. The candidate sample of high-z galaxies after mas/yr in NOMAD. We chose the threshold 50 mas/yr including the VISTA J and K photometry signi cantly because some SDSS QSOs are reported to have small reduced the amount of obvious low-z contamination as nonzero proper motions in the NOMAD catalog. After described in the following section, however, a new class inclusion of the VISTA J and K photometry, 5 of the of contaminant appeared in the candidate sample due 6 rejected objects became best- t by a low-z solution to erroneous VISTA photometry for some sources. To and therefore were rejected by our selection process au- 10 Stevans et al. tomatically. Results like this suggested the inclusion of measurements. The obviously high-z objects and the the VISTA data improved the delity of our sample. SDSS QSOs had roughly at rest-UV spectral slopes 0 0 0 (i.e., the i , r , and z uxes were comparable) while the 3.2.2. Cross-matching with X-ray catalog Stripe 82X obviously low-z objects had a clear spectral peak be- 0 0 In principle, AGNs can be distinguished from UV- tween the optical and mid-IR, speci cally the r i and 0 0 0 bright star-forming galaxies at high z by their X-ray i z colors were quite red while the z [3:6] was blue. emission as AGNs dominate X-ray number counts down We then selected three data "features" or measurements 17 2 1 to F = 1  10 erg cm s in 0.5-2 keV Chan- for the decision tree to choose from to predict the clas- dra imaging (Lehmer et al. 2012). In fact, Giallongo si cations of the training set: 1. The  of the best 0 0 et al. (2015) found 22 AGN candidates at z > 4 by tting high-z galaxy template, 2. The r i color, and measuring uxes in deep 0.5-2 keV Chandra imaging 3. The signal-to-noise ratio in the u band. We include at the positions of their optically-selected candidates. this signal-to-noise ratio because most obviously high- They required AGN candidates to have an X-ray detec- z sources had a signal-to-noise ratio of less than 3 (as 17 2 1 tion of F > 1:5 10 erg cm s . We attempted they should, as this lter is completely blue-ward of the to identify AGNs in our candidate sample by position Lyman break for the target redshift range) while the ob- matching with the 31 deg Stripe 82X X-ray catalog viously low-z source did not. We restricted the training (LaMassa et al. 2016). Unfortunately, the ux limit set to include only the obviously high-z objects (includ- 16 2 1 at 0.5-2 keV was 8:7  10 erg cm s , shallower ing SDSS classi ed QSOs) and obviously low-z objects. than the detection threshold used by Giallongo et al. To evaluate the decision tree performance we assigned (2015). Cross-matching the Stripe 82X X-ray Catalog 34% of the training set to a test set and left the test set with our candidate sample resulted in 8 matches within out of the rst round of training. The rst round of 7 or the matching radius used by LaMassa et al. (2016) training was validated using three-fold cross-validation to match the Stripe 82X X-ray Catalog with ancillary with a score of 94+/-2%. Three-fold cross validation data. Two of the X-ray-matched sources (m 19) are veri es the model is not over- tting or overly dependent SDSS AGNs, although one has two extended objects on a small randomly selected training (or validation) set. within the 7 matching radius, drawing into question The process of three-fold cross-validation involves split- the likelihood of the match. Another X-ray-matched ting the training sample into three sub-samples, train- 0 0 source (m 22) is a very red object (r i = 0.9 and ing the model independently on each combination of i [3:6] =3.4) without SDSS spectroscopy. The remain- two sub-samples while validating on the remaining sub- ing ve X-ray-matched sources are fainter (m 24) and sample, and then averaging the validation scores of all without SDSS spectroscopy, and two have large separa- the combinations. A validation score of 100% indicates tions (> 6 ) with their X-ray counterpart. All 8 X-ray- each sub-sample combination trained a model that suc- matched sources satis ed each of our selection criteria cessfully predicted the classi cations of every object in and made it into our nal sample. the remaining sub-sample. After validation, we applied the model to the test set and achieved a test score of 3.2.3. Insights from machine learning 91%. We then incorporated the test set into the train- To further understand and minimize the contamina- ing set and retrained the model. This second round tion in our preliminary sample we utilized two machine of training had a 3-fold cross validation score of 94 learning algorithms: a decision tree algorithm and a ran- 5%. The algorithm determined that classi cation of the dom forest algorithm. First, we tried the decision tree training set could be predicted at 92% accuracy using algorithm from the sci-kit learn Python package (Pe- 0 0 2 the r i color and the  only, where obviously high-z dregosa et al. 2011). Executing the decision tree algo- 0 0 2 objects have r i < 0.555 and  < 70.6. rithm involved creating a training set by classifying (via After performing the decision tree machine learning visual inspection) the 300 brightest objects as one of we tried the random forest algorithm by using the scikit- ve types: 1. obviously high-z galaxy or high-z AGN, learn routine RandomForestClassi er (Pedregosa et al. 2. obviously low-z galaxy, 3. indistinguishable between 2011). The random forest algorithm uses the combined high-z or low-z object, 4. SDSS spectroscopically iden- information of all the features of our dataset instead ti ed QSO, 5. spurious source. The classi cation was of using only the three well-motivated features that we driven by a combination of the shape of the optical SED, provided the decision tree algorithm. The random for- the  of the best tting high-z (z  4) galaxy tem- est algorithm ts a number of decision tree classi ers on plate, the di erence in  between the best tting high-z a randomly drawn and bootstrapped subset of the data galaxy template and the best tting low-z galaxy tem- using a randomly drawn subset of the dataset features. plate, SDSS spectral classi cation, and proper motion Galaxy and AGN Luminosity Functions at z = 4 with SHELA 11 The decision tree classi ers are averaged to maximize ac- duced the source uxes randomly to create a at dis- curacy and control over- tting. To provide the best clas- tribution of dimmed r -band mag spanning the range si cations to the algorithm we re-inspected the photo- of our candidate high-z galaxies (m = 18-26). We metry and the besting- tting EAZY template SEDs of assigned ux uncertainties to the dimmed uxes by ran- the brightest 311 objects (m < 22:9) and classi ed by domly drawing from the ux uncertainties of our can- eye each with a probability of being at high z, at low z, didate high-z galaxy sample. We then perturbed the and a spurious source. We also inspected the best- tting dimmed uxes simulating photometric scatter by draw- EAZY template at the redshift of the second highest ing ux perturbations from a Gaussian distribution with peak in the redshift PDF, which was usually between a standard deviation equal to the assigned ux uncer- z = 0:1 0:5. The features we provided the algorithm tainties. Since the VISTA J and K photometry from included all colors using the DECam and IRAC pho- the VICS82 Survey covered 85% of the total SHELA tometric bands, the  value of the best- tting EAZY survey area, 15% of our high-z galaxy candidates do template, and the u-band S/N ratio. We ran the random not have ux measurements in the VISTA bands. We forest algorithm using 1000 estimators (or decision trees) incorporated this property of our catalog into the mock with a max depth of two and balanced class weights. catalog by randomly omitting dimmed J and K uxes While the classi cation accuracy of the random forest at the same rate as the fraction of missing J and K was only marginally better than the decision tree, we uxes in our nal sample for a given m 0 bin. were able to determine the relative importance of the We created 4.8 million arti cially dimmed sources and features within the random forest. The most important ran them through EAZY and our selection criteria in an feature was the i [3:6] color, which was consistent identical manner as our real catalog. We found 30,594 with our prior observations of the low-z contaminants arti cially dimmed sources satis ed our z = 4 selec- phot having blue z [3:6] while the SDSS classi ed AGNs tion criteria. We inspected the undimmed and unper- had at or red z [3:6]. We compared the predictive turbed properties of these sources and found a range of 0 0 0 0 power of the random forest algorithm to that of a simple r i colors (0 < r i < 1:4), and the objects with the 0 0 0 i [3:6] color cut{where high-z candidates were required reddest r i colors contaminated at the highest rate. 0 0 0 to have an i [3:6] > 0:2{and found the simple color We de ned the contamination rate in the i-th r i cut to be just as predictive. We, therefore, elected to color bin as adopt the simple i [3:6] color cut as an additional dimmed,selected,i step in our selection process. After we applied this cut R = ; (7) our high-z object sample consisted of 3,833 objects. We N dimmed,i then inspected the brightest 3,200 brightest candidates where N is the number of dimmed dimmed,selected,i (m 0 < 24:0) and removed 61 spurious sources cutting 0 0 sources satisfying our z = 4 criteria in the i-th r i phot our candidate sample to 3,772 objects. The nal step in 0 0 color bin and N is the total number of dimmed our candidate selection process was an r i color cut, dimmed,i 0 0 sources in the i-th r i color bin. We created seven which was a result of the contamination test described 0 0 0 0 r i color bins spanning the range of r i colors in the following subsection. 0 0 recovered (0 < r i < 1.4) with each bin having width=0.2 mag. The contamination fraction, F , in our 3.3. Estimating Contamination Using Dimmed Real high-z galaxy sample was de ned as Sources R N P i We estimated the contamination in our high-z galaxy total;i (1R ) sample by simulating faint and low-S/N low-z interloper F = ; (8) galaxies following a procedure used by Finkelstein et al. (2015). We did this by selecting a sample of bright where R is the contamination rate de ned by Equation low-z sources from our catalog, dimming and perturbing 7, N is the total number of sources in our object total;i 0 0 their uxes, and assigned the appropriate uncertainties catalog in the i-th r i color bin with 0:1 < z < 0:6 phot to their dimmed uxes. This empirical test assumes that and N is the number of high-z candidates in our nal bright low-z galaxies have the same colors as the faint sample in a given redshift bin. We calculated F as a 0 0 0 lower-z galaxies. The sources we dimmed all had m = function of m and m and found contamination frac- r r i 14.3 - 18 mag, brighter than our brightest z =3.5-4.5 tions of between 0-20% generally increasing as a function candidate source, a best- tting z of 0.1 < z < of magnitude and not exceeding 25% in any magnitude phot phot 0.6, and no photometric ags in any optical band (see bin above our 50%-completeness limit (m 0 > 23:5). We section 2.2 for de nition of photometric ag). We re- learned we can improve our delity by implementing a 12 Stevans et al. Table 3. Summary of z = 4 Sample Selection Criteria 3.4. Completeness Simulations A crucial component of calculating the UV luminosity Criterion Section Reference function is measuring the e ective volume of the survey. S/N > 3:5 (3.1) The e ective volume measurement depends on quanti- S/N > 3:5 (3.1) fying the survey incompleteness due to image depth and PDF(z) dz > 0:8 (3.1) selection e ects. To quantify the survey incompleteness 2:5 4:5 PDF(z) dz > All other z = 1 bins (3.1) we simulated a diverse population of high-z galaxies with 3:5 i [3:6] > 0:2 (3.2.3) assigned photometric properties and uncertainties con- 0 0 r i < 1:0 (3.3) sistent with our source catalog and measured the frac- tion of simulated sources that satis ed our high-z galaxy selection criteria as a function of the absolute magnitude and redshift. The simulated mock galaxies were given properties 0 0 r i < 1:0 color cut, given that the objects with the 0 0 drawn from distributions in redshift and dust atten- reddest r i colors contaminated at the highest rate. uation (e.g., E[B-V]) while the ages and metallicities We then re-ran our contamination simulation with our were xed at 0.2 Gyr and solar (Z = Z ), respectively. nal set of selection criteria and found improvement by Because the fraction of recovered galaxies per redshift a few percentage points in two of our brighter magnitude bin is independent of the number of simulated galax- bins (m = 18.75, 20.5) and, again, contamination frac- ies per redshift bin as long as low-number statistics are tions of between 0-20% generally increasing as a function avoided, the redshift distribution was de ned to be at of magnitude and not exceeding 25% in any magnitude from 2 < z < 6, and the E(B-V) distribution was de ned bin brighter than our 50%-completeness limit as shown 0 0 to be log-normal spanning 0 < E(B-V) < 1 and peaking in Figure 6. The addition of the r i < 1:0 color cut re- at 0.2. Mock SEDs were then generated for each galaxy duced our high-z candidate sample size from 3,772 to the using pythonFSPS (Foreman-Mackey et al. 2014), a nal size of 3,740 with a median z of 3.8. The mea- phot 0 python package that calls the Flexible Stellar Popula- sured i -band magnitude and the redshift distribution tion Synthesis Fortran library (Conroy et al. 2009; Con- of our nal sample is shown in Figure 4. A summary of roy & Gunn 2010). We then integrated each galaxy our sample selection criteria is listed in Table 3. Bright 0 0 0 0 SED through our nine lters (DECam u , g , r , i , and (m 0 < 22) candidates are plotted in Figure 5 showing z ; VISTA J and K ; and IRAC 3.6 m and 4.5 m). the distinct color parameter space they occupy (along Each set of mock photometry was then scaled to have with SDSS spectroscopically classi ed AGNs) compared an r-band apparent magnitude within a log-normal dis- to the parameter space occupied by SDSS spectroscopi- tribution spanning 18 < m 0 < 27. This distribution cally classi ed stars in the SHELA eld. ensured we were simulating the most galaxies at the fainter magnitudes where we expected to be incomplete 4 3 and fewer at bright magnitudes where we expected to be 10 10 very complete. Flux uncertainties and missing VISTA 2 uxes were assigned in the same way as during the con- 2 tamination test in Section 3.3. Mock galaxies were not assigned photometric ags. This likely results in an in- completeness of only a couple percent as the fraction of all sources a ected by ags is small. All sources with a signal-to-noise ratio of greater than 3.5 in any one band, 20 22 24 26 3.0 3.5 4.0 4.5 m are agged in another band less than 2% of the time on i Redshift average (never exceeding 4%). In addition all images have 0.5% of all pixels agged on average (never ex- Figure 4. The i -band magnitude distribution (left) and the ceeding 2%). photometric redshift distribution (right) of our nal sample of z  4 candidates. Photometric redshifts are the redshifts The photometric catalog of 600,000 mock galaxies was where  is minimized for the all-template linear combina- then run through EAZY to generate z PDFs and our phot tion mode from the EAZY software (z a). We have found high-z galaxy selection was applied. The completeness high-z sources over a wide range of brightnesses and across the entire z =3.5-4.5 range. http://dfm.io/python-fsps/current/ Number Number Galaxy and AGN Luminosity Functions at z = 4 with SHELA 13 3.0 2 AGN (3. 2 ≤ z < 3. 5) AGN (3. 2 ≤ z < 3. 5) spec spec AGN (4. 5 ≤ z < 4. 8) AGN (4. 5 ≤ z < 4. 8) spec spec 2.5 2.0 1.5 1.0 −1 Final Candidates Final Candidates 0.5 (m ′ < 22) (m ′ < 22) i i Candidates omitted Candidates omitted ′ ′ b colors (m < 22) b colors (m < 22) i i −2 0.0 ′ ′ Stars (S/N > 100) Stars (S/N > 100) i i AGN (3. 5 ≤ z < 4. 5) AGN (3. 5 ≤ z < 4. 5) spec spec −0.5 −0.5 0.0 0.5 1.0 1.5 2.0 2.5 −0.5 0.0 0.5 1.0 1.5 2.0 2.5 ′ ′ ′ ′ r − i [mag] r − i [mag] 0 0 0 0 Figure 5. Left: Color-color plot showing DECam g r vs r i . Blue points correspond to bright (S=N > 100) sources classi ed as stars within SDSS DR14. Sources spectroscopic identi ed as QSO within SDSS DR13 at 3:2  z < 3:5, 3:5  z < 4:5, and 4:5  z < 4:8 are orange circles, black \x"s, and green circles, respectively. Bright candidate objects from our study are 0 0 0 0 shown as squares (see legend insert for color coding). Right: Color-color plot showing i [3:6] vs r i . The r i < 1:0 and the i [3:6] > 0:2 color selection criteria are denoted as the vertical dashed line and the horizontal dashed line, respectively. No candidates are hidden by the legend inserts. These plots illustrate how the inclusion of IRAC photometry breaks the optical color degeneracy of z = 4 sources and galactic stars. In the left panel where only optical colors are used the z  4 AGNs share the parameter space of galactic stars, while in the right panel where one color includes the [3:6] IRAC band there is a larger separation. was de ned as the number of mock galaxies recovered magnitude at rest-frame 1500 A (M ) using the fol- UV;i divided by the number of input mock galaxies, as a func- lowing formula tion of input absolute magnitude and redshift. Figure 7 shows the results of our simulation. We de ne the 50%- M 0 = m 5 log(d =10pc) + 2:5 log(1 + z); (10) UV;i i L completeness limit as the absolute magnitude where the area under the curve falls to less than 50% the areas un- where d is the luminosity distance in pc. The second der the average of the M 0 = 25 to M 0 = 28 and third terms of the right side are the distance mod- UV;i UV;i curves, which we nd to be at M = 22 (m = 24 UV;i ulus. for z = 4). Our UV luminosity function is shown as red diamonds in Figure 8. We note that we do not include the lumi- 4. RESULTS nosity function data points in bins M 0 -22 in our UV;i analysis as these bins are below our 50%-completeness 4.1. The Rest-Frame UV z=4 Luminosity Function limit as discussed in Section 3.4, where the complete- We utilize the e ective volume method to correct for ness corrections are unreliable due to the low S/N of incompleteness in deriving our luminosity function. The our data. This is con rmed by comparing to the HST e ective volume (V ) can be estimated as eff CANDELS results in these same magnitude bins, which are more reliable due to their higher S/N. In Figure 8 dV V (M 0 ) = C (M 0; z)dz; (9) eff i UV;i we also include the UV luminosity function of z = 4 dz UV-selected galaxies from deeper Hubble imagining by dV where is the co-moving volume element, which de- Finkelstein et al. (2015) as green squares. In the bin dz 0 0 pends on the adopted cosmology, and C (M ; z) is the where we overlap with this dataset (M = -22.5) UV;i UV;i completeness as calculated in Section 3.4. The integral both results are consistent, however, if the luminosity was evaluated over z = 3 5. function derived from Hubble imaging is extrapolated to To calculate the luminosity function, we convert the brighter magnitudes, it would fall o more steeply than apparent i-band AB magnitudes (m 0 ) to the absolute our luminosity function. While our luminosity function ′ ′ g r [mag] i [3. 6] [mag] 14 Stevans et al. P(z) only UV, i 1.0 1.0 P(z) + i − [3.6] -28.0 0 0 0 P(z) + i − [3.6] + r − i -27.0 50%-Completeness Limit 0.8 -26.0 0.8 -25.0 -24.0 0.6 0.6 -23.0 -22.5 -22.0 0.4 0.4 -21.5 0.2 0.2 0.0 0.0 3.0 3.5 4.0 4.5 5.0 18 19 20 21 22 23 24 25 Input z Apparent Magnitude (mag) Figure 6. The results of our contamination simulations Figure 7. The results of our completeness simulations, estimating the contamination fraction using dimmed real showing the fraction of simulated sources recovered as a func- sources, showing the contamination fraction, F , as a func- tion of input redshift. Each colored line is for the correspond- 0 0 tion of m (solid and dashed colored lines) and m (dotted ing 0.5 magnitude bin according to the legend. We de ne i r colored lines). The color of the lines represent the selection the 50%-completeness limit as the absolute magnitude curve criteria applied before F was calculated: blue for F after with an integral value of less than 50% the area under the only the z PDF selection cuts, purple for F after the average of the M = 25 to M = 28 curves. We nd phot 1450 1450 z PDF selection and the i [3:6] color cuts, and orange the 50%-completeness to be M = 22:0. phot 1450 for F after the z PDF selection, the i [3:6] color, and the phot 0 0 r i color cuts. The error bars indicate the standard devia- Hubble imaging and the SDSS AGN number densities, tion of the mean from subdividing our 4.8 million simulated our data can potentially provide a robust measurement sources into 20 subsamples. The nal contamination frac- of the bright-end slope of the star-forming galaxy lu- tion was found to be between 0-20% generally increasing as minosity function and the faint-end slope of the AGN a function of magnitude and not exceeding 25% in any mag- nitude bin brighter than our 50%-completeness limit m 0 luminosity function, which we explore in the following 23.5 (black solid line). section. 4.2. Fitting the Luminosity Function 0 0 declines from M = 22 to M = 24, it attens UV;i UV;i out at brighter magnitudes before turning over again. With our luminosity function in agreement with the As our sample includes all galaxies which exhibit a faint end of the AGN luminosity function from SDSS Lyman break, we expect that this attening is due to DR7 and the bright end of the star-forming galaxy lu- the increasing importance of AGN at these magnitudes. minosity function from CANDELS, we attempt to si- The large volume surveyed by SDSS data has led to multaneously t empirically motivated functions to both components. For the AGN component we use a DPL the selection and spectroscopic follow-up of AGNs at many redshifts. SDSS AGN studies have found the AGN function motivated by the AGN UV luminosity func- tion work on large, homogeneous quasar samples (e.g., UV luminosity function to exhibit a double power law (DPL) shape (e.g., Richards et al. 2006). In Figure 8 we Boyle et al. 2000; Richards et al. 2006; Croom et al. show as red "x"s the z = 4 AGN UV number densities 2009; Hopkins et al. 2007, and references therein). The derived from the SDSS DR7 catalog (Schneider et al. function form of a DPL follows 2010) by Akiyama et al. (2018), who select AGNs to M > 28:9. We can see that at the magnitudes UV;i (M ) = ; (11) where we overlap 27 . M 0 . 26, the agreement UV;i 0:4( +1)(MM ) 0:4( +1)(MM ) 10 + 10 is excellent with our data. The only di erence is at where  is the overall normalization, M is the char- M = 28, where our survey detects two quasars when acteristic magnitude, is the faint-end slope, and is the AGN luminosity function by Akiyama et al. (2018) the bright-end slope. would predict less than one in our volume. We attribute this di erence to cosmic variance. By combining our For the star-forming galaxy UV luminosity function we consider sepraratley both a Schechter function and a data with the star-forming galaxy number densities from DPL, as well as including magni cation via gravitational Contamination Fraction, F Completeness Galaxy and AGN Luminosity Functions at z = 4 with SHELA 15 -2 3.0 -3 2.0 1.0 0.0 -4 27 26 25 24 23 UV, i -5 -6 DPL+DPL Fit Galaxy DPL Component -7 AGN DPL Component DPL+Sch Fit Galaxy Schechter Component -8 AGN DPL Component This Work This Work -9 (<50% complete) Finkelstein+15 SDSS DR7 (Akiyama+18) -10 Ono+18 28 26 24 22 20 18 UV, i Figure 8. The rest-frame UV z = 4 luminosity function of star-forming galaxies and AGNs from the SHELA Field shown as red diamonds with Poisson statistic error bars. The open red diamonds are the luminosity function points in bins below our 50%-completeness limit as discussed in Section 3.4, where the completeness corrections are unreliable due to the low S/N of our data. We constrain the form of the luminosity function by including fainter galaxies from Hubble elds (Finkelstein et al. 2015; open black squares) and brighter AGNs from SDSS DR7 (Akiyama et al. 2018; black \x"s). For comparison, we overplot as gray circles the g -band dropout luminosity function from the > 100 sq. deg. HSC SSP by Ono et al. (2018), which shows lower number densities and larger error bars in the regime (M < 23:5 mag) where AGN likely dominate (see Section 1). UV;i Our measured luminosity function is consistent with these works where they overlap. Our two best- tting functional forms are shown, as discussed in Section 4.2. The t with the smallest  is the sum of two DPL functions (DPL+DPL Fit; red solid line), one for the AGN component (red dash-dotted line) and one for the galaxy component (red dashed line). The second best t (DPL+Sch Fit; blue densely dotted line) is comprised of a DPL function for the AGN component (blue dotted line) and a Schechter function for the galaxy component (blue dash-dot-dotted line). The absolute value of the residuals of the two ts are shown in the inset plot for a subset of the data in units of the uncertainty in each bin. The data favors the DPL+DPL Fit suggesting the M 0 = 23:5 mag bin is dominated by star-forming galaxies, though this is dominated by the observed UV;i number densities in just a few bins, thus it is dicult to rule out a Schechter form for the galaxy component. −1 −3 Φ [Number mag Mpc ] Absolute Residual (σ) 16 Stevans et al. lensing with both functions. The Schechter (1976) func- We account for Eddington bias in our tting routine. tion has been found to t the star-forming galaxy UV Rather than directly comparing the observations to a luminosity function well across all redshifts (e.g., Steidel given model, we forward model the e ects of Eddington et al. 1999; Bouwens et al. 2007; Finkelstein et al. 2015). bias into the luminosity function model, and compare The Schechter function is described as this \convolved" model to our observations. We do this by, for each set of luminosity function parameters, re- 0:4 ln(10) (M ) = ; (12) 0:4(MM ) alizing a mock sample of galaxies for that given func- 0:4( +1)(MM ) 10 10 e tion, where each galaxy has a magnitude according to where  is the overall normalization, M is the char- the given luminosity function distribution. The magni- acteristic magnitude, and is the faint-end slope. tude of each simulated object is perturbed by an amount We consider the e ects of gravitational lensing on the drawn from a normal distribution centered on zero with shape of the star-forming galaxy UV luminosity func- a width equal to the real sample median uncertainty in tion. Gravitational lensing can distort the shape and the corresponding magnitude bin. After perturbing, we magni cation of distant sources as the paths of photons then re-bin the simulated luminosity distribution and from the source get slightly perturbed into the line of this binned luminosity function is used to calculate the sight of the observer. Lima et al. (2010) showed that for that MCMC step. This is repeated for each step. this magni cation can contribute to a bright-end ex- We ran our MCMC tting algorithm twice. In both cess where the slope of the intrinsic luminosity function runs, we simultaneously t a DPL function to the AGN is suciently steep. A magni cation distribution for a component of the luminosity function while tting the given source redshift must be estimated by tracing rays galaxy component. In the rst run, we t a Schechter through a series of lens planes derived from simulations function to the galaxy luminosity function, while in the such as the Millennium Simulation as done by Hilbert second we t a second DPL function to the galaxy lu- et al. (2007). van der Burg et al. (2010) showed that minosity function. We burn each chain for 210 steps, a Schechter function corrected for magni cation can t which allows the chain to reach convergence for all free the bright-end of the luminosity function at z = 3 better parameters, veri ed by examining the parameter distri- than a Schechter t alone. They inspect the sources that butions in independent groups of 10 steps, which cease make up the excess and nd nearby massive foreground to evolve much past 10 steps. We then measure the galaxies or groups of galaxies that could act as lenses. posterior for each parameter from the nal 510 steps. We incorporate the e ects of gravitational lensing in our Our ducial values for each parameter are then derived tting by creating a lensed Schechter function parame- as the median and 68% central con dence range from terization following the method of Ono et al. (2018) who the posterior distributions. The best tting parameters adapts the method of Wyithe et al. (2011). We also for the two ts and their corresponding  values are produce a lensed DPL function. After performing our listed in Table 4. We over plot the ts to the luminosity simultaneous tting method, which we describe in the function data in Figure 8. To calculate the  for our following paragraph, we found there to be no di erence ts, we must compare the observed data to the lumi- in the best tting parameters of the ts including and nosity function ts after forward modeling the e ects of excluding the e ects of lensing. This is consistent with Eddington bias into our luminosity function ts just as the work of Ono et al. (2018) who found that taking we did during the tting process. The absolute residuals into account the e ects of lensing improves the galaxy of our \convolved" ts and the observed data are shown luminosity function t at z > 4 and not at z = 4 where in the insert in Figure 8. a DPL t is preferred. Therefore we do not consider the These tting results show that our data prefers the lensed parameterizations further. function that is a sum of two DPL functions, one for We employ a Markov Chain Monte Carlo (MCMC) the AGN component and one for the galaxy compo- method to de ne the posteriors on our luminosity func- nent, henceforth the DPL+DPL Fit. The DPL+DPL tion parameterizations. We do this using an IDL imple- Fit has a  = 42 over the entire range considered, mentation of the ane-invariant sampler (Goodman & 29 < M 0 < 17. The t that is a sum of a UV;i Weare 2010) to sample the posterior, which is similar DPL AGN component and a Schechter galaxy compo- in production to the emcee package (Foreman-Mackey nent, henceforth the DPL+Schetcher Fit, has a  = 71. et al. 2013). Each of the 500 walkers was initialized We investigate which t is preferred by the data using by choosing a starting position with parameters deter- the Bayesian information criterion for the two ts (Lid- mined by-eye to exhibit a good t, perturbed according dle 2007). The di erence in the Bayesian information to a normal distribution. We do not assume a prior for criterion value for our two ts can be de ned as: any of our free parameters. Galaxy and AGN Luminosity Functions at z = 4 with SHELA 17 Table 4. Fit Parameters for Luminosity Functions AGN Component Galaxy Component Fit Name Log  M Log  M (mag) [Faint End] [Bright End] (mag) [Faint End] [Bright End] +0:21 +0:4 +0:30 +0:21 +0:09 +0:16 DPL+DPL -7.32 -26.5 -1.49 -3.65 -3.12 -20.8 -1.710.08 -3.800.10 42 0:18 0:3 0:21 0:24 0:10 0:15 +0:58 +1:1 +0:18 +0:68 DPL+Sch -7.48 -26.7 -2.08 -3.66 -3.250.06 -21.30.06 -1.810.05  71 0:34 0:4 0:11 0:34 3 1 Note| 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). 1.0 2 2 SDSS Sample BIC =   + (k k ) ln N; (13) 2 1 2 1 Selected 0.8 2 2 2 where  is the  for the DPL+Schechter Fit,  is 2 2 0.6 the  for the DPL+DPL Fit, k and k are the num- 10 2 1 0.4 ber of t parameters in the DPL+Schechter Fit and the DPL+DPL Fit, respectively, and N is the number of 0.2 data points used during tting. We nd a BIC = 26 0 0.0 which suggests the DPL+DPL Fit is very strongly pre- 3.5 4.0 4.5 5.0 3.0 3.5 4.0 4.5 5.0 z z spec spec ferred over the DPL+Schechter Fit as BIC exceeds a signi cance value of 10 (Liddle 2007). Upon compar- ing the the ts' residuals we see the data points in the Figure 9. Left: The SDSS spectroscopic redshift distribu- tion of AGNs in SHELA (blue) overplotted with the same magnitude bins at M 0 = 24 mag and M 0 = 23 UV;i UV;i distribution for the sub-sample selected by our selection pro- mag are driving the preference for the DPL+DPL Fit. If cedure (red). Right: The di erential completeness fraction we exclude these two bins, the DPL+DPL Fit  = 41, of AGNs with SDSS spectroscopic redshifts (blue) and the the DPL+Schechter Fit  = 58, and the BIC = 14 expected di erential completeness fraction for objects with which indicates the DPL+DPL Fit is still strongly sta- a comparable M (green). Our completeness to AGNs is UV;i tistically preferred to the DPL+Schechter Fit. However, similar to what we expect from our simulations of galaxies given that the di erence between the ts is driven by except at z = 3:5 where we are less complete to AGNs due to 0 0 signi cant u and g ux driving the redshift probability dis- just a few data points, we do not believe we can rmly tribution functions to peak at redshifts below our selection rule out a Schechter form for the star-forming galaxy window (see Section 4.3 for details). component. 4.3. A Sample of Spectroscopically Con rmed AGNs Each match required a separation of less than 0.4" and the SDSS spectra to be un agged (i.e., ZWARNING = Given our method of tting the luminosity function 0). We found zero spectra with an SDSS classi cation of with a component for AGNs, we explored the SDSS Galaxy and 53 classi ed as AGN with spectroscopic red- spectral catalog for spectroscopically con rmed AGNs shifts (z ) greater than 3.2. The distribution of SDSS spec in SHELA and considered the e ectiveness of our se- z for this sample is shown in the left panel of Figure 9 spec lection procedure at recovering AGNs. This is crucial in blue. Of the 32 AGNs with 3:4 < z < 4:6, 23 were spec as our photometric redshift selection process did not in- selected by our z selection process, with seven of the phot clude templates dominated by AGN emission, though nine missed AGNs having 3:4 < z < 3:5. The z spec spec the strong Lyman break inherent in these sources should distribution for the AGN subsample selected by our se- still yield an accurate redshift. To con rm this, we lection procedure is also shown in the left panel of Figure queried the spectral catalog from SDSS (DR13; Albareti 9 in red. We used these samples to compute our di eren- et al. 2017) using the SDSS CasJobs website and cross- tial completeness of AGNs and compared it to our sim- matched the results with our entire DECam catalog. ulated completeness in the right panel of Figure 9. We found our completeness of spectroscopically con rmed http://skyserver.sdss.org/casjobs/ AGNs is consistent with our simulated completeness ex- Number Completeness 18 Stevans et al. cept in the 3:4 < z < 3:6 bin where we recovered only data, the contamination rate was signi cantly higher, two of the nine spectroscopically con rmed AGNs when although we acknowledge there are multiple di erences our simulations predicted we should recover 78. This in the selection techniques between these studies. di erence could be due to small number statistics. We At the bright end, the luminosity function from investigated why the seven spectroscopically con rmed Ono et al. (2018) extends to magnitudes brighter than AGNs in the 3:4 < z < 3:6 bin were not selected by our M 0 < 22:5 and falls between the bright end of UV;i method and found that these sources had photometry, our DPL galaxy component and our Schechter galaxy 0 0 particularly the u and g bands, preferred by galax- component. Ono et al. (2018) nd their luminosity ies templates at lower redshift in our z - tting code function is best t by a DPL with a steeper bright-end phot 0 0 EAZY. We attribute bright u and g uxes to the larger slope ( = 4:33) than our galaxy DPL component far-UV continuum levels of AGNs or strong Lyman- ( = 3:80). We also include the z = 4 galaxy lumi- emission as compared to non-AGN galaxies, which nosity function from the 4 deg ALHAMBRA survey by would weaken the Lyman- break in these sources. This Viironen et al. (2017) who used a z PDF analysis phot could imply signi cant leaking Lyman-continuum radia- to create a luminosity function marginalizing over both tion from these AGNs, which has implications on reion- redshift and magnitude uncertainties. Viironen et al. ization (e.g., Smith et al. 2016). The reliability of our (2017) nd volume densities at the bright-end larger selection procedure to recover AGNs across the major- than existing luminosity functions. In fact, their lu- ity of our redshift range of interest and the fact that our minosity function follows a Schechter function with a luminosity function is consistent with the AGN luminos- normalization o set of 0.5 dex above the Schechter ity function from SDSS suggests our incompleteness to form from Finkelstein et al. (2015). The cause of this AGNs is not substantial. di erence is unclear. 5. DISCUSSION 5.1.1. Comparison to Semi-Analytic Models 5.1. Comparison to z=4 Galaxy Studies Figure 11 shows the predicted z = 4 UV luminosity We compare our derived ts to the star-forming galaxy functions from Yung et al. (2018). These predictions luminosity function to measurements from the litera- come from a set of semi-analytic models (SAMs), which ture in Figure 10. At magnitudes fainter than M 0 > contain the same physical processes as the models pre- UV;i 22, the DPL galaxy component of the data-preferred sented in Somerville et al. (2015), but have been up- DPL+DPL Fit is very similar to the Schechter compo- dated and recalibrated to the Planck 2016 Cosmological nent of the DPL+Schechter Fit. These results are in parameters. We note that while these models include strong agreement with the luminosity function from the black hole accretion and the e ects of AGN feedback, CFHT Deep Legacy Survey by van der Burg et al. (2010) we do not examine the contribution to the UV luminos- at all magnitudes where they overlap. They are also ity function from black hole accretion here, and focus in strong agreement with the luminosity function from instead on the contribution from star formation. Their Hubble legacy survey data by Finkelstein et al. (2015) ducial model with dust (solid black line) assumes that which is included in our tting. The luminosity func- the molecular gas depletion time is shorter at higher gas tion from the Hubble legacy survey data by Bouwens densities (as motivated by observations), leading to an et al. (2015) and the luminosity function from Ono et al. e ective redshift dependence as high redshift galaxies (2018) using Subaru Hyper Suprime-Cam (HSC) data are more compact and have denser gas on average. On across 100 deg in the HSC Subaru strategic program a SFR surface density versus gas surface density plot, (SSP) are generally consistent with the results presented this model would show a steeper dependence of star for- here. The HSC SSP luminosity function and that from mation rate density on gas surface density than the clas- Bouwens et al. (2015) have volume densities of 2x sical Kennicutt-Schmidt relationship (e.g., Kennicutt & larger than the work by van der Burg et al. (2010) Evans 2012), above a critical gas surface density (for and Finkelstein et al. (2015) at magnitudes fainter than details see Somerville et al. 2015). We also show their M > 21. This factor is larger than the 10% un- model with a xed molecular gas depletion time, similar UV;i certainty expected due to cosmic variance in the Hubble to that used in Somerville et al. (2015), as seen in local elds, which cover 50x less area than the HSC SSP spiral galaxies (e.g., Bigiel et al. 2008; Leroy et al. 2008). (Finkelstein et al. 2015). We do not know the cause of The Yung et al. (2018) dust models assume that the this di erence, though we note that the HSC SSP selec- V-band dust optical depth is proportional to the cold gas tion was done with optical imaging only, and we found metallicity times the cold gas surface density. The UV in our study that without the addition of the IRAC attenuation is then computed using a Calzetti attenua- Galaxy and AGN Luminosity Functions at z = 4 with SHELA 19 -2 This Work -2 Finkelstein+15 -3 -3 -4 -4 -5 -5 Galaxy DPL Component Galaxy DPL Component -6 Galaxy Schechter Component -6 Galaxy Schechter Component Finkelstein+15 SAMs w/dust Bouwens+15 -7 -7 SAMs w/dust & fixed H t 10 2 dep 10 van der Burg+10 SAMs w/no dust Ono+18 -8 SAMs w/no dust & fixed H t -8 2 dep Viironen+17 10 10 24 23 22 21 20 19 18 17 25 24 23 22 21 20 19 18 M 0 M 0 UV, i UV, i Figure 10. The star-forming galaxies components of our Figure 11. The star-forming galaxy component of our ts ts to the total z = 4 rest-frame UV luminosity function to the total z = 4 rest-frame UV luminosity function com- compared to the data from other star-forming galaxy stud- pared to predictions from SAMs (Yung et al. 2018). Our ies See the legend insert for the list of compared works. Our measured total luminosity function (including star-forming DPL galaxy luminosity function is in agreement with lumi- galaxies and AGNs) is shown as red diamonds with Poisson nosity functions from the literature around the knee and to statistic error bars. The luminosity function from Finkel- fainter magnitudes, but has the shallowest bright-end slope stein et al. (2015) is shown as open black squares. At the ( = 3:80). faint end the luminosity functions from the SAMs with an evolving H gas depletion time (solid black line) and a xed H gas depletion time (dashed black line) have a higher nor- tion curve and a \slab" model (for details see Somerville malization than the observed luminosity function suggesting et al. 2015). The normalization of the dust optical depth stellar feedback in the models is too weak. At the bright (physically equivalent to the dust-to-metal ratio) is al- end the model with a xed molecular depletion timescale is lowed to vary as a function of redshift, and was adjusted robustly ruled out by either of our parameterized ts, sug- gesting that star formation scales with molecular gas surface to t the bright end of the luminosity from the previ- density, and thus is more ecient at z = 4 than today. ously published compilation of UV luminosity observa- tions from Finkelstein (2016). (It was not adjusted to t the new observations presented here). At the bright end, the Yung et al. (2018) model with At the faint end, the model predictions are insensitive dust is consistent with both of our ts to M > UV;i to the assumed star formation eciency, and mainly re- 22.5, lying closer to the Schechter t at brighter mag- ect the treatment of out ows driven by stars and su- nitudes. Interestingly, at these bright magnitudes both pernovae. The models have a higher normalization than the Schechter and DPL t rule out at high signi cance the observed luminosity function (which at these mag- the model with dust and xed molecular gas depletion nitudes comes from the CANDELS dataset), although time, indicating that star formation must scale non- the faint-end slope is similar. This could be caused by linearly with molecular gas surface density (or some re- stellar feedback in the models being too weak, as these lated quantity which evolves with redshift). This implies models were tuned to match the z = 0 stellar mass func- that star formation is more ecient at z = 4 than to- tion (Somerville et al. 2015). This suggests that stellar day. This is of course dependant on the dust model feedback is stronger/more ecient at z = 4 (i.e., mass in these simulations { if bright galaxies had no dust, loading rates are higher, or re-infall of ejected gas is then these model predictions (which include dust; mod- slower) than at z = 0 (also see White et al. 2015). This els without dust are shown for comparison) may not be was also seen by Song et al. (2016) when comparing the accurate. However, there is a variety of evidence that z = 4 stellar mass function to a number of similar mod- bright/massive UV-selected galaxies at these redshifts els { the observed stellar mass function also had a lower do contain non-negligible amounts of dust (e.g., Finkel- normalization at z = 4; as the stellar mass function stein et al. 2012; Bouwens et al. 2014), and there is a steepened from z = 4 to 8, this discrepancy weakened, not-insigni cant population of extremely massive dusty hence their conclusion of a weakening impact of feedback star-forming galaxies already in place (e.g., Casey et al. on the faint end with increasing redshift. 2014). Finally, Finkelstein et al. (2015) compared a sim- −1 −3 Φ [Number mag Mpc ] −1 −3 Φ [Number mag Mpc ] 20 Stevans et al. +0:18 ilar set of models to the CANDELS-only UV luminos- 2:08 ) predicts volume densities in agreement with 0:11 ity functions, nding that models with a di use dust Glikman et al. (2011) and Giallongo et al. (2015), who component only (e.g., no birth cloud) provided the best predict a signi cantly higher volume density of faint match to the data (albeit at fainter magnitudes than AGNs than the other studies. We note that while these considered here). studies have small numbers of AGN per magnitude bin, which may make the samples more sensitive to cosmic variance and false positives, Glikman et al. (2011) se- 5.2. Comparison to z = 4 AGN Studies lects candidates from a relatively large survey area of We compare our derived AGN luminosity function ts 3.76 deg and observes broad emission lines, indicative to measurements from the literature in Figure 12. At the of quasars, in every candidate spectra. Furthermore, bright end, our ts are consistent with previous studies Giallongo et al. (2015) requires candidates to have a where they overlap (27:5 < M < 25:5). The pre- UV;i signi cant X-ray detection in Chandra imaging. vious studies include a study by Richards et al. (2006) In summary, our observed z = 4 UV luminosity func- who used the z = 4 AGN SDSS DR3 sample and a tion is best t by the DPL+DPL Fit, but both the study by Akiyama et al. (2018) who selected z = 4 AGNs DPL+DPL and the DPL+Schechter ts show agreement from SDSS DR7. We convert the magnitudes M at i(z=2) with existing AGN luminosity function studies at the z = 3:75 from Richards et al. (2006) to M 0 at z = 3:8 UV;i bright end and around the knee. At the faint end the by adding an o set of 1.486 mag (Richards et al. 2006). DPL+DPL predicts smaller volume densities of AGNs At the faint end, we compare our ts to results from than other studies while the DPL+Schechter predicts studies that derive AGN luminosity functions from spec- among the largest volume densities at the faintest mag- troscopic observations of candidates selected via color nitudes. and size criteria (Glikman et al. 2011, Ikeda et al. 2011, and Niida et al. 2016). In addition, we include studies that rely on a z selection using deep multi wave- phot 5.3. Comparing Predictions of Our Two Fits - Is the length photometry in the COSMOS eld (Masters et al. UV luminosity function a double power law? 2012) and the CANDELS GOODS-S eld (Giallongo et al. 2015) where Giallongo et al. (2015) had the ad- Our two ts di er in two ways: 1.) The functional ditional criteria of requiring an an X-ray detection of form used to t the star-forming galaxy component is 17 2 1 F > 1:5 10 erg cm s in deep 0.5-2 keV Chan- a DPL in the DPL+DPL Fit and a Schechter in the dra imaging. The average redshift of these samples is DPL+Schechter Fit, and 2.) The component that ac- z = 4, slightly higher than our sample. For a consistent counts for the excess over an extrapolated Schechter at comparison with other works we scale the luminosity the bright end of existing star-forming galaxy luminos- functions of Glikman et al. (2011), Ikeda et al. (2011), ity functions. The DPL+Schechter accounts for the ex- Niida et al. (2016), and Masters et al. (2012) up by a fac- cess with a steeper faint end of the DPL AGN compo- 6:9 tor of 1.3 using the redshift evolution function (1+z) nent while the DPL+DPL Fit accounts for the major- from Richards et al. (2006). The sample used by (Gial- ity of the excess with the bright end of the DPL star- longo et al. 2015) had a redshift range of 4 < z < 4:5, so forming galaxy component. Thus the ts predict dra- we the scale the luminosity function by a factor of 1.8. matically di erent compositions of sources making up Finally, we consider the luminosity function of Akiyama the bin (M 0 = 23:5 mag) at the center of the ex- UV;i et al. (2018) derived from the Hyper Suprime-Cam Wide cess. In this bin the DPL+Schechter Fit predicts AGNs Survey optical photometry. outnumber galaxies by a 17:1 ratio (an AGN fraction The faint end of the AGN DPL component in our of 94%), while the DPL+DPL Fit predicts galaxies DPL+DPL Fit predicts volume densities at 26 < to outnumber AGNs by a 4.3:1 ratio (an AGN fraction M 0 < 23:5 about 0.3 dex lower than those found by of 18%). A simple experiment to test these predic- UV;i Richards et al. (2006), Ikeda et al. (2011), Masters et al. tions would be to use ground-based optical spectroscopy (2012), Niida et al. (2016), Akiyama et al. (2018). How- to follow-up a fraction of the 298 high-z candidates we ever, the luminosity function of Akiyama et al. (2018) nd in the M = 23:5 bin and count the fraction UV;i attens and falls towards our t at M 0  22. Our of spectra exhibiting broad emission lines and/or highly UV;i +0:30 faint-end slope ( = 1:49 ) is in agreement with ionized lines (e.g., N V, He II, C IV, Ne V) indicating ac- 0:21 the values found by these studies. The cause of the 0.3 cretion onto supermassive black holes. This experiment dex di erence is unclear. would also provide an independent measurement of our In the case of the AGN DPL component in our contamination fraction which we estimate via simula- DPL+Schechter Fit, the steeper faint-end slope ( = tions (Section 3.3). The measured fraction of AGNs can Galaxy and AGN Luminosity Functions at z = 4 with SHELA 21 -4 Richards+06 SDSS DR7 (Akiyama+18) HSC (Akiyama+18) -5 Glikman+11 -6 -7 -8 AGN DPL with galaxy DPL AGN DPL with galaxy Schechter Ikeda+11 -9 Niida+16 Masters+12 -10 Giallongo+15 28 26 24 22 20 UV, i Figure 12. The AGN components of our ts to the total z = 4 rest-frame UV luminosity function compared to the data from other AGN studies. See the legend insert for the list of compared works. The AGN luminosity function from the DPL+DPL Fit has number densities similar to existing luminosity function measurements at M < 25 and predicts relatively low number UV;i densities at fainter magnitudes. The DPL+Schechter Fit has an AGN faint-end slope that predicts number densities at the larger end of the range previously published. Table 5. UV Luminosity Densities provide strong empirical proof in favor of the DPL+DPL Fit or the DPL+Schechter Fit. AGN Component Galaxy Component One drawback of tting the total UV luminosity func- tion is the unknown contribution of composite objects. The total UV luminosity function likely contains a pop- Fit Name 1500 912 1500 912 ulation of composite objects with both AGN activity +0:7 +0:4 +6 DPL+DPL 2.0 1.2 184 0:4 0:3 7 and star-formation contributing to their observed UV +4 +2:6 +6 DPL+Sch 10 5.8 187 3 1:8 6 ux. This population may cause functions like a DPL 24 1 1 3 Note|All values are in units of 10 ergs s Hz Mpc . or Schechter function, which have been used to t AGN- only and star-forming-galaxy-only luminosity functions, to be poor ts to the total UV luminosity function. Spectroscopic follow-up as described in this section may aid in elucidating the frequency and impact of composite that of AGNs by a factor of 90. In the case of the objects. DPL+Schechter Fit where the galaxy component of the UV luminosity function is instead represented by a DPL 5.4. Rest-Frame UV Emissivities function, the galaxy  is greater than the AGN 1500 1500 by a factor of 19. In either case, galaxies are the dom- Here we calculate the rest-frame UV emissivities (also inant non-ionizing UV emitting population, though if known as speci c luminosity densities) and compare the AGNs do have a steeper faint-end slope (as would need output from AGNs and galaxies. We calculate the rest- to be the case if galaxies follow a Schechter form), then frame UV emissitvity of AGNs and galaxies by integrat- AGNs are non-negligible. ing each luminosity function from 30 < M 0 < 17. UV;i We convert the AGN  to a UV luminosity den- Results are shown in Table 5. In the case of DPL+DPL sity at 912 A ( ) and compare our results to other where the galaxy component of the UV luminosity func- studies by reproducing Figure 1 from Madau & Haardt tion is represented by a DPL function and the AGN (2015) in Figure 13. To convert  to  we as- 1500 912 component is represented by a DPL, galaxies produce a sume an AGN spectral shape of a DPL with a slope UV luminosity density at 1500 A ( ) greater than −1 −3 Φ [Number mag Mpc ] 22 Stevans et al. faint AGNs at z  4 Glikman et al. 2011; Giallongo et al. 2015. This would imply the AGN population could alone contribute enough hydrogen ionizing emission required to keep the universe ionized at z  4 and suggests AGNs may have played a signi cant role at early times. This scenario can be further constrained by the spectroscopic experiment proposed in the preceding sub-section. DPL+Sch Fit 6. CONCLUSIONS DPL+DPL Fit In this study, we measure the bright end of the rest- Madau & Haardt 15 frame UV luminosity function of z = 4 star-forming Hopkins+07 galaxies and the faint end of the rest-frame UV lumi- 0 1 2 3 4 5 6 7 nosity function of z = 4 AGNs. We use nine photo- 0 0 0 0 0 metric bands (u g r i z from DECam, J and K from VISTA, and 3:6 and 4:5 m from Spitzer /IRAC) cov- ering the wide area (18 deg ) SHELA Field to select Figure 13. The AGN hydrogen ionizing emissivity from 3,740 candidate z  4 galaxies via a photometric red- this work and others. The AGN emissivity predicted by the shift selection procedure. From simulations, we nd a DPL+DPL Fit and the DPL+Schechter Fit are represented relatively low contamination rate of interloping low-z as red triangles. The range they span is marked by the thick red line. The values from other works are inferred from Fig galaxies and galactic stars of 20% near our complete- 1 of Madau & Haardt (2015). The original sources of each ness limit (m  23) due in large part to the inclusion dataset are as follows: Schulze et al. (2009) (cyan pentagon), of IRAC photometry. Palanque-Delabrouille et al. (2013) (orange triangles), Bon- Our conclusions are as follows: giorno et al. (2007) (magenta circles), Masters et al. (2012) We combine our candidate sample with a sample of (black pentagons), Glikman et al. (2011) (blue square), and Giallongo et al. (2015) (green squares). The solid blue line bright AGNs from SDSS and fainter galaxies from is the functional form derived by Madau & Haardt (2015) deep Hubble imaging (including the HUDF and to coincide with the plotted observation from the literature. CANDELS) to produce a rest-frame UV luminos- The dotted green line shows the LyC AGN emissivity from ity function that spans the range -29< M 0 <- UV;i Hopkins et al. (2007). The  from our DPL+DPL Fit 17. This range contains both AGNs and star- is below the line by Hopkins et al. (2007) indicating AGNs forming galaxies several magnitudes above and contribute only a small fraction of the total ionizing back- below their respective characteristic luminosities, ground at z = 4, while the  from the DPL+Schechter t suggests AGNs would contribute signi cantly to the total thus we implement a tting procedure that simul- ionizing background. taneously ts the AGN luminosity function and the star-forming galaxy luminosity function with of = 1:41 (Shull et al. 2012) between 1000 A and independent functions. This simpli es the source selection process by not requiring a step for classi- 1500 A and a slope of = 0:83 (Stevans et al. 2014) between 912 A and 1000 A . We assume an AGN ioniz- fying objects as either an AGN or galaxy, which is commonly done with morphological criteria. We ing escape fraction of unity as is found for most bright nd the data is best t by our DPL+DPL Fit AGNs (Giallongo et al. 2015). Results are shown in Ta- which is a combination of a DPL function for ble 5. If we directly follow the work of Giallongo et al. the AGN component and a DPL function for the (2015) and assume a spectral shape with = 1:57 galaxy component. The DPL+DPL t is preferred (Telfer et al. 2002) between 1200 A and 1500 A and over the DPL+Shechter Fit, which is a combina- a slope of = 0:44 (Vanden Berk et al. 2001) be- tween 912 A and 1200 A , we nd  values that are tion of a DPL function for the AGN component and a Schechter function for the galaxy compo- 10% smaller. As shown in Figure 13, our two ts pre- dict values of  straddle existing values from observa- nent, and this excess over Schechter cannot be explained by the e ects of gravitational lensing. tions. The  predicted by our preferred DPL+DPL Fit (upward-pointed red triangle) falls below the line by We note that we cannot signi cantly rule out a Schechter form. Hopkins et al. (2007) and is too small to keep the IGM ionized at z  4. On the other hand, the DPL+Sch We compare our measured luminosity functions Fit predicts a  (downward-pointed red triangle) that to the literature and nd our DPL galaxy lumi- falls near the points from studies that found numerous nosity function is in agreement with luminosity −1 −1 −3 ρ [erg s Hz Mpc ] 912 Galaxy and AGN Luminosity Functions at z = 4 with SHELA 23 functions from the literature around the knee and Future work is needed to con rm the shape of the to fainter magnitudes while having the shallow- star-forming galaxy and AGN luminosity functions, es- est bright-end slope. The AGN luminosity func- pecially where they intersect. We discuss a simple ex- tion from the DPL+DPL Fit has number densities periment to measure the relative number densities of similar to existing luminosity functions at magni- AGNs and galaxy at luminosities where the respective tudes up to M = 25:5 mag while under pre- luminosity functions intersect. Spectroscopic follow-up UV;i dicting number densities by  0:3 dex at fainter of a sample of our z = 4 candidates in the M -23.5 UV;i magnitudes. The DPL+Schechter Fit has an AGN bin where our two ts predict di erent AGNs to galaxies faint-end slope that is among the steepest values ratios is underway. Imaging from space-based telescopes published. The shape of the galaxy bright end is such as HST or JWST would facilitate a robust mor- consistent with model predictions where star for- phological classi cation of our bright candidates. Other mation is more ecient at higher redshift due to possibilities for distinguishing AGNs including taking increased gas densities. deeper X-ray imaging in the eld and using JWST to measure mid-IR SEDs (Kirkpatrick et al. 2012). The authors acknowledge Raquel Martinez, Karl Geb- We measure  by integrating the rest-frame hardt, and Eric J. Gawiser for the interesting discussions UV luminosity function ts and nd that galax- we had and their suggestions which improved the quality ies dominate the production of non-ionizing ux of this research. We thank Neal J. Evans and William P. at z = 4 for both possible ts. Speci cally, galax- Bowman for useful comments on the draft of this paper. ies produce a factor of 90 more non-ionizing UV M. L. S. and S. L. F. acknowledge support from the Uni- output than AGNs according to the DPL+DPL versity of Texas at Austin, the NASA Astrophysics and Fit, while the DPL+Schechter Fit predicts galax- Data Analysis Program through grant NNX16AN46G, ies produce a factor of 19 more non-ionizing UV and the National Science Foundation AAG Award AST- output than AGNs. 1614798. The work of C. P. And L. K. is supported by the National Science Foundation through grants AST 413317, and 1614668. The Institute for Gravitation and We convert the AGN  to  and nd AGNs 1500 912 the Cosmos is supported by the Eberly College of Sci- do not produce enough ionizing radiation to keep ence and the Oce of the Senior Vice President for Re- the universe ionized at z = 4 by themselves if the search at the Pennsylvania State University. S. J., S. AGN is truly represent by the DPL component S. and J. F. acknowledge support from the University of in our DPL+DPL Fit. This suggests AGNs are Texas at Austin and NSF grants NSF AST-1614798 and not the dominant contributor to cosmic reioniza- NSF AST-1413652. R. S. S. and A. Y. thank the Downs- tion at earlier times. On the other hand, if the brough family for their generous support, and gratefully DPL+Schechter Fit is true, AGNs could alone pro- acknowledge funding from the Simons Foundation. duce the  needed to maintain the ionized state of the universe at z = 4 and perhaps at earlier times. Facilities: CTIO,SST(IRAC) REFERENCES Akiyama, M., He, W., Ikeda, H., et al. 2018, PASJ, 70, S34, Bertin, E., Mellier, Y., Radovich, M., et al. 2002, in Astronomical Society of the Paci c Conference Series, doi: 10.1093/pasj/psx091 Vol. 281, Astronomical Data Analysis Software and Albareti, F. D., Allende Prieto, C., Almeida, A., et al. 2017, Systems XI, ed. D. A. Bohlender, D. Durand, & T. H. Handley, 228 ApJS, 233, 25, doi: 10.3847/1538-4365/aa8992 Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846, Beers, T. C., Flynn, K., & Gebhardt, K. 1990, doi: 10.1088/0004-6256/136/6/2846 Astronomical Journal (ISSN 0004-6256), 100, 32 Bouwens, R. J., Illingworth, G. D., Franx, M., & Ford, H. 2007, The Astrophysical Journal, 670, 928 Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. doi: 10.1051/aas:1996164 2014, ApJ, 793, 115, doi: 10.1088/0004-637X/793/2/115 24 Stevans et al. |. 2015, ApJ, 803, 34, doi: 10.1088/0004-637X/803/1/34 Kirkpatrick, A., Pope, A., Alexander, D. M., et al. 2012, ApJ, 759, 139, doi: 10.1088/0004-637X/759/2/139 Boyle, B. J., Shanks, T., Croom, S. M., et al. 2000, Monthly Notices of the Royal Astronomical Society, 317, 1014 Kurucz, R. L. 1979, ApJS, 40, 1, doi: 10.1086/190589 Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, LaMassa, S. M., Urry, C. M., Cappelluti, N., et al. 2016, The Astrophysical Journal, 686, 1503 The Astrophysical Journal, 817, 172 Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, Lang, D., Hogg, D. W., & Mykytyn, D. 2016a, The Tractor: Astrophysical Journal, 345, 245 Probabilistic astronomical source detection and measurement, Astrophysics Source Code Library. Casey, C. M., Scoville, N. Z., Sanders, D. B., et al. 2014, ApJ, 796, 95, doi: 10.1088/0004-637X/796/2/95 http://ascl.net/1604.008 Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833, Lang, D., Hogg, D. W., & Schlegel, D. J. 2016b, AJ, 151, doi: 10.1088/0004-637X/712/2/833 36, doi: 10.3847/0004-6256/151/2/36 Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486, Le F evre, O., Cassata, P., Cucciati, O., et al. 2013, doi: 10.1088/0004-637X/699/1/486 Astronomy & Astrophysics, 559, A14 Lehmer, B. D., Xue, Y. Q., Brandt, W. N., et al. 2012, The Croom, S. M., Richards, G. T., Shanks, T., et al. 2009, Monthly Notices of the Royal Astronomical Society, 399, Astrophysical Journal, 752, 46 1755 Leroy, A. K., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2782, doi: 10.1088/0004-6256/136/6/2782 Cucciati, O., Tresse, L., Ilbert, O., et al. 2012, Astronomy & Astrophysics, 539, A31 Liddle, A. R. 2007, MNRAS, 377, L74, Eisenstein, D. J., Weinberg, D. H., Agol, E., et al. 2011, doi: 10.1111/j.1745-3933.2007.00306.x AJ, 142, 72, doi: 10.1088/0004-6256/142/3/72 Lima, M., Jain, B., & Devlin, M. 2010, MNRAS, 406, 2352, Finkelstein, S. L. 2016, PASA, 33, e037, doi: 10.1111/j.1365-2966.2010.16884.x doi: 10.1017/pasa.2016.26 Madau, P. 1995, ApJ, 441, 18, doi: 10.1086/175332 Finkelstein, S. L., Papovich, C., Salmon, B., et al. 2012, Madau, P., & Haardt, F. 2015, The Astrophysical Journal ApJ, 756, 164, doi: 10.1088/0004-637X/756/2/164 Letters, 813, L8 Finkelstein, S. L., Song, M., Behroozi, P., et al. 2015, Marchesini, D., van Dokkum, P. G., F orster Schreiber, ArXiv e-prints. https://arxiv.org/abs/1504.00005 N. M., et al. 2009, ApJ, 701, 1765, Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, doi: 10.1088/0004-637X/701/2/1765 J. 2013, PASP, 125, 306, doi: 10.1086/670067 Masters, D., Capak, P., Salvato, M., et al. 2012, The Astrophysical Journal, 755, 169 Gawiser, E., van Dokkum, P. G., Herrera, D., et al. 2006, ApJS, 162, 1, doi: 10.1086/497644 McGreer, I. D., Fan, X., Jiang, L., & Cai, Z. 2017, arXiv.org Geach, J. E., Lin, Y. T., Makler, M., et al. 2017, The Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, Astrophysical Journal Supplement Series, 231, 7 777, 18, doi: 10.1088/0004-637X/777/1/18 Giallongo, E., Grazian, A., Fiore, F., et al. 2015, Niida, M., Nagao, T., Ikeda, H., et al. 2016, The Astronomy & Astrophysics, 578, A83 Astrophysical Journal, 832, 208 Glikman, E., Djorgovski, S. G., Stern, D., et al. 2011, The Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713, Astrophysical Journal Letters, 728, L26 doi: 10.1086/160817 Goodman, J., & Weare, J. 2010, Communications in Ono, Y., Ouchi, M., Harikane, Y., et al. 2018, PASJ, 70, Applied Mathematics and Computational Science, Vol. 5, S10, doi: 10.1093/pasj/psx103 No. 1, p. 65-80, 2010, 5, 65, Papovich, C., Shipley, H. V., Mehrtens, N., et al. 2016, doi: 10.2140/camcos.2010.5.65 ApJS, 224, 28, doi: 10.3847/0067-0049/224/2/28 Hilbert, S., White, S. D. M., Hartlap, J., & Schneider, P. Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, 2007, MNRAS, 382, 121, Journal of Machine Learning Research, 12, 2825 doi: 10.1111/j.1365-2966.2007.12391.x Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. Hopkins, P. F., Richards, G. T., & Hernquist, L. 2007, The 2014, A&A, 571, A16, doi: 10.1051/0004-6361/201321591 Astrophysical Journal, 654, 731 Quadri, R., Marchesini, D., van Dokkum, P., et al. 2006, Ikeda, H., Nagao, T., Matsuoka, K., et al. 2011, The arXiv.org, 1103 Astrophysical Journal Letters, 728, L25 Richards, G. T., Strauss, M. A., Fan, X., et al. 2006, Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531, arXiv.org, 2766 doi: 10.1146/annurev-astro-081811-125610 Schechter, P. 1976, Astrophysical Journal, 203, 297 Galaxy and AGN Luminosity Functions at z = 4 with SHELA 25 Schla y, E. F., & Finkbeiner, D. P. 2011, The Szalay, A. S., Connolly, A. J., & Szokoly, G. P. 1999, The Astronomical Journal, 117, 68 Astrophysical Journal, 737, 103 Telfer, R. C., Zheng, W., Kriss, G. A., & Davidsen, A. F. Schneider, D. P., Richards, G. T., Hall, P. B., et al. 2010, 2002, The Astrophysical Journal, 565, 773 The Astronomical Journal, 139, 2360 van der Burg, R. F. J., Hildebrandt, H., & Erben, T. 2010, Shull, J. M., Stevans, M., & Danforth, C. W. 2012, The A&A, 523, A74, doi: 10.1051/0004-6361/200913812 Astrophysical Journal, 752, 162 Vanden Berk, D. E., Richards, G. T., Bauer, A., et al. 2001, Smith, B. M., Windhorst, R. A., Jansen, R. A., et al. 2016, The Astronomical Journal, 122, 549 ArXiv e-prints. https://arxiv.org/abs/1602.01555 Viironen, K., L opez-Sanjuan, C., Hern andez-Monteagudo, C., et al. 2017, arXiv.org, arXiv:1712.01028 Somerville, R. S., Popping, G., & Trager, S. C. 2015, White, C. E., Somerville, R. S., & Ferguson, H. C. 2015, MNRAS, 453, 4337, doi: 10.1093/mnras/stv1877 ApJ, 799, 201, doi: 10.1088/0004-637X/799/2/201 Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, Wyithe, J. S. B., Yan, H., Windhorst, R. A., & Mao, S. ApJ, 825, 5, doi: 10.3847/0004-637X/825/1/5 2011, Nature, 469, 181, doi: 10.1038/nature09619 Steidel, C. C., Adelberger, K. L., Giavalisco, M., Dickinson, Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., M., & Pettini, M. 1999, ApJ, 519, 1, doi: 10.1086/307363 Popping, G., & Dav e, R. 2018, ArXiv e-prints. Steidel, C. C., Giavalisco, M., Dickinson, M., & Adelberger, https://arxiv.org/abs/1803.09761 Zacharias, N., Monet, D. G., Levine, S. E., et al. 2004, in K. L. 1996, AJ, 112, 352, doi: 10.1086/118019 Bulletin of the American Astronomical Society, Vol. 36, Stevans, M. L., Shull, J. M., Danforth, C. W., & Tilton, American Astronomical Society Meeting Abstracts, 1418 E. M. 2014, The Astrophysical Journal, 794, 75

Journal

AstrophysicsarXiv (Cornell University)

Published: Jun 13, 2018

There are no references for this article.