EPOCHS. II. The Ultraviolet Luminosity Function from 7.5 < z < 13.5 Using 180 arcmin2 of Deep, Blank Fields from the PEARLS Survey and Public JWST Data

We present an analysis of the ultraviolet luminosity function (UV LF) and star formation rate density of distant galaxies (7.5 < z < 13.5) in the “blank” fields of the Prime Extragalactic Areas for Reionization and Lensing Science (PEARLS) survey combined with Early Release Science data from the CEERS, GLASS, and NGDEEP surveys/fields and the first data release of JADES. We use strict quality cuts on EAZY photometric redshifts to obtain a reliable selection and characterization of high-redshift (z > 6.5) galaxies from a consistently processed set of deep, near-infrared imaging. Within an area of 180 arcmin2, we identify 1046 candidate galaxies at redshifts z > 6.5 and we use this sample to study the UV LF in four redshift bins between 7.5 < z < 13.5. The measured number density of galaxies at z = 8 and z = 9 matches those of past observations undertaken by the Hubble Space Telescope (HST). Our z = 10.5 measurements lie between early James Webb Space Telescope (JWST) results and past HST results, indicating cosmic variance may be the cause of previous high density measurements. However, the number densities of UV-luminous galaxies at z = 12.5 are high compared to predictions from simulations. When examining the star formation rate density of galaxies at this period, our observations are still largely consistent with a constant star formation efficiency, are slightly lower than previous early estimations using JWST, and support galaxy driven reionization at z ≤ 8.


Introduction
One of the main design goals of the James Webb Space Telescope (JWST) is to expand the present redshift frontier and search for the very first galaxies that formed in the Universe.Major science goals are to learn more about the physics of early galaxy and star formation, the nature of dark matter, and to test our understanding of wider cosmology.A number of Guaranteed Time Observations (GTO) and Early Release Science programs were subsequently developed to use this new facility to conduct this search.These include, but are not limited to, a slew of early papers that have investigated potential distant galaxies (e.g., Castellano et al. 2022;Naidu et al. 2022;Finkelstein et al. 2023aFinkelstein et al. , 2023b;;Adams et al. 2023;Atek et al. 2023;Austin et al. 2023;Donnan et al. 2023;Harikane et al. 2023;Leung et al. 2023;McLeod et al. 2024;Willott et al. 2023;Yan et al. 2023b).
This search for the earliest galaxies in the Universe has a long history, dating back to the work of Hubble & Humason (1931) and advancing up until today.Samples of galaxies existing within the first 10% of the history of the Universe now number in the tens of thousands (Duncan et al. 2014;Bouwens et al. 2015;Finkelstein et al. 2015;Adams et al. 2021;Harikane et al. 2022b), many of which were found thanks to synergies between the Hubble Space Telescope (HST) and ground-based facilities such as Subaru and VISTA.Until the launch of JWST, the most distant galaxy candidates found had redshifts up to z ∼ 10 (e.g., Bouwens et al. 2011Bouwens et al. , 2015;;Finkelstein et al. 2015;McLeod et al. 2016;Morishita et al. 2018;Oesch et al. 2018;Salmon et al. 2018;Stefanon et al. 2019;Bowler et al. 2020) and claims of the detection of early galaxies stretch as far as z ∼ 13 (Harikane et al. 2022a).The spectroscopic confirmation of these sources has since followed, with examples at z > 7.5 including Oesch et al. (2015), Zitrin et al. (2015), Roberts-Borsani et al. (2016), andHashimoto et al. (2018), with the highest-redshift source prior to JWST was measured to exist at z = 10.957(GNz-11: Oesch et al. 2016;Jiang et al. 2021).
Over the course of the year of 2022, JWST has been successfully launched, calibrated, and yielded its first data releases, with thanks to our colleagues working for various national and international organizations (Rigby et al. 2023).As a result, a revolution in our understanding of the early Universe is now underway and the first photometric data are now available and being studied in some detail (e.g., Treu et al. 2022;Finkelstein et al. 2023a, C. J. Conselice et al. 2024, in preparation).The first JWST data were released in 2022 July, triggering rapid efforts to identify and follow up candidate galaxies at z > 11 (Castellano et al. 2022;Naidu et al. 2022;Finkelstein et al. 2023a;Adams et al. 2023;Atek et al. 2023;Donnan et al. 2023;Harikane et al. 2023;Yan et al. 2023aYan et al. , 2023b)).However, there are inconsistent samples of galaxies being found between the various studies published to date.There are multiple instances of some studies reporting candidate z > 11 galaxies that either have a low-z solution in another work or are found to be insufficiently detected to be considered robust (Bouwens et al. 2023).These results are likely due to the intricacies of reducing data from a new and complex instrument combined with the varied photometricredshift codes and approaches being undertaken.
The NIRSpec instrument has already begun the job of confirming large numbers of high-redshift galaxies (Bunker et al. 2023b;Curtis-Lake et al. 2023;Tang et al. 2023;Wang et al. 2023).However, it has also shown that some of these candidates are actually at lower redshifts.Recently, work by Arrabal Haro et al. (2023a) has shown a candidate z = 16.4 galaxy is a galaxy at z = 4.9 with strong emission lines.This demonstrates that it is important to be extra careful when dealing with candidate galaxies exhibiting potentially extreme properties.Beyond just redshifts, emission lines and their strengths are key to understanding fundamental galaxy properties.For example, the work of Endsley et al. (2023) has shown that emission lines can influence photometry such that 2 dex differences in stellar mass estimations were possible for one of the originally proposed massive sources presented in Labbé et al. (2023).Further to this, Cosmic Evolution Early Release Science Survey (CEERS) NIRSpec observations analyzed in Kocevski et al. (2023) have shown one such massive z = 7-9 galaxy candidate was really a z ∼ 5 active galactic nucleus (AGN).
Beyond identifying galaxies, one of the first, fundamental measurements that can be conducted with a high-redshift sample of galaxies is the ultraviolet luminosity function (UV LF).This is the measure of the number density of galaxies as a function of UV luminosity, typically measured between 1400 and 1600 Å.The evolution of the UV LF as a function of redshift depends on a number of galaxy properties that can subsequently be used to trace how the galaxy population changes with time.These include star formation rates (e.g., Bower et al. 2012;Paardekooper et al. 2013), AGN emission (e.g., Aird et al. 2015;Giallongo et al. 2019;Matsuoka et al. 2019;Niida et al. 2020;Adams et al. 2021;Harikane et al. 2022b), and dust obscuration (e.g., Bowler et al. 2015;Clay et al. 2015;Tamura et al. 2019;Bowler et al. 2020).UV emission from galaxies is also responsible for a key phase of the Universe's evolution, a process known as reionisation (see, e.g., Robertson et al. 2015;Finkelstein et al. 2019, and references therein).There is still debate with regards to whether AGN, massive galaxies, or the more numerous, fainter galaxies provide greater contributions toward the budget of ionizing photons (Madau & Haardt 2015;Yoshiura et al. 2017;Bosch-Ramon 2018;Hassan et al. 2018;Parsa et al. 2018;Dayal et al. 2020;Naidu et al. 2020).
In the first year, several LFs of distant galaxies have been measured using new JWST data.This includes studies such as Finkelstein et al. (2023a), Bouwens et al. (2023), Castellano et al. (2023), Donnan et al. (2023), Leung et al. (2023), McLeod et al. (2024), Pérez-González et al. (2023), and Finkelstein et al. (2023b).In general, these studies probe the LF in the rest-frame UV up to z ∼ 12, but there are some cases of conflicting results.There are several reasons for this, including (1) the limited survey volumes used to measure the UV LF, due to the few early observations available from JWST and (2) potential systematics in previous measurements, such as contaminant lower-z galaxies.What is needed is a combination of different data sets that samples a large area and in which the data have been reduced and analyzed in a consistent fashion.
In this paper we conduct such a consistent reduction and analysis of the deepest Near Infrared Camera (NIRCam) pointings that have been undertaken in the first few months of operations of JWST.This data set includes observations from the Prime Extragalactic Areas for Reionization and Lensing Science (PEARLS; PIs: R. Windhorst andH. Hammel, PIDs: 1176 and2738;Windhorst et al. 2023), the Cosmic Evolution Early Release Science Survey (CEERS, PID: 1345, PI: S. Finkestein; see also Bagley et al. 2023a), the Grism Lens Amplified Survey from Space (GLASS, PID: 1324, PI: T. Treu;Treu et al. 2022), the Next Generation Deep Extragalactic Exploratory Public (NGDEEP) survey (PID: 2079, PIs: Finkelstein, Papovich, and Pirzkal; Bagley et al. 2023b), the JWST Advanced Deep Extragalactic Survey (JADES, PID: 1180, PI: D. Eisenstein; Eisenstein et al. 2023), and the Early Release Observations (ERO) of SMACS J0723.3-7327(SMACS 0723 hereafter; Pontoppidan et al. 2022).We describe each of these surveys in more detail in Section 2.1.This collection is called the EPOCHS data set version 1 (C.J. Conselice et al. 2024, in preparation), which will be used to study the properties of the first galaxies that formed in the Universe.We use these images and photometry to conduct a photometric search for candidate galaxies with redshifts z > 6.5, using the deep, near-infrared coverage now available to provide strong constraints on the Lyman break λ rest < 1216 Å (Guhathakurta et al. 1990;Steidel & Hamilton 1992;Steidel et al. 1996).This work presents measurements of the UV LF and star formation rate density (SFRD) using this new sample.Future works will target other physical properties of these galaxies, such as their star-forming properties (D.Austin et al. 2024, in preparation) and stellar masses (T.Harvey et al. 2024, in preparation).
The outline of this paper is as follows.In Section 2, we describe the imaging data used along with the reduction processes followed to obtain robust photometry.In Section 3 we present our spectral energy distribution (SED) modeling and selection procedures, assess their reliability, compare to other samples, and show some highlights of unique objects.In Section 4 we use our sample to measure the UV LF in four redshift bins in the range 7.5 < z < 13.5.In Section 5 we discuss these results in the context of the UV LF's evolution and the implications on the SFRD in the early Universe.We summarize our findings and conclusions in Section 6.We assume a standard cosmology with H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7.All magnitudes listed follow the AB magnitude system (Oke 1974;Oke & Gunn 1983).

The Surveys and Fields Used
The data we use for this analysis originate from several sources, including the Early Release Science programs of CEERS, GLASS, and NGDEEP alongside the ERO of SMACS 0723.We also use the current PEARLS Survey GTO fields of El Gordo, MACS J0416.1-2403(MACS 0416 hereafter), and the North Ecliptic Pole Time Domain Field (NEP-TDF; Windhorst et al. 2023).These data primarily include observations taken with NIRCam (Rieke et al. 2005(Rieke et al. , 2008(Rieke et al. , 2015(Rieke et al. , 2023a)).In order to generate a data set which has been consistently handled, we reprocess all of the NIRCam imaging from their lowest-level, raw form obtained from the Mikulski Archive for Space Telescopes (MAST) database.
A more general discussion of these fields and their suitability for finding high-redshift galaxies is presented in C. J. Conselice et al. (2024, in preparation).While the public fields of CEERS, NGDEEP, GLASS, and SMACS have been extensively used in the past, this is not the case for the PEARLS data sets.For the three lensing fields of SMACS 0723, MACS 0416, and El Gordo, one of the two NIRCam modules is positioned to cover the lensing cluster.The second module is located 3′ to the side, with its precise location determined by the position angle at which the observations were taken.While we reduce both modules in these fields, we do not include the cluster regions in this study so as to simplify the later analysis.The inclusion of these clusters would require the use of strong gravitational lensing models and consideration of contamination from intracluster light.A more detailed look into lensed objects located behind the clusters is being coordinated as a parallel effort to this work (Carleton et al. 2023;Diego et al. 2023;Frye et al. 2023;Kamieneski et al. 2023;R. Bhatawdekar et al. 2024, in preparation).

PEARLS
For this work, we include recent GTO observations from the PEARLS survey (PI: R. Windhorst andH. Hammel, PIDs: 1176 and2738) survey (Windhorst et al. 2023).These currently consist of two fields targeting gravitationally lensing galaxy clusters and one blank field.The two clusters we use are MACS 0416 and El Gordo.The blank (noncluster) field is located near NEP-TDF (Jansen & Windhorst 2018).This field consists of four sets of partially overlapping pairs of NIRCam pointings, with each pair of NIRCam pointings referred to as a "spoke."As a GTO survey, PEARLS has proprietary time on its data, but the first set of NIRCam imaging of the NEP-TDF were made publicly available immediately.El Gordo and the first epoch of MACS 0416 are also now public.For each of the two clusters, the cluster proper is centered in one NIRCam module, with the second NIRCam module placed in a random, 3 arcmin distance, parallel location depending on the date and angle of the observation undertaken.MACS 0416 was observed three times at different epochs and position angles, generating a deep, single module image of the cluster and three small parallels.This work utilizes all three epochs of the MACS 0416 observations, meaning there are three NIRCam modules worth of data offset from the primary lensing cluster.All PEARLS observations used eight NIRCam photometric bands, F090W, F115W, F150W, F200W, F277W, F356W, F410M, and F444W.Within NEP-TDF, deep HST/Advanced Camera for Surveys (ACS) imaging in the F606W filter is also present through GO-15278 (PI: R. Jansen) and GO-16252 and GO-16793 (PIs: R. Jansen and N. Grogin) between 2017 October 1 and 2022 October 31.Mosaics of these data, astrometrically aligned to Gaia/DR3 and resampled on 0 03 pixels, were made available prepublication by R. Jansen & R. O'Brien (2024, private communication).

SMACS 0723
Observations of the SMACS 0723 galaxy cluster were undertaken as part of the ERO program (PID: 2736, PI: K. Pontoppidan; Pontoppidan et al. 2022) For this study, we utilize the six-band NIRCam photometry taken in the F090W, F150W, F200W, F277W, F356W, and F444W broadband filters.While we include high-z candidates from these observations in our final catalogs, we do not use this field when measuring the UV LF within this paper.While a number of high-redshift galaxies have been previously identified in this field (Yan et al. 2023b;Adams et al. 2023;Atek et al. 2023), the lack of the F115W photometric band leads to significant photo-z scatter in the redshift range of 7 < z < 10.NIRSpec observations of four, high-redshift sources have shown that photo-zs were systematically overestimated, largely due to the lack of constraining power on the Lyman break positioning (Carnall et al. 2023;Trussler et al. 2023).The use of this field in measuring the UV LF and other population based statistics is subsequently at risk of scatter between redshift bins due to these uncertainties.

GLASS
The GLASS observational program focuses primarily on the Abell 2744 galaxy cluster with a selection of JWST instrumentation (PID: 1324, PI: T. Treu;Treu et al. 2022).But in parallel to these observations, the GLASS program has generated one of the deepest NIRCam imaging sets publicly available.GLASS contains two overlapping parallel NIRCam observations in seven filters (F090W, F115W, F150W, F200W, F277W, F356W, and F444W).This field has already provided several strongly detected high-redshift candidates up to z = 12.5 (Castellano et al. 2022;Naidu et al. 2022;Castellano et al. 2023).This study makes use of both epochs of GLASS NIRCam observations.A wider reduction of the Abell 2744 region (including the UNCOVER program) will be included in future analysis (see Bezanson et al. 2022, for details about UNCOVER).

CEERS
This study also makes use of both primary observing runs (2022 July and 2022 December) of CEERS (PID: 1345, PI: S. Finkelstein; see also Bagley et al. 2023a).This consists of 10 NIRCam pointings mosaiced over the Extended Groth Strip (EGS: Groth et al. 1994) with seven photometric bands (F115W, F150W, F200W, F277W, F356W, F410M, and F444W).This field provides the single largest area used in our study at 64.15 arcmin 2 after masking.The field does lack the F090W band and so we include HST CANDELS imaging in the F606W and F814W filters that was reduced by the CEERS team (their HDR1) following Koekemoer et al. (2011) in order to cover this bluer wavelength range.

NGDEEP
This study makes use of the NGDEEP survey (PID: 2079, PIs: S. Finkelstein, Papovich, and Pirzkal;Bagley et al. 2023b;Leung et al. 2023) and our reduction is updated relative to the initial work conducted in Austin et al. (2023).In particular, flatfield improvements in the NIRCam pipeline have led to great improvements in the red images, though wisp subtraction remains the main limiting factor in the blue images.NGDEEP consists of NIRCam imaging that was run in parallel to NIRISS spectroscopy of the Hubble Ultra Deep Field (HUDF).The NIRCam imaging covers part of the HUDF-Par2 parallel field and consists of six broadband filters, F115W, F150W, F200W, F277W, F356W, and F444W, all with average depths of m AB > 29.5.Due to the lack of data below 1 μm, we add in existing HST data from F606W and F814W through the use of the most recent GOODS-S mosaic (v2.5) from the Hubble Legacy Fields team (Illingworth et al. 2016;Whitaker et al. 2019).

JADES
The first data release (DR1) of JADES (PID: 1180, PI: D. Eisenstein; Hainline et al. 2024b;Bunker et al. 2023a;Eisenstein et al. 2023) was conducted in 2023 June, including full mosaics using the pmap1084 calibrations (Rieke et al. 2023b).The release consists of six overlapping NIRCam pointings from their DEEP tier which uses the filters F090W, F115W, F150W, F200W, F277W, F335M, F356W, F410M, and F444W.The field location is around the HUDF in GOODS-S.In this work, we rereduce this data using our own pipeline for consistency with our other fields.We find that our reductions have average depths that are around 0.1 mag shallower than the official reductions, with one small region of the field affected by residual wisping in the F150W and F200W bands.We employ the same HST/ACS F606W mosaic of the wider GOODS-S region as used in our NGDEEP analysis.

The Image Reduction Pipeline
All of the uncalibrated lower-level JWST data products ("uncal.fits") in this study have been processed following our modified version of the JWST official pipeline.This is a similar process to that used in Ferreira et al. (2022) and Adams et al. (2023), but with minor changes based on improvements developed over the first year of handling JWST data.The full pipeline can be summarized as follows: (1) we use v1.8.2 and CRDS v1084, which contained the most up-to-date NIRCam calibrations at the time of writing and includes the third round of in-flight calibrations.We note that v1084 leads to substantial improvements in the flat-fielding of the red NIRCam bands over the previous major version v0995.(2) We subtract templates of wisps, artefacts present in the F150W and F200W imaging, between stage 1 and stage 2 of the pipeline.(3) After stage 2 of the pipeline, we apply the 1/f noise correction derived by Chris Willott.24 (4) We skip the sky subtraction step from stage 3 of the pipeline.Instead, we perform background subtraction on each NIRCam frame between stage 2 and stage 3 ("cal.fits"files), allowing for quicker assessment of the background subtraction performance and fine-tuning.This consists of an initial constant/flat background subtraction followed by a two-dimensional background subtraction using photutils (Bradley et al. 2022).(5) After stage 3 of the pipeline, we align the final F444W image onto a Gaia-derived world coordinates system (WCS; Gaia Collaboration et al. 2018, 2023) using tweakreg, part of the DrizzlePac python package, 25 and then match all remaining filters to this derived WCS, ensuring the individual images are aligned to one another.We then pixel match the images to the F444W image with the use of reproject (Hoffmann et al. 2021). 26he final resolution of the drizzled images is 0 03 pixel −1 .We note that this reduction differs from the official PEARLS and CEERS reductions in Windhorst et al. (2023) and Bagley et al. (2023a).We compare the photometry extracted from the different reduction pipelines in Section 2.4.

Source Photometry and Cataloging
With the final mosaics of the field completed, we carry out source identification and extraction.For this process, we use the code SExtractor (Bertin & Arnouts 1996).We run in dual-image mode with a weighted stack of the red wide-band images (F277W, F356W, and F444W) used for object selection, from which we carry out forced aperture photometry for multiband measurements.This photometry is calculated within 0 32 diameter circular apertures and we include an aperture correction derived from simulated WebbPSF pointspread functions (PSFs) for each band used (Perrin et al. 2012(Perrin et al. , 2014)).This diameter was chosen to enclose the central/ brightest 70%-80% of the flux of a point source, but small enough to avoid contamination.This enables us to balance the use of the highest-signal pixels when calculating fluxes, while avoiding the dependence on a correction derived from the PSF model that is as high (or higher) than the measurement made.
We calculate the depth of our final images by placing circular apertures in "empty" regions of the image.An "empty" aperture is one where no preexisting source in the image, as identified in our source extraction, is located within 1″ of the central coordinate of our aperture.These apertures are used to derive an average depth for each field as well as calculate local depths across each field.The final photometric errors for each source are calculated by using the normalized median absolute deviation (NMAD; Hoaglin et al. 1983) of the nearest 200 empty apertures to calculate the local depth.This process is required in order to generate realistic photometric errors, which SExtractor is known to underestimate.The average field depths in each photometric band are presented in Table 1.
Each of our fields is masked by eye.These masks cover diffraction spikes, the few remaining snowballs, regions of high-intensity intracluster medium (in the NIRCam modules containing any foreground cluster), and a buffer around the edges of the images.The edges of the images are typically shallower due to the dithering patterns used by JWST observations.The total amount of unmasked area used in this study is listed alongside the average depths of each field in Table 1.

Comparing Image Processing Pipelines
At the time of writing, the PEARLS team has developed an official internal pipeline, dubbed ProCess, which has been used to generate official NIRCam reductions for each of its fields using the ProFound software package (Robotham et al. 2018;A. Robotham et al. 2024, in preparation).The key difference is ProFound conducts its own implementation of background subtraction and wisp removal (see Appendix A of Windhorst et al. 2023;Robotham et al. 2023, for more details and comparisons).Similarly, the CEERS collaboration has also generated a slightly different pipeline and publicly released reductions of four NIRCam pointings in their v0.5 and v0.6 data releases (Bagley et al. 2023a).To assess if these reduction procedures result in systematic differences when measuring the final photometry, we comparison our photometry with that extracted using SExtractor with the same settings on each of the images made by the PEARLS and CEERS teams.It should be noted that the CEERS images are produced using pmap0995, and so we compare pmap0995 version of our CEERS images to theirs in order to assess differences in the pipelines properly; no other changes other than CRDS version are used.
To conduct our comparisons, we first generate a catalog of galaxies selected in each band and crossmatch to the one derived from the F444W image.We then measure the separation of the matches to check if each of our images are correctly aligned with each other.It is important that the alignment between these images is very good for extracting correct colors with forced photometry.We find that the typical source centroid offset is less than 0 02 across the fields analyzed, or less than two-thirds of a pixel.
After this, we compare both aperture and total (AUTO) photometry extracted from our final F444W selected catalogs to those extracted from the PEARLS team's reductions.We calculate the mean magnitude offset of sources between apparent magnitudes of 22 < m < 26 in each photometric band, ensuring we are using bright, but unsaturated, sources.We find that our photometry is accurate relative to other reductions, with mean offsets of less than 0.03 mag in the blue photometric bands and 0.01 mag in the red photometric bands.Similarly, the astrometry is consistent to less than 0 07.Our reductions currently only use Gaia stars in the field to calibrate the WCS, while the PEARLS reductions are matched to Gaia-aligned catalogs generated from wider-field, ground-based imaging.The PEARLS reductions thus use a larger set of secondary targets calibrated to the Gaia system when constructing the JWST astrometry and are expected to have higher performance compared to our own.However, this finding shows that this offset is at the <2 pixel level.The comparison to the official CEERS collaborations reductions reaches the same conclusion.Future astrometric improvements, e.g., using complementary ground-based data, will be implemented in the second generation of the EPOCHS sample in order to aid spectroscopic follow up.Note.Depths are calculated by placing random apertures in regions of the image that are empty in the final segmentation maps used to make our catalogs.We perform this process on each of the fields used in this study.The NEP-TDF, El Gordo, and MACS 0416 fields are components of the wider PEARLS survey, while SMACS 0723, GLASS, CEERS, JADES DR1, and NGDEEP are public observations.Source catalogs use local depths calculated with the NMAD of the nearest 200 empty apertures, deriving final photometric errors on a source-by-source basis.The four NEP-TDF spokes and 10 CEERS pointings have mean depths consistent with one another to within 0.1 mag, with the exception of CEERS-P9 which has additional exposures in F115W and F444W that increase the depths in these two bands by more than 0.2 mag.Where applicable, HST depths are also shown.For the two fields within GOODS-S (NGDEEP and JADES), the HST data have two tiers of depth, and we show the average depth for each tier.The bottom row shows the total area used in square arcminutes, after the masking was applied and cluster modules omitted.

A Robust Sample of Ultrahigh-redshift Sources
With the imaging and cataloging completed, we fit our NIRCam-derived SEDs for all sources identified in order to derive photometric redshifts for each source.We conduct this process with the use of EAZY-py, the python implementation of the EAZY photometric-redshift code (Brammer et al. 2008).We include the newly developed templates presented in Larson et al. (2023b).These templates expand upon one of the default template sets, which use 12 templates generated using the Flexible Stellar Population Synthesis code (Conroy & Gunn 2010) to include galaxies which are bluer and have stronger emission lines, which we expect are more appropriate for z > 8 galaxy SEDs.We also test the SED templates used in the JADES works (Hainline et al. 2024b), but find that the Larson et al. (2023b) templates are slightly better performing when compared to current spectroscopic results.We also tested the use of a second rapid SED fitting code LePhare (Arnouts et al. 1999;Ilbert et al. 2006) with the use of Bruzual & Charlot (2003) templates, but find that emission line strengths often struggled to fit medium-band excesses which could lead to the rejection of legitimate sources when quality cuts are applied.Some photo-z codes have mechanisms for fine-tuning the zero-points of the photometric bands in order to optimize the results against a provided spectroscopic sample.We do not employ these techniques because each of the individual chips that make up the NIRCam modules (eight in the blue and two in the red) have their own independent calibrations and photometric zero-point offsets.We would thus require each of these zero-point modifications to be made on a chip-by-chip basis as opposed to the final mosaic.The number of objects with spectroscopic redshifts within the area of a single chip is inherently small due to the small field of view that each covers.As a result, the zero-point modifications we could apply may be biased toward certain galaxy colors, depending upon the type of spectroscopically confirmed galaxies within each module.Following discussion with members of the community, residual zero-point errors were anticipated to be approximately 5%.We subsequently implemented a minimum 10% error on the measured photometry in order to account for potential zeropoint issues within the NIRCam reduction pipeline in addition to other systematics such as minor reduction imperfections and errors in the template SED shapes themselves.
To select a robust sample of high-redshift galaxies, we employ the following criteria: 1. 3σ detection in all bands blueward of the expected Lyman break.2. 5σ detection in the first two bands directly redward of the Lyman break, and 2σ detection in all other redward bands, excluding F410M.If the galaxy appears only in the long-wavelength NIRCam photometry (i.e., an F200W or higher dropout), we increase the requirement to a 7σ detection to minimize the potential for selecting red camera artefacts, which can appear to look like noisy 16 < z < 18 sources.for the best-fitting SED to be classed as robust (good).
5. Δχ 2 4 between high-z and low-z EAZY runs (where the maximum redshift is set to 6).This ensures that the high-z solution is much more statistically probable.6.If the 50% encircled flux radius (FLUX_RADIUS parameter in SExtractor) is smaller than the FWHM of the PSF in the F444W band, then we require that Δχ 2 4 between the best-fitting high-z galaxy solution and the best-fitting brown dwarf template (see below for further discussion).7. 50% encircled flux radius is 1.5 pixels in the longwavelength, wide-band NIRCam photometry (F277W, F356W, and F444W).This avoids detecting hot pixels in the long-wavelength detectors as F200W dropouts.
Targets must pass all of the above criteria to fall within our parent sample of high-z galaxies.We limit our parent high-z sample to having a redshift above z > 6.5.

Identifying and Removing Brown Dwarfs
Despite the wavelength coverage provided by JWST, local brown dwarf stars can be misidentified as high-redshift galaxies and enter high-z samples if not properly accounted for.In order to address this, we conduct SED fitting of brown dwarf templates from the Sonora Bobcat model set (Marley et al. 2021;Hainline et al. 2024a).We then remove galaxies identified as being both better fit by a brown dwarf template and have a size consistent with the FWHM of the PSF in the F444W band.This process removed 56 objects from our sample.The majority of these (39 objects) had galaxy photo-zs of 6.5 < z < 7.5.This procedure is only run on objects initially identified as being robust high-z candidates and not on every single source, so it does not amount to a complete brown dwarf sample though this may be explored in future work.As a sanity check, we fit the objects identified as brown dwarfs in the CEERS field by Hainline et al. (2024a), which all did not pass all of our quality checks to enter our high-z sample in the first place, and find that we would have still identified these objects as potential brown dwarfs.

Testing the Selection Procedure
To validate the photometric redshifts, we compare our photo-z measurements to catalogs of spectroscopic redshifts available in these fields.Within SMACS 0723, we compare to the 10 objects covered by NIRSpec in Carnall et al. (2023).In El Gordo and MACS 0416 we compare to the MUSE observations reported by Caminha et al. (2017Caminha et al. ( , 2023)).Within CEERS, we combine catalogs obtained from the DEEP2 (Newman et al. 2013), 3D-HST (Momcheva et al. 2016), and MOSDEF/MOSFIRE (Kriek et al. 2015) programs to construct a catalog of 1500+ spectroscopic redshifts.A summary of our findings is presented in Figure 1.
For the SMACS 0723 cluster, we find that 6/10 sources from Carnall et al. (2023) have photometric-redshift solutions in our catalogs within 15% of the spectroscopic solution.Three of these outliers are high-redshift sources, with photometric redshifts z phot > 9 while the spectroscopic redshift is around z ∼ 7.6.While correctly identified as being high redshift, the lack of F115W leads to large uncertainties on the Lyman break positioning.Without this constraint, photo-z codes favor fitting the observed F444W excess of these sources with a Balmer break at z ∼ 9.5, as opposed to [O III]/Hβ line emission at z ∼ 7.6 (see Trussler et al. 2023, for a detailed discussion).
While the data available in this field have proven their capabilities in correctly identifying high-z from low-z targets, these broad uncertainties in the redshift range of 7 < z < 10 make their use in the measurement of the UV LF inherently risky due to the low number statistics combined with the potential for sources to scatter between redshift bins.
Within El Gordo, we test our photo-z procedure by matching our photometric sample to 402 spectroscopic redshifts obtained by the MUSE instrument (Caminha et al. 2023).The sample contains a large number of cluster members at z ∼ 0.87, as well as a wide range of galaxy redshifts from 0 < z < 6.We obtain crossmatches for 238 sources with confidence flags indicating a secure redshift.Our EAZY photometric procedure obtains 218/238 (92%) redshifts.Similar to El Gordo, we match our catalogs of MACS 0416 to MUSE spectroscopic catalogs (Caminha et al. 2017).We find 261 crossmatches, of which 235/261 (90%) have photometric redshifts within 15% of the spectroscopic redshift.
From the CEERS survey itself, we compare to 19 high-z spec-z galaxies (Arrabal Haro et al. 2023b;Tang et al. 2023) and find a 79% spec-z recovery rate.Of the failures, two are fainter than our 5σ limits, one is a partially blended source, and one is still classed as high-z but our redshift was too high (z = 10.9 instead of the z = 8.88 spec-z).We successfully recover a number of significant high-z sources, such as the the AGN candidate at z = 8.68 identified in Larson et al. (2023a) with z phot = 8.65 and "Muse's galaxy" at z = 11.44 with z phot = 11.90 (Finkelstein et al. 2023a) In total, 11/15 identified high-z galaxies enter our final sample, with two galaxies in our masking near image edges and two more below our detection limits.The CEERS team implements smaller extraction apertures (0 2 diameter) than in this study (0 32), enabling fainter sources to be identified for spectroscopic follow up, though with increased risk of influence of bad pixels and reliance on correct modeling of the PSF.Consequently, it is a common occurrence for CEERS selections to identify sources which are approximately 4σ in our work.
Finally, we test the setup of our photometric-redshift codes on the early photometry and spectroscopy released from the JADES collaboration (Curtis-Lake et al. 2023; Robertson et al. 2023).Using the DR1 PRISM catalog, there are 144 spectroscopic crossmatches with our catalogs.Of these, 97/124 (78%) unmasked objects have photometric redshifts within 15% and all 15 objects at z > 6.5 have the correct redshifts.12/ 15 meet our robust criteria, with those rejected failing our detection thresholds.Further investigation identifies that two of these three objects lie in a small region of poor wisp subtraction, leading to shallow local depths.

Source Detection
We carry out completeness tests on our data using a simple simulation to extract artificial sources that are inserted into the images.We insert galaxies with a Sérsic index of n = 1, and absolute luminosities ranging from −16 to −24 into the F444W image used for our selection.Half-light radii of the artificial sources are varied and follow the distribution of theoretical predicted sizes of distant galaxies from the BlueTides simulation (Feng et al. 2016;Wilkins et al. 2017;Marshall et al. 2022).The maximum half-light radius of our simulated galaxies was placed at 1000 pc.With the exquisite resolution provided by JWST, galaxies at ultrahigh redshifts are often slightly larger than the PSF, which can be a key factor that must be considered when trying to understand completeness.
The simulations we carry are are performed by sweeping the magnitude range in 0.2 mag intervals and with redshift intervals of 0.5, placing 1000 artificial galaxies in the field of interest for each interval, and then running SExtractor with the same settings as for our unmodified images.A size-luminosity relation of r ∝ L 0.5 was assumed throughout (Grazian et al. 2012), with a reference radius of 800 pc at M UV = −21.This allows us to determine how efficiently our setup detects galaxies at various redshifts, and what fraction enter our catalogs.This is necessary for using our galaxy counts to produce a LF.This procedure is carried out for each individual field.

Selection Procedure
To test our photometric redshifts and selection procedure further, we employ the use of simulated photometric catalogs generated by the JAGUAR semianalytic model (Williams et al. 2018).This model is selected due to its performance in replicating colors of known high-z galaxies (Wilkins et al. 2023).In order to have reasonable numbers of luminous The CEERS data include HST F606W and F814W imaging due to its lack coverage below 1 μm.Despite being largely limited to NIR wavelength range, photo-zs generally perform well, with the greatest spread around 2 < z < 3. The blue line shows the ideal case of the photometric redshifts matching the spectroscopic redshifts, while the red lines show 15% offsets in 1 + z, often used to describe the performance of photometric redshifts.
high-z objects, we use five realizations of the simulation, producing around 550 arcmin 2 of simulated area (10 times that of any one field used in this work).From these simulations, we add photometric errors and scatter the true photometry in accordance to the average depths for each of our survey fields.We then run our photometric-redshift and selection pipeline on the resultant catalogs to assess the completeness and contamination (which we define as sources at z < 5) that fall into luminosity bins (Δm = 0.25) in the photometric band containing the rest-frame UV per redshift interval of Δz = 1.
We find that our selection criteria are the dominant limitation in completeness as opposed to source detection.This is largely because we utilize a stack of the deepest bands in the reddest wavelengths for detection but employ 5σ cuts in the blue bands in our selection process.As is anticipated, the completeness is high (85%-95%) for the most luminous sources and rapidly plummets as the 5σ criterion is reached.Of interest is the contamination of low-z sources and how this behaves with apparent magnitude.For bright sources, this is unsurprisingly very low (0%-5%).However, we find that contamination begins to rise just before completeness drops off, with upwards of ∼10% levels of contamination in the 0.25 mag wide magnitude bin just before the 5σ limits are reached.We show an example of this simulations applied to the JADES survey design in Figure 2.
For the purpose of measuring the UV LF in this work, we only use the parameter space identified as being greater than 50% complete in our simulations.This is to avoid depending on using corrections that are as large, or larger, than any measurements being undertaken.We additionally do not consider the use of a field in the UV LF when the contamination level is greater than 10%.This conservative limit is selected (and a cut implemented as opposed to correcting for it) because the contaminating sources themselves are poorly understood and their relative distribution within the JAGUAR simulation may not replicate reality.We therefore cut at the point we start to see an uptick in contamination.Ultimately, obtaining spectroscopy of contaminating sources, in addition to real high-z sources, will be key to understanding what sources really are falling within the selections and how numerous they may be.

The Final Sample
Applying this procedure to our catalogs for each field, we obtain a final sample of 1046 high-redshift candidates (z > 6.5) over an area of 181 arcmin 2 .Of these, 413 have best-fitting redshifts between 7.5 < z < 13.5.This is one of the largest samples (in area and raw numbers) of high-z, JWST-selected candidates to date.The full list of candidates will be presented in the companion paper (C.J. Conselice et al. 2024, in preparation).To provide some context for this work, we present a selection of SED fits for some "highlight" sources, as well as some random sources, in Appendix B.

The Most UV-luminous Sources
The parent galaxy sample contains over 1000 sources located within the first billion years of the history of the Universe.Here, we highlight a few examples of extreme sources used in the measurement of the UV LF at z > 7.5.
The known AGN candidate at z = 8.68 identified in Larson et al. (2023a) remains the most intrinsically luminous source in our compilation of JWST fields with M UV = −22.28.The next brightest source at z > 7.5 is identified as NEP-3:15371 in our catalogs with M UV = −21.34and z = 9.00 ± 0.1.NEP-TDF contains a relatively large number of UV-luminous sources at z ∼ 8.1-8.3 including five sources with M UV < −21.0, which may be indicative of an overdensity.At z > 10, the brightest source is NEP-3:18221 with z = 11.5 ± 0.4 and M UV =−21.03, which is followed by the two sources located within the GLASS field at z = 10.7 and z = 12.4 identified in Naidu et al. (2022) and Castellano et al. (2022).This indicates how fortuitous the GLASS galaxy candidates were, given this study probes nearly 20 times the area and identifies only a single intrinsically brighter source at z > 10.SED examples for these luminous galaxies are presented in Appendix B.  which sources appear in both our catalogs and those generated by these works, where the matching radius used is 0 3.A summary of our findings is presented in Table 2.We find that on average, only about 20%-40% of sources appear across multiple early studies, a conclusion also reached by Bouwens et al. (2023) for their own sample.For the sources that are selected by other works and not in our own, we find that about half of them do have a primary high-redshift solution and half of them are fainter than our detection limits.Since the shape of the UV LF means that there are more galaxies at fainter luminosities, the scattering of sources above and below these detection limits can be significant.This is especially so when considering different studies that use different apertures for measuring photometry (resulting in different limiting magnitudes; see, e.g., Finkelstein et al. 2023a, who use 0 1 radius apertures).We note one main exception for our comparisons with early studies is Bouwens et al. (2023), who select 10 luminous objects at z ∼ 6.5 which have a Lyman break located within the F090W JWST band and hence fail our criterion of having a fully clean band beyond the break.The numbers in Table 2 are consequently not as first appears for this reason.More detailed comparisons are presented in Appendix A.

Comparisons to Samples from Other Studies
Perhaps reassuringly, more recent studies using newer calibrations and implementing improved knowledge of high-z selection and contamination management have greater degrees of overlap.Our strongest sample overlap is with the work by McLeod et al. (2024), with 75% of their sources entering our final sample and with 88% of sources agreeing on primary photo-z solutions.Also recently, the CEERS team published their official list of highz candidates (Finkelstein et al. 2023b).While only 33% enter our final sample, the primary cause of this is the aperture size choice.This means there are a large fraction of objects which lie just below our 5σ limits and this is the primary reason they do not make our sample, although we tend to agree (70%) that these sources have strong high-z solutions.

Background
Using our final list of candidate galaxies with fitted SEDs and an understanding of our completeness, we construct an initial UV LF at z > 8 from these JWST data.To conduct this measurement, we use the 1/V max method (Rowan-Robinson 1968; Schmidt 1968).The maximum observable redshift for each source z det is calculated by iteratively shifting the bestfit SED of each source in small steps of δz = 0.01 and convolving with the filter sets until it no longer meets our restframe UV detection criteria.The object can thus be detected within a comoving volume (V max ), which is the volume contained within a section of a spherical shell between z z z min max < < .Here, z min is the lower boundary of the redshift bin being considered and z max is either the maximum bound of the redshift bin or z det if it is found to be lower.We also experiment with the Lynden-Bell C − method (Lynden-Bell 1971; Woodroofe 1985; Wang et al. 1986;Efron & Petrosian 1992) and find results that are generally consistent to within 0.5σ.
To calculate the M UV value for each galaxy, we convolve a 100 Å wide top-hat filter centered on the rest-frame 1500 Å for the best-fitting SED (Finkelstein et al. 2015;Adams et al. 2021;Finkelstein et al. 2023a).These SEDs are generated using aperture-corrected photometry, which can underestimate the total flux for resolved sources.We therefore apply a simple correction factor using the ratio between the aperture and total photometry (MAG_AUTO) measured by SExtractor in the photometric band containing the rest-frame 1500 Å.
For the measurement of the UV LF, we discard the SMACS 0723 field.With the absence of the F115W photometric band, there is a much greater degree of uncertainty on candidate redshifts between 8 < z < 10 (see, e.g., Trussler et al. 2023).Since we are not using the NIRCam modules containing lensing clusters in this study, the SMACS 0723 field provides minimal cosmic volume, but comes with an elevated risk of objects scattering between redshift bins.

The UV LF Construction
Implementing the above, the rest-frame UV LF (Φ(M)) is calculated using: where ΔM is the width of the magnitude bins and C i,f is the completeness correction for galaxy i depending on its location within field f.We only use fields and redshift ranges in this calculation if their completeness is greater than 50% and the contamination is less than 10%.These ensure that our measurements are not dominated by any derived corrections (similar to the aperture choice when doing PSF corrections).
The magnitude bin widths used are 0.5 mag, increasing to 1.0 mag for the faintest bins.
The measured uncertainty of the LF is given by: For small numbers of objects (N < 5), the above formula overestimates the Poisson error (leading to instances of, e.g., 1 ± 1).In this regime, we directly calculate the Poisson confidence interval using an improved, but still conservative, estimator based on the χ 2 distribution 1990).We measure the UV LF in four primary redshift bins: 7.5 < z < 8.5, 8.5 < z < 9.5, 9.5 < z < 11.5, and 11.5 < z < 13.5.Taking into account the completeness criteria of each field, these four redshift bins have a total of 166, 59, 27, and 9 galaxies respectively.We note that our full sample of galaxies selected is larger than that used to measure the UV LF.This is due to the nature of the shape of the UV LF, where there is a greater number density of fainter sources occupying the <50% completeness regime.The final results are presented in Table 3.

Cosmic Variance
Pushing the redshift frontier, the implications of cosmic variance on the measurement of the UV LF are unclear (Dawoodbhoy et al. 2023).To obtain a ballpark figure for the systematic errors introduced by cosmic variance, we utilize the cosmic variance calculations introduced in Driver & Robotham (2010).We calculate the variance from each individual field using its approximate dimensions (e.g., we treat all 10 pointings of CEERS as a single large field) and combine them together using a volume-weighted sum in quadrature (Moster et al. 2011).Across our UV LFs, the impact of cosmic variance is estimated to be approximately 17%-21% at the bright end, and up to 42% at the faint end (where there is only the NGDEEP and JADES fields).Cosmic variance is thus a significant contributor to the overall uncertainty on the measured UV LFs, even when combining several, well-separated JWST fields.We add this error in quadrature with the Poisson error and these are included in the errors reported in Table 3.

Fitting the Measured Ultraviolet Luminosity Function
To fit our UV LFs, we use both the Schechter (Schechter 1976;Sch) and double power-law (DPL) functional forms.We use the Levenberg-Marquardt (LM) approach to minimize the χ 2 and obtain the best-fitting parameters.We then compare this result to one obtained using a full Markov Chain Monte Carlo (MCMC) method using emcee (Foreman-Mackey et al. 2013) with wide and uniform priors.The results of these fits are presented in Table 4.Our observations do not probe enough area to provide constraints on the UV LF much brighter than the "knee" at M UV < −21.5.To assist in our fitting, we include the use of the ground-based, UltraVISTA-derived data points presented in Donnan et al. (2023).The MCMC methodology is found to produce larger, but likely more realistic, errors than the simpler LM approach.Overall, the addition of these extra data points does not change the fitted results, but does give us better accuracy in their measurements.
Due to the few data points present at z = 10.5 and z = 12.5, we find that leaving all variables free leads to unconstrained models (see, e.g., the DPL fit for z = 10.5 in Table 4).We thus fix the faint-end slope for z = 10.5 and all but the normalization for z = 12.5 for both the Sch and DPL functional forms.For the faint-end slope we fix to α = −2.1 and the other parameters are fixed to the best solution from the previous redshift bin.While we provide the various fit values in Table 4, greater depths and volumes are required in order to constrain properly the UV LF at z 12 and generate realistic model fits.

The Measured z = 8 and z = 9 Ultraviolet Luminosity Function
Since we are using only blank-field NIRCam data in this study, the 7.5 < z < 8.5 redshift bin does not include any benefits which can be provided by gravitational lensing.As a result, the JWST observations in this redshift range are more comparable to the observations previously undertaken by the HST and presently do not extend our knowledge of the faint end (see e.g., Livermore et al. 2017;Ishigaki et al. 2018;Bhatawdekar et al. 2019, for lensing assisted studies using the HST).This measurement thus serves as a sanity check against more established LF measurements that lie comfortably within the capabilities of both JWST and HST.
We present our JWST UV LF measurements at z = 8, which use 166 sources, in Figure 3 and find they are systematically high relative to past HST results, particularly at the bright end.Investigations reveal that NEP-TDF contributes a disproportionate amount to these bins.Of interest, is that the brightest z ∼ 8 candidates are largely located within the third spoke of NEP-TDF observations and may be indicative of a potential Table 3 The Measured Rest-frame Ultraviolet Luminosity Function and its Error Margin in the Redshift Bins 7.5 < z < 8.5, 8.5 < z < 9.5, 9.5 < z < 11.5, and 11.5 < z < 13.5 Column (2) shows the number density of objects and column (3) shows the uncertainties in the number density, which are calculated with Equation (2) and summed in quadrature with the derived cosmic variance.The "NoNEP" list of values refers to the 7.5 < z < 8.5 results when the suspected overdense NEP field is omitted.
overdensity.A more detailed clustering analysis of the field is left to future work.Given the NEP-TDF may have an overdensity, we have produced two sets of measurements of the z = 8 UV LF, one with and one without this field.
When the field is omitted, the excess is largely removed and the UV LF well replicates the previous HST observations from    et al. (2021).This validates both the JWST observations and methodologies employed within this study and indicates that there is no significant tension between JWST-selected samples of z = 8 galaxies and those derived using HST data.The dynamical range of luminosity presently probed by this data set is insufficient to constrain the shape of the UV LF strongly.We therefore utilize the faint-end data points (M UV > −19) from McLure et al. ( 2013) to enable higher-quality fits of the UV LF at this redshift.These extra data do not significantly change the fits we obtain, but reduce their uncertainty.The inclusion of gravitational lenses is likely required for JWST-based UV LFs to push fainter at z ∼ 8.The z = 9 UV LF is presented in Figure 4 and is measured using 59 sources.Here, the combined depth and area covered by these JWST programs enable our observations to overlap with the results from the ground-based VISTA programs, which cover square degrees of the sky (Donnan et al. 2023).We find that our measurements are in close agreement with other JWST-based studies and are slightly higher than recent HST observations by Bouwens et al. (2021).The bright end of the UV LF is elevated due to the inclusion of the AGN candidate from Larson et al. (2023a) and this brightest data point is discounted when it comes to model fitting.Compared to Harikane et al. (2023), the total area in this study is a factor of 2 higher and the majority of our area is deeper by upwards of half a magnitude in the blue NIRCam bands, which contain the rest-frame UV emission at these redshifts.Our effective volumes at lower luminosities are thus significantly greater, providing a larger number of sources (59 at z = 9 versus 13) and minimizing cosmic variance.We provide measurements that not only agree with the Harikane et al. (2023) results, but build upon them by reducing the uncertainties.Toward the faint end, we find a slightly lower number density compared to the observations of Donnan et al. (2023).This is possibly due to the SMACS 0723 photo-z uncertainties and the tendency for SED fitting codes to use Balmer breaks to explain any F444W excess.Our fits match the ultradeep MIRI GTO Our z = 10 UV LF, presented in Figure 5, uses a total of 27 sources in its measurement.This redshift range has sparked some initial discussion within the observing community regarding the measured number density of sources and the implications for the UV luminosity density.The work of Oesch et al. (2018), using HST derived results, observed a significant evolution of the UV LF compared to that measured at z = 8 and z = 9.However, Donnan et al. (2023) andMcLeod et al. (2024), using JWST, found a much greater number density of galaxies over the same luminosity range.It is presently unclear if the difference between the JWST-based results and HSTbased results are the consequence of selection procedures or cosmic variance.Many early JWST UV LF measurements have the CEERS survey contribute a large proportion of their volume (see the discussion in Willott et al. 2023) and if this was overdense, it could explain these observed differences.In this work, the CEERS field provides about 36% of our cosmic volume and our number densities are found to be slightly smaller than in other JWST-based works.The inclusion of even more area to these depths (e.g., the full JADES survey and the PRIMER survey) and the inclusion of gravitational lenses will ultimately be required in order to complete our picture of the z = 10.5 galaxy population.
While our sample uses only nine galaxies to measure the 11.5 < z < 13.5 UV LF (see Figure 6), its constraints can already provide us with an insight into the evolution of the very first galaxies.At z ∼ 12.5, we find that our measured number density agrees well with the observations of Donnan et al. (2023), Harikane et al. (2023), McLeod et al. (2024), and Finkelstein et al. (2023b), which sample similar redshift bins.Of note is the high number density of luminous UV sources, which implies a mild evolution in the number density of such sources.This extends the findings of Bowler et al. (2020), who find the bright end of the UV LF evolves mildly to higher redshifts (8 < z < 10) compared to galaxies at the faint end.Bouwens et al. (2022) find a number density of sources that is approximately an order of magnitude greater than the three other studies used for comparison.As noted within that work, the volume used is small, but the use of medium-band filters enables tighter photometric-redshift constraints.It is possible the disagreement is the consequence of large-scale structure, since the medium-band survey covered only 4.6 arcmin 2 .These number densities are high relative to the expectation from hydrodynamical simulations, yet a selection of these sources have already been spectroscopically confirmed (e.g., Maisie's galaxy and those from JADES).More spectroscopy and more area will be required in order to prove the UV LF evolves slower than anticipated at z > 12.

Pushing toward the First Galaxies
To date, there has been a small selection of z ∼ 16 JWST candidates presented (Yan et al. 2023b;Austin et al. 2023;Donnan et al. 2023;Harikane et al. 2023).An initial front runner for an extreme redshift galaxy was a z = 16.4 candidate in the CEERS survey.The first spectroscopy of this z ∼ 16 candidate was recently attempted as a part of a DDT proposal (PID:2750, PI: P. Arrabal Haro).Initial analysis shows that the z = 16.4 candidate is in fact a dusty, extreme line emitter at z ∼ 4.9 (Arrabal Haro et al. 2023a), as predicted independently by Zavala et al. (2023) and Naidu et al. (2022).In this work, this z = 16.4 candidate is initially identified as being very high-z, but later rejected as a z = 4.7 solution, close to the now known specz, and has a χ 2 comparable to the initial high-z solution.
Our work identifies a total of 22 candidate z > 13.5 galaxies meeting our criteria using the template set from Larson et al. (2023b).Of these 22 galaxies, only six also meet the criteria when employing the JADES templates (Hainline et al. 2024b).These are z = 13.9 and z = 15.8 galaxies in NGDEEP (see also Austin et al. 2023;Leung et al. 2023), a z = 14.1 candidate in NEP (though a by-eye examination reveals a suspiciously low F410M measurement), z = 14.6 and z = 18.1 candidates in CEERS pointing 9, and finally a z = 16.1 galaxy in JADES DR1.We do not produce a UV LF measurement using these extreme candidates as their reliability is presently unknown.Either spectroscopy or additional medium-band coverage will be required in order to increase the confidence that these systems are real.Spectroscopic studies between 10 < z < 13.5 have shown photometric selections of galaxies are effective overall and so we limit our UV LF measurements to this regime (Arrabal Haro et al. 2023a;Bunker et al. 2023b;Curtis-Lake et al. 2023;Wang et al. 2023).

Comparing the Measured Ultraviolet Luminosity Function and its Evolution to Theory
We present the time evolution of our DPL fits to the UV LF in Figure 7.To provide better context for the UV LFs presented in this study, and the time evolution observed, we compare our  observations and fits to predictions from a selection of recent simulations.These include Bluetides (Feng et al. 2016;Wilkins et al. 2017), Delphi (Dayal et al. 2014(Dayal et al. , 2022)), DREAM (Drakos et al. 2022), Thesan (Kannan et al. 2022), the Santa Cruz semianalytic model (Yung et al. 2024), SHARK V2.0 (Lagos et al. 2023), and FLARES (Lovell et al. 2021;Vijayan et al. 2021;Wilkins et al. 2023).The distribution of UV LFs generated by these simulations agree more with each other and with the observations at lower redshifts (z = 8), but show an increasing spread toward higher redshifts.Primary trends are that the Delphi and FLARES simulations tend to predict greater numbers of UV-faint galaxies (−18 < M UV < −17) with Thesan also jumping up in number density at M UV ∼ −16, while the Santa Cruz semianalytic model predicts fewer of these galaxies relative to the other models.Simulations begin to diverge in their predictions at z 10.However, the observational errors also increase in this regime, matching the spread of simulation predictions.This means that we are presently not yet able to favor one physical model over another confidently.
At the most extreme redshifts probed (z ∼ 12.5), we find that the predicted number density of UV-luminous sources (M UV ∼ −20.5) varies by up to 2 dex between simulations.It is immediately apparent that much greater volumes at such redshifts are required in order to constrain fully the model predictions in this early epoch.However, current measurements indicate that the number density of UV-luminous sources are at the upper end of the prediction ranges, which agrees most strongly with the FLARES simulation.

The Star Formation Rate Density in the Early Universe
Using the best-fit DPL models to our UV LF measurements, we integrate the LFs as faint as M UV −17 to obtain the luminosity density of UV light.This measurable has previously been used as a proxy for the SFRD.To do so, we replicate the procedure of Harikane et al. (2023), which uses the relations derived in Madau & Dickinson (2014).Here we use the simple conversion factor of k = 1.15 × 10 −28 M e yr −1 (erg s −1 Hz −1 ) -1 .The results of this integral and conversion are presented in Figure 8 alongside recent estimations from the literature.We find that the SFRD estimation at lower redshifts (z = 8-9) agree well between past studies.Of significant interest is our results at higher redshift.Our data point at z = 10.5 is lower than the results of many other JWST-based studies and is well matched to the extrapolation of theoretical models, such as the constant star-forming efficiency model from Harikane et al. (2022a).This is in agreement with the recent results from the CANUCS survey (Willott et al. 2023).At the highest redshifts (z = 12), our results show a high UV luminosity density.However, this excess relative to simple models is only significant to 1σ and is consistent with the findings from Willott et al. (2023).Consequently, significant deviations from simple models (e.g., constant star-forming efficiency) may not necessarily be required in order to explain these excesses.Ultimately, larger samples of galaxies at z > 12 and more spectroscopy will be required in order to limit the impact of cosmic variance and increase confidence in the selection procedures.
If the evolution of the UV luminosity density can be described with a simple linear evolution (as used in Donnan et al. 2023)   The UV luminosity density also provides insight into the process of reionisation.We use the rest-frame UV slope (β UV ), which is β UV = −2.4± 0.3 for our sample at z = 8-10 (D.Austin et al. 2024, in preparation; see also Cullen et al. 2023), to estimate the sample average escape fraction of photons following the prescription derived in Chisholm et al. (2022).This provides an average escape fraction ∼ 11%, with bluer galaxies as high as ∼25% and redder galaxies as low as 5%.We then apply this to the critical UV luminosity density required for the galaxy population to drive runaway reionisation (where ionization happens at a faster rate than recombination).We follow the prescription derived in Shull et al. (2012) to determine at what redshift this critical luminosity density is achieved.We assume an 11% escape fraction, a gas temperature of 2 × 10 4 K, and a redshift dependent clumping factor (C H ) following a power law of C H = (2.9)[(1+ z)/6] −1.1 (Shull et al. 2012).Applying this to our LFs, integrated down to M UV −10, we find that the critical luminosity density is achieved around z = 8.The faint end of the UV LF may well not extend to such faint luminosities and there has been debate regarding if it turns over at M UV ∼ −14 (Bouwens et al. 2017;Livermore et al. 2017;Atek et al. 2018;Zhang et al. 2022).In this situation, the critical luminosity density will be reached at lower redshifts of z ∼ 6-7.Repeating this calculation following the prescription described in Madau et al. (1999), Bolton & Haehnelt (2007), and Duncan & Conselice (2015), and a constant clumping factor of C H = 3, we find the critical density is often 40%-50% lower than with Shull et al. (2012), leading to a slightly earlier reionisation.This is a population average result, as individual systems will likely begin reionising their surroundings at earlier times.We will consider reionisation in more detail using our sample in future work.

Summary and Conclusions
In this paper we present a detailed analysis of the first results of searching for and characterizing distant galaxies with the JWST using both PEARLS GTO and public ERO data.Our total sample, which we call the EPOCHS sample, spans an area of 180 arcmin 2 and reaches near-infrared depths greater than m ∼ 28.5, forming one of the largest JWST-based search for high-redshift galaxies to date.Applying a set of strict selection criteria to the data, we identify over 1000 candidate highredshift galaxies (z > 6.5).
We use this sample to measure the UV LF of galaxies at z > 7.5.Other papers in this series will examine the individual properties of these galaxies, including their star formation rates and their stellar masses as well as morphologies and structures.Our main conclusions can be summarized as follows.
1. We find that the UV LF measured using JWST at z ∼ 8 and z ∼ 9 agrees well with past observations with HST, though there is a notable excess of z = 8 galaxies in the NEP-TDF part of the PEARLS survey.This shows that there is no significant tension between the results obtained from the photometric selection of high-z galaxies in the wavelength/redshift regime where both telescopes can independently and comfortably find highredshift systems.2. We observe a moderate number density of z = 10 sources that lies between early JWST results and past HST results.This may indicate that early JWST results which were reliant on small volumes or were dominated by the CEERS survey area may have been subject to cosmic variance effects (see also Willott et al. 2023).3. We find that the measured number density of sources at z ∼ 12 matches recently reported values using JWST, but lies on the upper end of the range predicted by a selection of different simulations.This indicates that there may be more UV-luminous sources at this epoch than what models predict.4. We find that the SFRD matches well with measurements from past studies covering the same redshift range as ours, similar to the LF.The z = 10.5 SFRD is a perfect match to expectations from a constant star-forming efficiency.Our UV LFs find a comparable SFRD compared to past studies at z > 12 and it remains slightly higher than expectations from a constant star formation efficiency.However, it must be noted that sample sizes at such high redshifts remain very small, despite the collated nature of our EPOCHS sample.More deep observations are needed in larger areas to address this point in more detail.5. We find that simulations of the high-redshift Universe tend to agree more with each other at lower redshifts (z = 8-9), but however diverge beyond HST's redshift frontier (z = 10-12).Presently, observational constraints are insufficient to distinguish which model is more favorable.However, the initially high number densities of UV-luminous sources with M UV ∼ −20.75 at z ∼ 12.5 are up to 2 dex higher than expected from some of the models.
Future improvements to these UV LF estimations are already on the horizon.Once the COSMOS-Web and PRIMER surveys are completed, the area/depth combinations will bridge the deep, pencil-beam surveys currently conducted by JWST with the wide-area VISTA surveys that have been previously undertaken.This will pin down the shape around the "knee" of the UV LF at these high redshifts.On the faint end, the completed NGDEEP and JADES programs will provide deeper and wider blank fields that reach fainter and increase sample sizes, allowing us to constrain fully the faint-end slope (α).Additionally, the large number of gravitational lenses that are being targeted by PEARLS, UNCOVER, CANUCS, and others will provide an alternative methodology to these ultradeep fields in order to find and study the more "typical" galaxies that exist at these early times.To avoid limitations from cosmic variance, it is likely that many of these programs will need to be combined in order to provide the most complete picture of the early Universe.Beyond JWST, the LSST, Euclid, and Roman telescopes will provide shallower, but much wider-area data sets (in the 10+ square degrees regime), which will enable us to constrain properly the number density of the most The integration limits go down to M UV < −17.We also show a selection of results from other recent studies covering similar redshift ranges.Small shifts of 0.1 in redshift have been applied to aid readability.The two thin solid lines show the upper and lower bounds of the the dark matter halo evolutionary models used in Oesch et al. (2018).The blue dashed line shows the constant star formation efficiency model of Harikane et al. (2023).While some early studies have made a measurement at z > 14, these estimations were largely based on the z = 16.4 candidate in the CEERS field, which later was found to be at z = 4.9.Our 11.5 < z < 13.5 bin has been placed at its average sample redshift of z = 12.luminous systems and identify the first unobscured AGN to flare up in the early Universe.
The JADES CEERS and UNCOVER teams have now provided some of the first spectroscopic confirmation of galaxies at z > 11 (Arrabal Haro et al. 2023aHaro et al. , 2023b;;Bunker et al. 2023b;Curtis-Lake et al. 2023;Wang et al. 2023).However, knowledge regarding the success rates of these measurements is still incomplete (e.g., number of sources where spectra provide ambiguous results due to insufficient lines/breaks, the number of sources found to really be lower-z, etc.).Increasing the number of ultrahigh-z sources with spectroscopic confirmation (or rejection) will be a key step in validating early findings using photometric approaches, and will ultimately be how a robust final UV LF, and other physical properties, will be measured at these early times.
All of the raw JWST data used in this work can be accessed via MAST at doi:10.17909/5h64-g193.The authors note that some data from PEARLS are within its proprietary period at the time of writing, but these will all be made accessible over 2024.Full catalogs of sources will be made public alongside the release of C. J. Conselice et al. (2024, in preparation) in early 2024. 2 < z < 5 (22 sources) and redshifted 1.6 μm bump features of z < 1 galaxies (34 sources).
Within the GLASS field, we recover both of the initially reported high-z objects that overlap between the Naidu et al. (2022) and Castellano et al. (2022) studies.In our work, these two high-redshift candidates have solutions of z = 10.2-10.9 and z = 12.3-13.6.The study by Castellano et al. (2022) reports an additional five candidates, of which only one (GHZ4) has a matching redshift solution in our work and which meets all of our selection criteria.The remaining sources have visible breaks, but have low-z solutions of 2.4 < z < 2.8 in our catalogs, indicating these breaks are potentially Balmer breaks instead of Lyman breaks.
In the first four pointings of CEERS that were taken in mid-2022 (P1, P2, P3, and P6), we select 12 out of the full sample of 26 sources found by Finkelstein et al. (2023a).An additional six sources have photo-z solutions which agree with the CEERS team but are not selected in our final robust sample.For these six, four are fainter than our detection limits, one is masked, and one has a comparable-quality low-z solution.Recently, Finkelstein et al. (2023b) published an updated catalog of CEERS targets.This more comprehensive list contains 88 objects, of which we find 87 are present in our catalogs.A total of 29 sources are selected as robust candidates in our work and 62 sources have agreeing high-z solutions.For the 33 unselected sources, 10 are masked, 22 are fainter than our selection limits, and one source has a comparable-quality low-z solution.The CEERS team utilize a smaller aperture size of diameter 0 2 in their work, resulting in slightly increased depth at the cost of increased dependence on fewer pixels and larger PSF corrections.As a result, we find there are many sources that may be reasonable candidates, but lie at the 4σ level in our work.
The work by Harikane et al. (2023) presented a sample of 21 galaxies across the GLASS, SMACS 0723, and CEERS fields.We find matches to all objects and find only four enter our final sample, with an additional two agreeing on high-z solutions.There are six galaxies where we find better fits with a Balmer break and seven galaxies where a 1.6 μm bump at z < 1 can be fit.Within the GLASS field, only the luminous objects at z = 10 and z = 12 and the additional source GHZ4 in Castellano et al. (2022) overlap between our catalogs.Following their original naming system, the two high-z galaxies which we do not consider robust enough to enter our sample are GL-z9-2 and GL-z9-9.GL-z9-2 is rejected because of a very broad redshift PDF in our work.GL-z9-9 is found to fall below our detection limits.
The work by Donnan et al. (2023) also covered the fields of GLASS, SMACS 0723, and CEERS.Of the 43 sources crossmatched to our sample, we find 28 sources have agreeing high-z solutions and 16 of these reach our final sample.There are three sources with 1.6 μm bump solutions at z < 1 and six sources with Balmer Break solutions at 2 < z < 5. Of the high-z sources that were later rejected as unreliable, there were three that were masked, two had a very broad redshift PDF, six were fainter than our detection limits, and one had a comparablequality low-z solution.
Comparing to Bouwens et al. (2023)ʼs independent analysis of SMACS 0723, GLASS, and CEERS, we find 20/33 objects have photo-z solutions that agree but only 6/33 objects enter our final sample.A total of four objects are better fit with a 1.6 μm bump at z < 1 and nine objects are fit by a Balmer break at 1.5 < z < 4. Of the objects with agreeing primary photo-z solutions, we find two are masked, 10 sources are z ∼ 6.5 with the Lyman break not completely exiting the bluest (F090W) band we utilize in SMACS 0723 and GLASS, one object is fainter than our limits, and two objects have comparablequality low-z solutions.
Finally, we compare to a similar work conducted in McLeod et al. (2024), who also compiled a variety of JWST survey programs.Of the 32 objects which were crossmatched between our samples, we find 28 have matching photo-z solutions and 24 reach our final sample.Since McLeod et al. (2024) implement conservative 7σ detection limits on their sample, we find that the rate of objects scattering below our own 5σ limits is much lower than in other studies.This ultimately means the sample presented in McLeod et al. (2024) is largely recovered in our own sample.We find one object is better fit by a 1.6 μm bump at z < 1 (their GLASS-Z11-1481) and two objects are better fit by a Balmer break at 2 < z < 4 (their CEERS-NE-z10-3069854 and CEERS-2-3-z12-8536).The remaining object we disagree upon is CEERS-SW-z11-26592, which we identify at z = 9.5 instead of z = 11.4.For the four rejected objects where we agree on the primary photo-z solutions, we do not select them because all four fall within our masked regions.Of these four, only one would be rejected for another reason (GLASS-Z14-33570 has a comparable low-z solution), indicating the other three sources may be valid.

Appendix B Spectral Energy Distributions of Highlight Sources
Within this section we show in Figure 9 the best-fitting SEDs for the UV-luminous objects described in Section 3.3.1.Alongside these we also present the SED of the highestredshift, spectroscopically confirmed source in the sample from the JADES field (Curtis-Lake et al. 2023).In addition, we also present a selection of two random candidate galaxies from each redshift bin used in measuring the UV LF in Figure 10.The full list of sources, including their coordinates, fluxes, and a selection of other properties will be presented in full as a downloadable catalog within C. J. Conselice et al. (2024, in preparation).).All SEDs show in orange the best-fit template, in blue the best-fit low-z template (z < 6), the significance of the detection just above the primary x-axis, 1σ upper limits if the measurement is less significant, 1″ cutouts in each filter used, and the accompanying photo-z PDF/statistics.In green, we show an example of an attempt to fit a brown dwarf with a JADES galaxy spectrum.
Figure 10.As with Figure 9, but for two random galaxies from each of our four redshift bins used for measuring the UV LF.The top row is z ∼ 8, and the following rows are z ∼ 9, z ∼ 10.5, and z ∼ 12.5, respectively.
ensure the majority of the redshift probability density function (PDF) is located within the primary peak.4.

Figure 1 .
Figure 1.Diagnostic plots showing the performance of our EAZY photometric-redshift procedure when applied to our NIRCam photometry.Data from a CEERS Director's Discretionary Time (DDT) proposal in the CEERS field (Arrabal Haro et al. 2023a) is shown as well as other CEERS sources from Tang et al. (2023) and Arrabal Haro et al. (2023b).The CEERS data include HST F606W and F814W imaging due to its lack coverage below 1 μm.Despite being largely limited to NIR wavelength range, photo-zs generally perform well, with the greatest spread around 2 < z < 3. The blue line shows the ideal case of the photometric redshifts matching the spectroscopic redshifts, while the red lines show 15% offsets in 1 + z, often used to describe the performance of photometric redshifts.

Figure 2 .
Figure 2.An example of the selection completeness simulation run using the JAGUAR semianalytical model with the survey setup of JADES.The color bar indicates the fractional completeness recovered and the contamination fraction of z < 5 objects entering the sample.This forms the completeness factor C i which is used in Equation (1).Bins with contamination above 10% or completeness less than 50% are not used in the measurement of the UV LF as they are considered unreliable.Other fields show similar results but shifted to brighter luminosities on account of their shallower depths.
For the MCMC, we show the median and standard error obtained from the posterior.The individual MCMC sample with the highest probability matches the best-fitting parameters of the LM methodology.For the z = 8 redshift bin, our results do not probe as faint as some previous HST studies, we thus include the faintest three data points fromMcLure et al. (2013) to improve the constraints on the faint end.Values presented with an asterisk indicate the parameter is fixed to this value in the fitting procedure.McLure et al. (2013), Finkelstein et al. (2015), and Bouwens (PID: 1283, PIs: H. Nørgaard-Nielsen and G. Östlin) parallel observations conducted by Pérez-González et al. (2023).Most recently, Finkelstein et al. (2023b) conducted a full search of high-z galaxies in the CEERS field.We find close agreement with their results, with only slightly lower number densities in our observations across the overlapping luminosity range.4.6.The Measured z = 10 and z = 12 Ultraviolet Luminosity Function

Figure 6 .
Figure 6.The UV LF at 11.5 < z < 13.5 as measured in this study (large red points).We compare our measurements to a selection of recent observations (left) and predictions from simulations (right).We display here the best-fitting Sch and DPL functional forms but emphasize that uncertainties on these fits are large.Observational data are from Finkelstein et al. (2023a, 2023b), Bouwens et al. (2023), Donnan et al. (2023), Harikane et al. (2023), and Pérez-González et al. (2023).
(2023) and much steeper than the derived slope from McLeod et al. (2024; −0.23 ± 0.04).If we compare our observations to the predicted SFRD of Harikane et al. (2023), which is based on a constant star formation efficiency, we find our results up to z = 12 are consistent with such an evolution.Tension with such a model presently arises if the z ∼ 15-16 samples of galaxies are introduced.

Figure 7 .
Figure 7.The measured time evolution of the UV LF.Here, we present the best-fitting DPL functions for our four redshift bins alongside each other.We find that the bright end of the UV LF has relatively minor to negligible evolution within our errors (agreeing with Bowler et al. 2020), while the faint end decreases in normalization toward higher redshifts.

Figure 8 .
Figure8.The measured UV luminosity density and associated SFRD, calculated by integrating UV LFs generated by randomly sampling the MCMC distribution for the DPL fits.The integration limits go down to M UV < −17.We also show a selection of results from other recent studies covering similar redshift ranges.Small shifts of 0.1 in redshift have been applied to aid readability.The two thin solid lines show the upper and lower bounds of the the dark matter halo evolutionary models used inOesch et al. (2018).The blue dashed line shows the constant star formation efficiency model ofHarikane et al. (2023).While some early studies have made a measurement at z > 14, these estimations were largely based on the z = 16.4 candidate in the CEERS field, which later was found to be at z = 4.9.Our 11.5 < z < 13.5 bin has been placed at its average sample redshift of z = 12.

Figure 9 .
Figure 9. SEDs for a selection of sources noted in this work.The naming of NEP-TDF sources identifies which spoke number (epoch) in which the object was identified.In this case, all NEP sources presented are from the third spoke (NEP-3) taken in 2023 February from PID 2738.Top left: the source with the secondbrightest intrinsic M UV in our sample, after the AGN identified in Larson et al. (2023a), with M UV = −21.34.Top right: the candidate source at z > 10 with the brightest M UV .This is at z = 11.5 and has M UV = −21.03.Bottom left: an example of one of the five bright sources (M UV < −21) at 8.1 < z < 8.3 in the NEP-TDF.Bottom right: the SED of the highest-redshift spectroscopic source in the sample from the JADES team (Curtis-Lake et al. 2023).All SEDs show in orange the best-fit template, in blue the best-fit low-z template (z < 6), the significance of the detection just above the primary x-axis, 1σ upper limits if the measurement is less significant, 1″ cutouts in each filter used, and the accompanying photo-z PDF/statistics.In green, we show an example of an attempt to fit a brown dwarf with a JADES galaxy spectrum.

Table 1
Average 5σ Depths for Point Sources in 0 32 Diameter Apertures

Table 2
The Overlap of Sources in Our Final Catalogs Compared to Those of Other Studies that Have Analyzed the Same Fields

Table 4
The Best-fit Parameters for the Schechter and Double Power-law Functional Forms for Both the Levenberg-Marquard Approach and the Markov Chain Monte Carlo Approach