The galaxy UV luminosity function at z ≃ 11 from a suite of public JWST ERS, ERO, and Cycle-1 programs

ABSTRACT

We present a new determination of the evolving galaxy ultraviolet (UV) luminosity function (LF) over the redshift range 9.5 < z < 12.5 based on a wide-area (>250 arcmin²) data set of JWST NIRCam near-infrared imaging assembled from 13 public JWST surveys. Our relatively large-area search allows us to uncover a sample of 61 robust z > 9.5 candidates detected at ≥8σ, and hence place new constraints on the intermediate-to-bright end of the UV LF. When combined with our previous JWST + UltraVISTA results, this allows us to measure the form of the LF over a luminosity range corresponding to four magnitudes (M₁₅₀₀). At these early times we find that the galaxy UV LF is best described by a double power-law function, consistent with results obtained from recent ground-based and early JWST studies at similar redshifts. Our measurements provide further evidence for a relative lack of evolution at the bright-end of the UV LF at z = 9–11, but do favour a steep faint-end slope (α ≤ −2). The luminosity-weighted integral of our evolving UV LF provides further evidence for a gradual smooth (exponential) decline in co-moving star-formation rate density (ρ_SFR) at least out to z ≃ 12, with our determination of ρ_SFR(z = 11) lying significantly above the predictions of many theoretical models of galaxy evolution.

1 INTRODUCTION

Although JWST has not yet completed its first year of science operations, it is already transforming our view of the young Universe. Despite its many achievements, the Hubble Space Telescope (HST) was unable to advance the search for early galaxies significantly beyond redshifts z ≃ 10, due to its limited near-infrared wavelength coverage (λ < 1.6 μm; e.g. Coe et al. 2013; Ellis et al. 2013; Oesch et al. 2016). By contrast, the exquisite near/mid-infrared imaging now being provided by the NIRCam instrument on-board the larger and colder JWST has already pushed the redshift frontier out to z ≃ 13 (Finkelstein et al. 2022b; Naidu et al. 2022; Robertson et al. 2023), with candidate high-redshift galaxies already being uncovered at redshifts as extreme as z = 16–17 (Harikane et al. 2023b).

The photometric selection of high-redshift galaxy candidates from deep imaging not only supplies targets for spectroscopic follow-up to determine their astrophysical properties, but also enables the statistical study of the evolution of the galaxy population provided sufficiently large and robust samples can be assembled. In particular, a key goal is to chart the evolution of the galaxy luminosity function (LF) back to early times for comparison with theoretical predictions. At very high redshifts, near/mid-infrared imaging enables the evolving rest-frame ultraviolet (UV) galaxy LF to be determined, the luminosity-weighted integral of which gives the evolution of (co-moving) UV luminosity density, which can then be converted into the evolution of cosmic star-formation rate density, ρ_SFR. In large part due to the impact of HST, there is now a general consensus over the evolution of ρ_SFR out to z ≃ 8. However, before the advent of JWST, the lack of significant galaxy samples at higher redshifts resulted in considerable disagreement over the evolution of star-formation rate density at earlier times. Specifically, some studies concluded in favour of a rapidly steepening decline in ρ_SFR with increasing look-back time beyond z ≃ 8 (e.g. Oesch et al. 2013, 2018), while others presented evidence for a continued smooth, more gradual evolution implying the existence of substantial star-formation activity at still earlier times (McLeod et al. 2015; McLeod, McLure & Dunlop 2016).

This debate/disagreement over the very high-redshift evolution of the UV LF, and hence ρ_SFR, was largely a consequence of the limitations of pre-JWST facilities mentioned above. Specifically, at redshifts z ≥ 9 the WFC3/IR camera on HST rapidly runs out of filters capable of sampling continuum emission beyond the Lyman break at λ_rest = 1216 Å. For most of the deep extragalactic imaging surveys performed with HST only F125W and F160W imaging are available for this purpose, at least for the majority of the CANDELS (Grogin et al. 2011) area. The addition of the F140W filter for such surveys as the HUDF (Beckwith et al. 2006; Ellis et al. 2013), CLASH (Postman et al. 2012), and the Frontier Fields (Lotz et al. 2017) was key for selecting z = 9 galaxies in McLeod et al. (2015), McLeod, McLure & Dunlop (2016) as it provided a third filter with high-resolution imaging that probes beyond the Lyman break. The Spitzer space telescope (e.g. S-CANDELS; Ashby et al. 2015) and the deepest ground-based K_s −band imaging [e.g. Very Large Telescope (VLT)’s HUGS; Fontana et al. 2014] helped to extend the wavelength baseline beyond λ = 1.6 μm. Using the additional wavelength coverage of Spitzer, studies such as Roberts-Borsani et al. (2016) and Finkelstein et al. (2022a) provided constraints on the rest-frame optical continuum of some brighter high-redshift galaxies. The increased constraints on the rest-frame optical from Spitzer also allows more accurate stellar mass estimates at high-redshift, as demonstrated by Stefanon et al. (2021b) who were able to constrain the galaxy stellar mass function at z = 6–10 by leveraging HST-selected samples with exceptionally deep IRAC data from the GREATS program (Stefanon et al. 2021a). While Spitzer had been incredibly useful for expanding our understanding of high-redshift galaxy properties, particularly for bright galaxies, this redder non-HST imaging lacked the resolution and the depth to benefit searches for fainter galaxies at z > 8.

JWST has effectively unlocked the longer wavelength near/mid-infrared coverage previously only provided by ground-based facilities and Spitzer, but now with enormous improvements in both angular resolution and depth. Several early, public JWST NIRCam surveys have reached 5σ limiting depths of 28.5–29.0 mag in filters spanning the wavelength range λ ≃ 1–5 μm, already probing several magnitudes deeper than typical Spitzer surveys, and with diffraction-limited imaging down to λ ≃ 1 μm delivering near-infrared imaging of much improved resolution compared to HST’s WFC3/IR. This combination of depth and resolution has already been transformative in the selection of ultra-high redshift galaxies. Encouragingly, many of the z > 10 candidates photometrically selected from the early NIRCam imaging are already being confirmed spectroscopically with JWST NIRSpec (Arrabal Haro et al. 2023a, b; Harikane et al. 2023a).

While several theoretical models of galaxy evolution appear to predict a rapid decline in ρ_SFR at z > 8 (e.g. Mason, Trenti & Treu 2015; Yung et al. 2019), early JWST studies have provided evidence to the contrary. Donnan et al. (2023a) conducted a wide-area ground-based search for bright z = 8–10 galaxy candidates using COSMOS/UltraVISTA DR5 (McCracken et al. 2012), in combination with a search for fainter galaxies at these redshifts with the JWST CEERS (Finkelstein et al. 2023), GLASS-parallel (Treu et al. 2022) and SMACS0723-73 (Pontoppidan et al. 2022) public surveys. The resulting UV LF determination was found to be consistent with the earlier HST-based works of McLeod et al. (2015) and McLeod, McLure & Dunlop (2016), with ρ_SFR displaying an extended smooth decline at z > 8 and indeed continuing to display a log-linear relation with redshift until at least z ∼ 12. At the bright-end, the number densities of z = 8–10 galaxy candidates found using the latest COSMOS/UltraVISTA ground-based data suggested that the UV LF follows a double power-law functional form, as previously found by Bowler et al. (2020). Although the exceptionally bright z = 16 candidate CEERS-93316 first reported in Donnan et al. (2023a) has now been spectroscopically confirmed to be a lower redshift interloper at z = 4.9 (Arrabal Haro et al. 2023a), this is a particularly pathological case in which a peculiar combination of very strong emission lines conspired to masquerade as the blue continuum above the Lyman-break. Another potential z = 16 candidate behind Stephan’s Quintet may yet suggest that significant star-formation activity is still taking place at z ≥ 14 (Harikane et al. 2023b).

As well as the discovery of extremely high-redshift galaxies at z > 12 (Naidu et al. 2022; Atek et al. 2023; Donnan et al. 2023a; Harikane et al. 2023b), JWST has been yielding surprisingly large numbers of bright z ≃ 10 galaxies, with a particularly high number density reported by Castellano et al. (2023) in the area around the Abell 2744 cluster. Such densities are potentially up to a factor ten times greater than expectations based on previous observations and theoretical models of galaxy evolution. This would challenge our existing understanding, with many existing theoretical models of galaxy evolution having to invoke more efficient star formation or less dust attenuation at the highest redshifts to match these observational results. Alternatively, these early observations may be biased to a particular population of young and rapidly star-forming galaxies (Mason, Trenti & Treu 2023). There is thus an open question as to how accurately the early JWST-selected galaxy samples represent the true number densities of bright z = 10 galaxies, with existing high-redshift galaxy searches typically restricted to small areas.

To address this issue, the study presented here aims to exploit all relevant available JWST NIRCam Early Release Observations (ERO), Early Release Science (ERS), and public Cycle-1 programs to conduct a systematic wide-area search for extreme-redshift galaxy candidates with high signal-to-noise (≥8σ). The resulting combined survey, covering a raw area of ≃ 260 arcmin², provides a new opportunity to uncover statistically significant samples of galaxies in the redshift range 9.5 < z < 12.5, and hence yield improved constraints on the number densities of bright to intermediate M₁₅₀₀ galaxies over this redshift range. By combining the new results with our previously derived constraints on the number densities of fainter galaxies (Donnan et al. 2023a), we are hence able to determine the evolving UV LF and ρ_SFR at z = 10–12, deep into the redshift range where the results derived using HST had diverged. Moreover, our wide-area search may also uncover other examples of extremely high-redshift (z > 12) galaxy candidates, such as those found in the ERO surveys (e.g. Atek et al. 2023; Harikane et al. 2023b).

The structure of this paper is as follows. In Section 2 we describe each of the data sets included in our high-redshift galaxy search. Section 3 describes the construction of our multiwavelength catalogues and candidate selection. We discuss our candidates and draw comparisons with other recent JWST searches in the literature in Section 4. In Section 5, we use our high-redshift galaxy sample to determine the form of the UV LF over the redshift range 9.5 < z < 12.5 and how it evolves from z ≃ 8 to z ≃ 12. In Section 6 we integrate the UV LF in order to determine the star-formation rate density and how it evolves with cosmic time. Finally, in Section 7 we summarize our conclusions. Throughout the paper, we adopt a cosmology with Ω₀ = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹Mpc⁻¹, and use the AB magnitude system (Oke 1974; Oke & Gunn 1983).

2 DATA

2.1 Imaging data

We begin with a description of the various data sets employed in this study. Each of the data sets were downloaded as ‘rate’ products, either from the Canadian Astronomy Data Centre (CADC) or STScI MAST data base, and processed through the Primer Enhanced NIRCam Image Processing Library (pencil; Magee in preparation) software, which is a customized version of the JWST pipeline version 1.6.2. This version of the pipeline includes additional routines for the removal of snowballs and wisps, as well as background subtraction and correcting for 1/f noise striping. Given the different timeframes in which the various data sets were taken, released, reduced, and analysed, there are very slight differences in the CRDS context depending on when the data set was released, although the pmap is always at least pmap 0989, which included the most recent zero-point corrections at the time of writing. Where applicable, HST data sets were downloaded as their high-level science products from the STScI MAST, aligned to the JWST imaging and registered using swarp (Bertin et al. 2002). The exception is for CEERS, where HST EGS imaging in F606W and F814W was provided in a data release by the CEERS team (see Koekemoer et al. 2011).

The global 5σ limiting magnitudes for each of the different surveys and in each filter can be found in Table 1. These global depths were measured in 0.35 arcsec-diameter apertures and corrected to total assuming a point-source correction, following the method outlined in Donnan et al. (2023a). Note that the global depth in a given survey is subject to potentially significant spatial variations in local depth, owing to differences in exposure time, source-crowding (particularly for the Quintet, Cartwheel, and cluster fields) or sensitivities of the NIRCam modules. We note that these quoted depths are somewhat shallower than those of other studies who use smaller apertures, e.g. 0.2 arcsec-diameter such as Bagley et al. (2023a) in the data release for CEERS. While the adoption of larger apertures is a more conservative approach, it is still instructive to verify that our reductions are of a similar quality and depth to those used in other studies. We therefore tested that using an equivalent aperture size would yield similarly deep limiting magnitudes. Our resulting depths were consistent to within ≃ 0.1 mag of the quoted depths by Bagley et al. (2023a). We also include the effective survey area for each of the fields, accounting for de-lensing where applicable.

We illustrate the distribution in depths across the various surveys in Fig. 1. Assuming a flat UV spectral slope (f_ν), we derive a rest-frame M₁₅₀₀ limiting magnitude for each survey, corresponding to a z = 11 galaxy at the 8σ apparent magnitude limit in the F200W filter. Using the UV LF determined in Donnan et al. (2023a), we also show the predicted average number of z ∼ 11 galaxies found per NIRCam pointing (≃ 9.6 arcmin²).

2.1.1 SMACS 0723

In this study, we cover three ERO (Pontoppidan et al. 2022) surveys: the SMACS 0723 cluster, Stephan’s Quintet and the Cartwheel galaxy. All three of these surveys are comprised of imaging in six NIRCam filters spanning |$0.9-4.4\ \mu$|m: F090W, F150W, F200W, F277W, F356W, and F444W.

The first ERO field we will discuss is the SMACS 0723 cluster field (PID 2736, PI Pontoppidan; Pontoppidan et al. 2022). This field has already seen extensive study by various groups in search of high-redshift candidates since the 2022 July data release (e.g. Adams et al. 2023b; Atek et al. 2023; Donnan et al. 2023a; Bouwens et al. 2023; Harikane et al. 2023b). This survey is covered by two NIRCam modules, with one module centred on the cluster field and the other on a ‘blank’ parallel field. Due to lensing by the foreground galaxy cluster, which serves to magnify objects in the background, there is also a reduction in the effective area of the search.

To avoid uncertain lensing volume corrections and areas with significantly shallower depths than the global depth, we remove the portion of the image worst affected by intracluster light and enhanced crowding. We then determine the magnification of sources and the remaining effective survey area by utilizing the lensing model from Pascale et al. (2022). Although the raw survey area of the SMACS 0723 field is 11.5 arcmin², after excluding the cluster and accounting for lensing this is reduced to an effective area of 6.4 arcmin². This area is calculated with the assumption that any part of the field-of-view that is not covered by a gravitational lensing map has a magnification μ = 1.2, which is what we adopt for all cluster + parallel pair surveys to take into account that the ‘blank’ field is still subject to magnification (see Section 5).

Over the cluster portion of the field there is ancillary HST data from the RELICS survey (PI Coe; Coe et al. 2019), with ACS imaging in F435W, F606W, and F814W as well as WFC3/IR imaging in F105W, F125W, F140W, and F160W. We include the HST ACS data in order to provide wavelength coverage in the optical portion of the spectrum. Although the 5σ limiting magnitudes of the HST ACS data are significantly shallower than the JWST data, the data still serves to verify any non-detections short-ward of the Lyman break. The 5σ depths (0.35 arcsec-diameter apertures corrected to total) of the ACS F435W and F814W imaging are approximately 26.8 mag, with the F606W deeper at 27.7 mag.

2.1.2 Stephan’s Quintet

The Stephan’s Quintet field was also imaged as part of the ERO (PID 2732; PI Pontoppidan; Pontoppidan et al. 2022). To date, this field has had one published search for high-redshift galaxies by Harikane et al. (2023b), which yielded one galaxy candidate at z ≃ 12 and another candidate at z ≃ 16. This field has the same NIRCam filter set as SMACS 0723, but lacks sufficiently deep ancillary HST data. The Quintet field has significant variations in depth across the imaging owing to crowding and enhanced background light. The overall area of the field is approximately 46.7 arcmin², which reduces to an effective area of 40.9 arcmin² after removing the Quintet.

2.1.3 Cartwheel galaxy

The Cartwheel galaxy, or PGC 2248, was imaged as part of the ERO (ERO PID 2727; PI Pontoppidan). Once again, the complement of NIRCam filters for this field is the same as for SMACS 0723 and Stephan’s Quintet, and there is no useful ancillary HST ACS imaging. Similarly to the Quintet field, the depth across the imaging varies significantly, even when masking out the Cartwheel galaxy itself. The field spans ∼5.7 arcmin², which reduces to 4.2 arcmin² after masking out the Cartwheel galaxy and the two large foreground galaxies to the north-east (NE).

2.1.4 WHL 0137

The WHL0137 cluster, known for the lensed ‘sunrise arc’ galaxy and the z = 6.2 star candidate ‘Earendel’, was imaged as part of Cycle 1 GO PID 2282 (PI Coe; see Welch et al. 2022; Bradley et al. 2023). Imaging was taken over a single NIRCam pointing, with one module on the cluster and the other on a parallel blank field. The filter set includes the south-west (SW) channel F090W, F115W, F150W, and F200W bands, as well as the LW channel F277W, F356W, F410M, and F444W bands. The raw survey area of this field is approximately 10.5 arcmin², which reduces to an effective area of 5.9 arcmin² when accounting for the gravitational lensing. This data set has been the subject of a previous search for high-redshift galaxies by Bradley et al. (2023), which yielded a sample of four z = 8–10 candidates. Our current study covers the first epoch of observations from 2022 July.

As with SMACS 0723, this field has ancillary HST RELICS (Coe et al. 2019) data over the cluster region, which we use both in spectral energy distribution (SED) fitting and to check against any potential low-redshift interlopers in our final sample. The 5σ depths (0.35 arcsec-diameter apertures, corrected to total) of the ACS F435W and F814W imaging are approximately 26.7 mag, with the F606W imaging significantly deeper at 27.4 mag.

The lensing maps we utilize for WHL 0137 were those prepared for the RELICS survey: Zitrin-LTM (see Broadhurst et al. 2005; Zitrin et al. 2015a), glafic (Oguri 2010), and lenstool (Jullo et al. 2007; Jullo & Kneib 2009).

2.1.5 MACS 0647

The MACS 0647 galaxy cluster was previously one of the HST CLASH clusters (Postman et al. 2012). The field in which MACS 0647 resides is known to contain a triply lensed z = 11 galaxy candidate, MACS0647-JD1 (Coe et al. 2013), which was one of the highest-redshift galaxy candidates in the pre-JWST era. MACS 0647 was imaged as part of JWST GO proposal 1433 (PI Coe) in order to follow-up this object, with a cluster + parallel pair of NIRCam modules imaged in F115W, F150W, F200W, F277W, F356W, and F444W. With this imaging, taken in 2022 September, the object was re-analysed and found to be a potential merger system at z ∼ 10 (Hsiao et al. 2023). Most recently, the object was spectroscopically confirmed to lie at z = 10.17 (Harikane et al. 2023a).

To complement the JWST photometry, we also employ the HST ACS optical imaging from the CLASH survey, which covers the cluster module and has a 5σ depth of ≃27.0 mag. The raw area of the MACS 0647 field was found to be 10.5 arcmin², which reduces to 5.1 arcmin² after de-lensing, and the removal of un-useable regions due to the presence of particularly bright stars.

We use the lensing maps that were produced as part of CLASH, i.e. the Zitrin-NFW and Zitrin-LTM-Gauss models (Zitrin et al. 2015a).

2.1.6 RXJ 2129

Another HST CLASH cluster, the RXJ 2129 field was imaged in 2022 October as part of Director’s Discretionary Time (DDT, PID 2767; PI Kelly). One NIRCam module was imaged in F115W, F150W, F200W, F277W, F356W, and F444W, although the exposure times vary significantly between filters, with the F356W imaging being particularly deep. The complement of HST imaging matches that of MACS 0647. The lensing maps also include Zitrin-LTM-Gauss and Zitrin-NFW but with an additional model from Caminha et al. (2019). This field was initially selected as a DDT proposal to follow up a supernova at z = 1.5, but has also yielded a very faint, yet highly magnified (μ ∼ 20) galaxy confirmed with spectroscopy to be at z = 9.51 (Williams et al. 2023b). The raw area of this field is 5.3 arcmin², with a de-lensed effective area of 3.1 arcmin².

2.1.7 GLASS

As part of JWST PID 1324 (PI Treu; Treu et al. 2022; see also Paris et al. 2023), a pair of NIRCam pointings were taken in parallel to NIRISS observations of the Abell 2744 cluster field. The first epoch of observations was taken in 2022 July, and has been the subject of numerous searches for high-redshift galaxies, including the discovery of an exceptionally bright z = 12 galaxy candidate GLASS-z12 (Castellano et al. 2022; Naidu et al. 2022). The second epoch of observations was taken in 2022 November, and has yielded yet further exciting results, including the discovery of an over-density of bright z = 10 galaxy candidates (Castellano et al. 2023).

We include both epochs of observations in our present study. The NIRCam filter set comprises F090W, F115W, F150W, F200W, F277W, F356W, and F444W. One NIRCam module benefits from increased exposure time from repeated observation, and there is an additional two NIRCam modules with significant overlap, producing a total raw search area of 13.3 arcmin².

Some of the GLASS parallel field is covered by the Furtak et al. (2023) lensing model of the Abell 2744 cluster, which was released as part of UNCOVER Data Release 1. Where there is coverage we adopt this lensing map, and set a floor of μ = 1.2 for the rest of the area that is not covered. This gives us an effective de-lensed area of 10.1 arcmin².

2.1.8 North Ecliptic Pole Time Domain Field

The North Ecliptic Pole Time Domain Field (NEP TDF) field was imaged in 2022 September as part of the PEARLS survey (PID 2738, PI Windhorst; see Windhorst et al. 2023). NIRCam imaging was taken in F090W, F115W, F150W, F200W, F277W, F356W, F410M, and F444W. The majority of this data set is proprietary, however one of the pointings in the NEP TDF was made publicly available immediately for use by the community. Although the footprint of the F150W and F356W imaging is larger, we focus on the ≃ 10.5 arcmin² overlap between all of the available NIRCam filters.

2.1.9 Cosmic Evolution Early Release Science Survey

The Cosmic Evolution Early Release Science Survey [CEERS; ERS 1345, PI Finkelstein; see Bagley et al. (2023a) for imaging release] is an early release science program covering approximately 100 arcmin² of the CANDELS Extended Groth Strip (EGS) field. This field has been studied extensively since the first epoch of observations was taken in 2022 July (e.g. Finkelstein et al. 2022b, 2023 ; Bouwens et al. 2023; Donnan et al. 2023a). This study has also already seen the spectroscopic confirmation of several high-redshift galaxies at z = 8–11, including ‘Maisie’s galaxy’ (Finkelstein et al. 2022b; Arrabal Haro et al. 2023a, b).

In this study, we include an analysis of both the first and second epochs of observations, taken in 2022 July and December, respectively. JWST imaging is available in seven NIRCam bands: F115W, F150W, F200W, F277W, F356W, F410M, and F444W. We also employ HST ACS imaging to provide additional wavelength coverage in the optical. As noted earlier, a reduction of the HST imaging in ACS F606W and F814W has been made publicly available by the CEERS team (HDR1; see Koekemoer et al. 2011). Additional F435W imaging was available from the UVCANDELS survey (GO 15647, PI Teplitz; DOI: 10.17909/8s31-f778).

2.1.10 UNCOVER

The UNCOVER survey (Cycle 1 2561; PI Labbe; Bezanson et al. 2022) is an ultra-deep survey of the Abell 2744 cluster field covering a raw area of ≃ 29 arcmin² (≃ 12 arcmin² after delensing and masking) to a depth in excess of 29 mag (0.35 arcsec-diameter apertures). The NIRCam filter set is F115W, F150W, F200W, F277W, F356W, F410M, and F444W. The lensing map used for this data set was the Furtak et al. (2023) model that was publicly released as part of UNCOVER’s Data Release 1. We also employ this lensing map for the regions of GLASS and DDT 2756 where there is coverage.

Abell 2744 was imaged with HST as one of the Hubble Frontier Fields (Lotz et al. 2017). This survey includes among the deepest HST observations in the ACS and WFC3/IR filters, with depths of 28.2–28.6 mag. Although the ACS observations cover less than half of the UNCOVER survey area, we include the ACS F435W, F606W, and F814W data in order to extend our wavelength coverage into the optical.

2.1.11 DDT 2756: Abell2744

The Abell 2744 cluster was imaged with NIRCam as part of DDT program 2756 (PI Chen), a program designed to follow-up a z = 3.47 supernova. Two epochs of observations in six NIRCam bands (F115W, F150W, F200W, F277W, F356W, and F444W) were taken in 2022 October and December. Although the Abell 2744 module imaging was included in the UNCOVER reduction, each epoch includes an additional parallel NIRCam module. These two partially overlap, but together provide another 6.6 arcmin² of effective search area. This field was also found by Castellano et al. (2023) to contain a very bright z ≃ 10 galaxy candidate, with reported M_UV ≃ −21.60. Ultimately, the object has been spectroscopically confirmed to lie at z = 9.31, and is likely an interacting system (Boyett et al. 2023). That being said, an independent search over both epochs of data for this survey is clearly desirable to discover if there are other exceptionally bright candidates at similarly extreme redshifts.

2.1.12 PID 1063 Nircam Flats: J1235

As part of NIRCam commissioning (PID 1063, PI Sunnquist), exceptionally deep imaging was taken of a field centred at RA = 12^h35^m48^s, Dec = +4°55′45′. The filter set included is also more extensive than the other JWST surveys: F070W, F090W, F115W, F150W, F200W, F277W, F300M, F356W, and F444W. The coverage across this rich filter set is not homogeneous, and so we restrict our analysis to an area of ≃ 9.9 arcmin² covered by all of the available filters. This field was recently studied as part of a search for quiescent galaxies by Valentino et al. (2023).

2.1.13 JEMS

Deep medium band observations of the HUDF were taken as part of the JEMS survey (ID 1963; PI Williams; see Williams et al. 2023a). This included SW F182M and F210M imaging along with LW imaging in F430M, F460M, and F480M. The total 5σ global depths are ≃ 29 mag for the SW channel filters and ≃ 28.5 mag for the LW channel filters.

Using this data set, Donnan et al. (2023b) uncovered a sample of six high-redshift candidates, including three spectroscopically confirmed z > 10 galaxies (Curtis-Lake et al. 2023; Robertson et al. 2023). Among these is UDFj-39546284, a z ≃ 12 galaxy first discovered a decade previously with HST imaging (Bouwens et al. 2011; Ellis et al. 2013). Although we do not search the JEMS data set in this study, we include the six candidates found in Donnan et al. (2023b) in our LF analysis in Section 5.

3 SAMPLE SELECTION

3.1 Catalogue construction

For each of the fields, we construct multiwavelength catalogues using source-extractor (Bertin & Arnouts 1996) in dual-image mode. As our study is focused on selecting 9.5 < z < 12.5 galaxies, the position of the Lyman break will move through the F150W filter as we probe to higher redshift, with the intergalactic medium (IGM) absorption having attenuated half of the flux in the F150W band by z = 11. Hence, in order to facilitate the detection of galaxies beyond a Lyman break of 1.5 μm, we use the F200W imaging as our primary detection band. We create additional catalogues using each of the broadband LW filters (F277W, F356W, F444W) as detection images, in order to supplement this first catalogue with any sources that may have been missed. For the three fields lacking F115W coverage (SMACS 0723, Cartwheel, Quintet), we cannot search as effectively for z < 11 galaxies and a F150W detection image would not add significant value. Hence for consistency across all of the fields, we do not include F150W detection images in this study. We also adjust our selection criteria in these three fields to account for the reduced filter coverage (see Section 3.3).

Prior to any aperture photometry, we homogenized the point spread function (PSF) of all of the NIRCam images to match that of the F444W image, in order to minimize any colour systematics in subsequent analyses. Aperture photometry was performed on these PSF-homogenized images with fluxes measured in 0.35 arcsec-diameter apertures. This aperture size is designed to contain close to 75 per cent of point-source flux for the F444W PSF. We note that this choice of aperture size is reasonable given the sizes of z = 9–10 galaxies found previously in the literature, e.g. Morishita et al. (2018) and Naidu et al. (2022). Photometric uncertainties were calculated object-by-object by taking the median absolute deviation of the nearest 150–200 blank-sky 0.35 arcsec-diameter apertures, and scaling to σ by multiplying by 1.4826, following the method by McLeod et al. (2021). Aperture fluxes across all filters were then corrected (typically at the ≃1–2 per cent level) to 75 per cent of point-source total to take into account any remaining inhomogeneities in the curves of growth between filters after PSF homogenization.

For each of our initial catalogues, we required a 5σ detection in the detection image, as well as an additional 3σ detection in at least one other band to reduce the number of spurious detections.

3.2 SED-fitting and photometric redshifts

We performed SED fitting on all of the sources within our initial catalogues using the SED fitting code LePhare (Arnouts et al. 1999; Ilbert et al. 2006). The SED library that we used was Bruzual & Charlot (2003; BC03), using a Chabrier (2003) initial mass function (IMF), a Calzetti et al. (2000) dust attenuation law, and the IGM absorption prescription from Madau (1995). These SED templates included declining star-formation histories with various τ values from 0.1 to 15 Gyr, and metallicities of 0.2 Z_⊙ and Z_⊙. We included emission lines, and the A_V was allowed to range from 0 to as high as 6 to clean the sample of any extremely reddened low-redshift interlopers masquerading as high-redshift galaxies.

We measured rest-frame M₁₅₀₀ mag using a 100 Å-wide tophat filter on the best-fitting SED. Given the remarkably high resolution of JWST’s NIRCam imaging, even the highest-redshift galaxies are often not well approximated by point-sources, and so require larger corrections to total than point-source corrections. Hence, rather than simply scaling to point-source total flux, we correct the M₁₅₀₀ mag by scaling to the FLUX_AUTO (Kron 1980) value measured by source-extractor, and we add a further 10 per cent correction to account for the light outside of the Kron aperture. This further correction is motivated by an analysis performed on stacks of our extended source z ≥ 9 candidates (see Section 4), where we found that FLUX_AUTO under-estimated the flux versus a large 1.2 arcsec-diameter aperture by |$\simeq 10{{\ \rm per \, cent}}$|⁠.

3.3 Selection of high-redshift candidates

We initially selected candidates with a photometric redshift solution z_phot ≥ 8.5 and a goodness of fit |$\chi ^{2}_{\nu }\le 10$|⁠. Only candidates considered ‘robust’, with Δχ² ≥ 4 between the primary and secondary photo-z solutions, were retained for the final sample. We also further refined our sample selection to require a ≥8σ detection in at least one of the F150W, F200W, and F277W images, to ensure that we were only analysing high signal-to-noise objects. Due to the IGM absorption of any flux short-ward of the Lyman break at z ≤ 9.5, we further required a non-detection (≤2σ) in the F090W band, as well as in F115W and the HST ACS imaging if available. The condition of F115W<2σ precludes any complete statistical investigation of the 8.5 < z < 9.5 population, given that a z ≃ 9.0 galaxy can still be detectable in F115W. Consequently, we confine our later LF analysis to z > 9.5, but include any 8.5 < z < 9.5 candidates in our tabulated sample.

As a validation exercise, we double-checked any high-redshift candidates did not present a preferred low-redshift solution when performing SED fitting using the code eazy (Brammer, van Dokkum & Coppi 2008), including the Pegase (Fioc & Rocca-Volmerange 1999) templates. The eazy photometric redshift solution was determined allowing the templates to vary over the redshift range 0 < z < 25.

As well as performing the flux cuts to remove any significant detections shortward of the break, we inspected the F090W, F115W, and ACS imaging, in order to verify that our candidates do not have any marginal detections.

The Cartwheel, SMACS 0723, and Quintet fields lack the F115W filter, and this creates additional challenges in determining the redshift of any candidates in those fields. As the F090W transmission ends at 1.0 μm, any F090W-dropout followed by a significant F150W detection can essentially lie anywhere between z ∼ 7.2 and z ∼ 12.0, assuming it is not a low-redshift interloper. This issue is mitigated where there is a F115W non-detection, which helps narrow down the redshift solution to z ≳ 9. We therefore included an extra criterion for these fields, whereby we required either a non-detection (≤2σ) in F150W, or a red F150W–F200W≥0.75 colour. Although this shifts our selection window for these three fields to above z ≃ 11.5, these criteria prevent the final sample from being contaminated by z ≃ 7–8 interlopers.

4 HIGH-REDSHIFT CANDIDATES

In Tables 2 and 3 we list basic properties for the candidates uncovered across the various fields. We also provide a colour-coded illustration of the distribution of intrinsic M₁₅₀₀ with redshift, and the distribution of F200W 0.35 arcsec-diameter aperture magnitude with redshift for our final sample in Fig. 2. In Figs 3 and 4 we present some example SEDs and postage stamps for a sub-set of the final candidates.

In Appendix  A, we provide a field-by-field analysis of the high-redshift candidates that we have selected, where we draw comparisons with other works in the literature. Although we do not search the HUDF for candidates in this particular study, we performed a search for z > 9.5 galaxies in Donnan et al. (2023b) very recently using the JEMS (Williams et al. 2023a) data. We include the six robust F182M-selected candidates found over one NIRCam pointing in our UV LF determination in Section 5.

5 UV LF

In this section, we present our determination of the UV LF at 9.5 ≤ z ≤ 12.5. We first describe our methods for assigning survey areas, effective volumes, completeness corrections, and uncertainties in galaxy number densities. We then calculate and fit the UV LF, before comparing our new results with relevant determinations from the literature.

5.1 Survey area definition

As this study is conducted over many inhomogeneous fields, care must be taken in order to define an accurate overall survey area. In particular, the fields with clusters are subject to the effects of gravitational lensing. At the faintest luminosities, the differences between lensing models has been shown to have a potentially large impact on the faint end of the UV LF, with the high magnification regime (μ > 10) more subject to systematic uncertainties (see e.g. Bouwens et al. 2017). Moreover, the cluster regions are typically significantly shallower than the survey global depth due to enhanced background from intracluster light and source crowding. For each of our cluster fields, we therefore mask out these shallow cluster regions. Any candidates positioned on these masks are excluded from the UV LF calculation.

For the cluster surveys, we subtract the clusters in a fairly simple way. We first measure the global depth σ_g across the whole field of view, as well as local depths across the survey pixel by pixel, and produce a depth map as shown in Fig. 5. We then mask out all pixels that have a local depth that is shallower than 2 × σ_g, removing the cluster and any regions particularly impacted by large foreground galaxies or stars. We then re-calculate σ_g as our final global depth for the survey. We adopt the same strategy for the fields that are hampered by large foreground objects, e.g. the Cartwheel and Quintet fields.

We de-lens each of the cluster surveys using the lensing models described in Section 2.1. To measure effective survey areas, we calculate the magnification map taking the mid-point of our redshift bin, z = 11. However, for each candidate, we calculate μ by generating a magnification map using its z_phot. Where there are multiple lensing models considered (i.e. for WHL 0137, MACS 0647, and RXJ 2129), we perform this calculation for each lensing model, and then take the median μ and de-lensed area. For SMACS 0723 and Abell2744 (UNCOVER, GLASS, and DDT 2756) we simply use the magnification and areas calculated using the single lensing model.

For RXJ 2129, we only have one module on the cluster; however, for the three surveys with cluster + parallel survey configurations, the lensing maps typically do not to extend to the parallel module, and so we lose the ability to measure the magnification μ. For cluster + parallel pairs, and for the GLASS and DDT-2756 parallels to Abell 2744, we adopt a constant value of μ = 1.2 where we do not have any lensing map coverage. This value is motivated by the typical μ values found around the edge of the lensing maps of the various surveys. We note that adopting an alternative choice of μ = 1.1 reduces the overall number densities at the 1–2 per cent level, as the parallel regions comprise <15 per cent of the overall survey area.

5.2 Completeness

An accurate determination of the UV LF requires that the incompleteness of the galaxy sample is accounted for. This becomes of particular importance close to the detection limit of the data. Due to our adoption of a more stringent 8σ detection limit throughout, this becomes less of an issue in this survey. The completeness as a function of apparent magnitude in the F200W imaging (F182M for the JEMS data) is determined using the method described in Donnan et al. (2023a). Synthetic galaxies, with a Sérsic profile of n = 1 and half-light radius 0.4 kpc, i.e. comparable to the typical profile and sizes of F115W-dropouts found in Morishita et al. (2023) at similar M₁₅₀₀, were injected into the F200W imaging across each field from m_F200W = 25–30 in steps of 0.05 mag. At each step, 1000-2000 sources (depending on the area of the field) are injected and the rate of recovery measured, producing the completeness as a function of the apparent magnitude in F200W. This process is repeated 50 times for each field. We correct the effective volumes of our candidates by their completeness as determined in this procedure.

In a study of the first Hubble Frontier Fields, Oesch et al. (2015) showed the impact of including shear on completeness simulations, finding that large μ can affect the relative completeness of sources close to the completeness limit. Having removed much of the area around the cluster regions, the effective area of the lensed fields with significant magnification probed in our study presents a modest fraction of the total search area (<10 per cent). This, and the fact that we are conservative in only selecting SNR > 8 sources, allows us to be confident that the impact of shear on our overall completeness and effective survey volumes is negligible.

5.3 The z  = 11 UV LF

To determine the UV LF for our 9.5 ≤ z ≤ 12.5 sample we use the 1/V_max estimator (Schmidt 1968). We only include objects for which our F200W detection is greater significance than 8σ_g, where σ_g is the global depth for a given survey. We split our sample into four number density bins between M₁₅₀₀ = −21.80 and M₁₅₀₀ = −19.55. At the bright-end, we scale the M₁₅₀₀ = −22.57 ground-based LF bin from Donnan et al. (2023a) to match our redshift range, add in the area we searched, and include it as a fifth bin. In Table 4 we present our number densities ϕ in each M₁₅₀₀ bin, along with their errors. To calculate the uncertainty in our number densities, we sum in quadrature the uncertainty found through a bootstrap analysis (σ_b) and cosmic variance (σ_CV) as estimated using the Trenti & Stiavelli (2008) cosmic variance calculator. One of the key advantages of using many non-contiguous fields for this study is that the relative uncertainties across all the bins are significantly reduced.

In the left-hand panel of Fig. 6 we present our determination of the UV LF at z = 11, along with determinations across a host of other studies in the literature. There have been numerous LF measurements based on early CEERS epoch 1 and GLASS data, including Harikane et al. (2023b), Bouwens et al. (2023), and Finkelstein et al. (2023). In the range M₁₅₀₀ = −21.5 to −19.5, we find good agreement with the general consensus in the literature, although we note that they are all over slightly different selection windows: the Bouwens et al. (2023) data points are based on z ≃ 9–11 galaxies, Donnan et al. (2023a) is based on a 9.5 < z < 11.5 sample and Finkelstein et al. (2023) is over 9.5 < z < 12.0. Although we recover all of the candidates in GLASS from Castellano et al. (2023), our number densities at M₁₅₀₀ ≤ −20.5 are significantly lower. This suggests that the GLASS footprint does cover an overdensity, which is smoothed out by having a much larger overall search area in this study. The faintest M₁₅₀₀ bin is in concordance with our earlier work in Donnan et al. (2023a) and indeed with earlier HST z = 10 results from McLeod, McLure & Dunlop (2016). However, the excess in the M₁₅₀₀ < = − 20.5 bins presents a clear departure from HST-based determinations of the UV LF (e.g. McLeod, McLure & Dunlop 2016; Oesch et al. 2018).

In order to increase the overall dynamic range, we incorporate the faintest two bins of the UV LF as determined in Donnan et al. (2023a), encouraged by the excellent agreement between our faintest bin and their bin at a similar M₁₅₀₀. As well as sharing a similar set of reductions, SED fitting codes, and selection criteria, the median redshift of our 9.5 ≤ z ≤ 12.5 sample is z = 10.7, which is close to the mid-point of their bin at z = 10.5. The amalgamation of these two studies allows us to determine a UV LF spanning approximately 4 mag in seven bins, from which we can determine the functional form.

We proceed to fit a double power law, found by Bowler et al. (2014, 2020) and Donnan et al. (2023a) to be the preferred functional form at the highest redshifts, rather than a Schechter (Schechter 1976) function as found in earlier HST studies (e.g. McLeod et al. 2015; McLeod, McLure & Dunlop 2016; Bouwens et al. 2021). We fit the seven number density bins both with a free-fit and with a fixed faint-end slope. The tabulated double power law fits are presented in Table 5, and the combined LF and functional forms are shown in the right-hand panel of Fig. 6.

We find the free-fit double power law is steep at α = −2.59 ± 0.35; however the errors attached to each of the parameters are large, particularly for M^⋆. For the fixed fit, we use α = −2.35, which is the z = 10 slope predicted by the prescription for the evolving double power law as determined by Bowler et al. (2020). It is also the mid-point between the fixed α = −2.10 slope utilized by early JWST studies at similar redshifts such as Donnan et al. (2023a, b) and Harikane et al. (2023b), and the free-fit value determined in this study. The fit with fixed α mitigates the uncertainties of the free fit significantly, and so we adopt this as our fiducial fit for the functional form of the UV LF at z = 11 for the rest of this paper. Encouragingly, our fiducial double power law fit appears to be in good agreement with the faintest bin from Pérez-González et al. (2023), who probe exceptionally deep imaging of HUDF-par2 as part of their MIRI deep survey of the UDF (PID 1283, PI Oestlin).

It should be noted that alternative fits utilizing a range of different fixed faint-end slope are almost indistinguishable over the luminosity range spanned by this study. At z = 8, Donnan et al. (2023a) fitted a faint-end slope α = −2.04 ± 0.29, leveraged by the excellent HST-based constraints of the faint-end which reach M₁₅₀₀ = −17 (McLure et al. 2013). Although our free-fit double power law suggests some steepening between the z = 8 faint-end slope and that of our z = 11 LF, fitting our data with fixed α = −2.04 also presents an acceptable fit over the magnitude range spanned by this study. This does leave the evolution of the faint-end slope at z > 8 an open question, despite all of the recent advances brought with this first tranche of JWST surveys. With deeper data sets such as NGDEEP (ID 2079, PI Finkelstein; Bagley et al. 2023b), studies will be able to constrain α at z > 10 to the accuracy made possible at z = 8 with HST.

As well as drawing comparisons to recent JWST-based determinations of the UV LF, it is instructive to compare to studies probing the bright-end (M₁₅₀₀ < M^⋆) of the UV LF at high-redshift. In the left-hand panel of Fig. 7 we plot our fiducial LF but also include a compendium of studies probing the bright-end of the LF, including a suite of ground-based studies (Stefanon et al. 2019; Bowler et al. 2020; Kauffmann et al. 2022; Donnan et al. 2023a) and combined HST Legacy fields from Bouwens et al. (2023). We also include a comparison of our double power law fit to those determined by Donnan et al. (2023a) in the right-hand panel.

It is clear from Fig. 7 that our results suggest there is very little evolution in the bright end of the UV LF between z = 9–11, a feature previously reported by Bowler et al. (2020) between z = 8–10. In the left-hand panel, it can be seen that the z = 9–10 literature points at M₁₅₀₀ < −21.50 scatter symmetrically around our z = 11 double power law fit. Moreover, in the right-hand panel, it can be seen that at M₁₅₀₀ < −20.50 there is little to distinguish our z = 11 fit from the z = 8–10.5 fits derived by Donnan et al. (2023a).

A relatively un-evolving bright-end has also been reported by Finkelstein et al. (2022a) over the CANDELS fields (see also Finkelstein & Bagley 2022). Although there have been some surprisingly high number densities of bright z > 10 galaxies reported in early JWST studies (Castellano et al. 2022, 2023), reports of high number densities of bright z > 8 galaxies had become commonplace even pre-JWST, with studies reporting high number densities with wide-area pure parallel HST studies (Rojas-Ruiz et al. 2020; Bagley et al. 2022), albeit with a greater potential risk of contamination due to limited filter sets.

The interpretation for greater numbers of bright z > 8 galaxies has been discussed previously in Bowler et al. (2020), who suggest that the double power law shape reflects a reduction in dust attenuation and/or inefficient mass quenching. Another explanation could lie in the star-formation efficiency (SFE), which may be higher at the bright-end to produce the double power law shape at high-redshift (e.g. Harikane et al. 2022). After the early JWST studies reporting large number densities of bright z > 8 galaxies, Mason, Trenti & Treu (2023) presented predictions of z > 8 LFs under the assumption of 100 per cent SFE. They demonstrated that a 100 per cent SFE scenario yields LFs several orders of magnitude higher than current results, i.e. observations thus far are fully consistent with current cosmological frameworks. They conclude that our view may be biased in that we are predominantly witnessing young galaxies that are undergoing bursts of star-formation, rather than a representative population.

There is also the potential for significant field-to-field variation resulting in skewed number densities at the bright-end, where number statistics are generally low. For example, (Finkelstein et al. 2022a) found the CANDELS EGS field to be over-dense, with 7/11 of their bright 8.5 < z < 10.5 candidates across the five CANDELS fields coming from that field alone. Two of these candidates were spectroscopically confirmed with MOSFIRE observations to lie at z = 8.7, EGSY8p7 by Zitrin et al. (2015b) and EGS_z910_44164 by Larson et al. (2022). The EGS overdensity at z = 8.7 could give rise to ionising bubbles (Larson et al. 2022). Further evidence for an ionising bubble was provided by the discovery of an overdensity of z ≃ 9 candidates found near EGSY8p7 using CEERS JWST NIRCam imaging (Whitler et al. 2023). Also recently, the over-density found in GLASS by (Castellano et al. 2023) causes their number densities to lie ∼3 × higher than other JWST studies.

This study mitigates some of the field-to-field variation with numerous non-contiguous pointings over a relatively large area, while also bridging some of the disconnect between the ground-based and early JWST results at z > 9.5. Yet more progress at the bright-end of the UV LF can be gained through using wide-area surveys such as PRIMER (ID 1837; PI Dunlop) and COSMOS-Web (ID 1727, PI Casey; Casey et al. 2023).

Before moving on to the star-formation rate density, we briefly discuss a higher redshift bin generated using our 12.5 < z < 14.5 candidates. The number densities are tabulated in Table 4. As we currently have no constraints on the number densities of 12.5 < z < 14.5 galaxies at M₁₅₀₀ ≤ −20, we cannot yet comment on the evolution of the bright-end of the UV LF at these redshifts. However, it is nevertheless instructive to discern how the the number densities evolve between our redshift bins. We hence fit a double power law to our z ∼ 13.5 LF, fixing α, β, and M^⋆ to their fiducial z = 11 values. This is illustrated in Fig. 8. By allowing ϕ^⋆ to float, we find a factor of four decrease in ϕ^⋆ with respect to our fiducial z = 11 LF, i.e. ϕ^⋆ = (5.14 ± 2.40) × 10⁻⁶. Future wide-area searches for bright z > 12.5 galaxies will be required in order to further determine if the lack of evolution in the bright-end of the UV LF extends to such early times.

6 STAR-FORMATION RATE DENSITY

We proceed to calculate the luminosity-weighted integral for our functional fits to the UV LF in order to determine the UV luminosity density, ρ_UV. The integral limit that we choose is M₁₅₀₀ = −17, as has been convention with previous HST-based studies. Our ρ_UV values are tabulated along with the UV LF functions in Table 5. To convert this to a star-formation rate density, ρ_SFR, we use the conversion factor |$\cal {K}_{\rm UV} = 1.15 \times 10^{-28}$| M_⊙ yr⁻¹/erg s⁻¹ Hz⁻¹ (Madau & Dickinson 2014). The free-fit and fixed α = −2.35 double power law cases yield very similar values of ρ_UV. We adopt |$\mathrm{\log }_{10}{(\rho _{\mathrm{UV}})} = 25.15^{+0.13}_{-0.14}$| as our fiducial value.

We also include a tentative measurement of ρ_UV at z = 13.5, based on our 12.5 < z < 14.5 sample. Integrating our z = 13.5 LF fit yields |$\mathrm{\log }_{10}{(\rho _{\mathrm{UV}})} = 24.54^{+0.16}_{-0.25}$| at z = 13.5.

In the left-hand panel of Fig. 9 we present our determination of ρ_UV along with those from other observations in the literature. Overplotted is the previous HST-based results from Oesch et al. (2018) and their rapidly declining ρ_SFR(z) relation, following the evolution of the dark matter halo mass function. We also include the constant SFE model from Harikane et al. (2022), which follows ρ_SFR ∝ (1 + z)^−0.5, and the log-linear relation we previously determined in Donnan et al. (2023a),

As can be seen from the left-hand panel, our new determination of ρ_UV, and hence ρ_SFR, at z = 11 is found to be slightly higher than, although still consistent with, our previous estimates from Donnan et al. (2023a, b). With our significantly larger area spread over numerous additional fields, we once again find evidence in support of a smoothly declining ρ_SFR at least to z = 12. Although tentative, our z = 13.5 data point also appears to be consistent with the results from Donnan et al. (2023a), Harikane et al. (2023b), and Pérez-González et al. (2023), and suggests the smooth decline in ρ_SFR may yet continue beyond z = 12.

Given that the ρ_UV can be sensitive to the assumed faint-end slope, we want to ensure that our result is robust against the assumed α slope, particularly as the current JWST data are still unable to robustly constrain this parameter to the same level seen at z = 8 with HST. Hence, we also explore the impact of adopting α = −2.04, which was the slope found at z = 8 by Donnan et al. (2023a) and still presents an acceptable fit to our UV LF at z = 11. We find that |$\mathrm{\log }_{10}{(\rho _{\mathrm{UV}})}(\mathrm{z} = 11) = 25.09^{+0.14}_{-0.14}$|⁠, which is in good agreement with the log-linear relation from Donnan et al. (2023a). Our conclusion of a smoothly declining ρ_SFR(z) between z = 8–11 is therefore robust against the assumed value of α.

Pre-JWST, the rate that ρ_SFR declined at z > 8 was a matter of debate, with Oesch et al. (2018) advocating a rapid decline based on a lack of z = 10 galaxies found among HST Legacy fields. This was a conclusion shared by Ishigaki et al. (2018) when integrating to a limit of M₁₅₀₀ = −17, the limit of HST’s capabilities, although they noted that this depends on the choice of M₁₅₀₀ limit, and that a limit of M₁₅₀₀ = −15 resulted in a smoothly declining ρ_SFR. Many theoretical models of galaxy evolution also follow such a trend of rapidly declining ρ_SFR(z) beyond z > 8. In the right-hand panel of Fig. 9, we include a comparison between our results and those of theoretical models of galaxy evolution, where ρ_UV has been determined by integrating to the same M₁₅₀₀ limit. Our new results lie above a number of theoretical models that follow more closely the rapidly declining relation, such as Yung et al. (2019), FiBy (Paardekooper, Khochfar & Dalla 2013), Millenium TNG (Kannan et al. 2023), and Mason, Trenti & Treu (2015). However, the Behroozi & Silk (2015), Ferrara, Pallottini & Dayal (2023), and FLARES (Wilkins et al. 2023) models appear to be entirely consistent with our z = 11 results.

As well as lying above the Oesch et al. (2018) rapidly declining ρ_SFR relation we also lie above the constant SFE model from Harikane et al. (2022). Although the constant SFE model is in close agreement with Harikane et al. (2023b) and Donnan et al. (2023a) up to z = 9, the observations begin to diverge from the model at higher redshifts. This is discussed in detail by Harikane et al. (2023b), who show that increasing the SFE at z > 10, or adjusting the IMF to be more top-heavy, can account for such discrepancies. Alternatively, as was mentioned earlier, Mason, Trenti & Treu (2023) suggest that observations thus far with JWST may be biased towards young galaxies with copious star-formation, drawn from the upper envelope of the |$M_{1500}-\mathcal {M}_{\mathrm{halo}}$| relation, and that a more representative sample may be found when probing fainter apparent magnitudes, e.g. m₂₀₀ > 30.

On the other hand, a study by Bouwens et al. (2023) suggests that perhaps the star-formation rate density may yet lie higher than the results of Harikane et al. (2023b), Donnan et al. (2023a), and this study. They compute ρ_SFR when considering a more inclusive (but still referred to as ‘solid’) sample of high-redshift candidates, rather than solely the most ‘robust’ candidates, and find a significantly higher ρ_SFR. Their study highlights the requirement for spectroscopy to confirm or refute the high-redshift nature of candidates. That said, it is worth noting that early JWST NIRSpec results from e.g. Curtis-Lake et al. (2023), Roberts-Borsani et al. (2023), Bunker et al. (2023), and Arrabal Haro et al. (2023a) suggest that z = 9–11 candidate selection with photometric redshifts has been largely successful.

7 CONCLUSIONS

We have presented a search for high signal-to-noise z ≥ 9.5 galaxy candidates across a suite of 12 JWST surveys spanning a raw survey area of ∼260 arcmin². We uncover 61 z ≥ 9.5 candidates detected at 8σ, of which 18 are at z ≥ 11.5. The exceptional brightness of many of these candidates makes them ideal for follow-up spectroscopy with NIRSpec.

With the inclusion of six additional candidates found in our earlier study, Donnan et al. (2023b), we proceed to construct the UV LF over 9.5 < z < 12.5. When combined with faint-end constraints from Donnan et al. (2023a), and UVISTA-based constraints at the bright-end, our z = 11 LF spans 4 mag of dynamic range in M₁₅₀₀. Our LF is consistent with previous JWST-based number densities in the literature at similar redshifts, but with tighter constraints owing to our wider search area. At the bright-end, our LF is consistent with a lack of evolution between z = 9–11, similar to results from Bowler et al. (2020), who found such trends between z = 8–10. This lack of evolution suggests that we are entering a regime where mass quenching efficiency is lower or that there is reduced dust attenuation (Bowler et al. 2020).

At the faint-end, our α measurements suggest a steep α, which we fix to α = −2.35 for our fiducial LF. With the present data, we are not yet able to robustly determine how α evolves between z = 8–11, given that a fit with α = −2.04 (as measured at z = 8 by Donnan et al. 2023a) still provides an acceptable fit.

By integrating our LF, we arrive at the luminosity density and hence star-formation rate density. Our ρ_UV at z = 11 lies slightly above, but is still consistent with, our previous log-linear ρ_UV(z) relationship from Donnan et al. (2023a). Our ρ_SFR at z = 11 suggests a continued smoothly declining ρ_SFR to at least z = 12, with tentative z = 13.5 measurements based on a handful of z > 12.5 candidates further reinforcing this. This result lies above both the rapidly declining ρ_SFR suggested in pre-JWST studies (e.g. Ishigaki et al. 2018; Oesch et al. 2018) and the constant SFE model from Harikane et al. (2022), as well as numerous theoretical models of galaxy evolution.

There are numerous potential interpretations for the high number densities of bright z ≃ 10 galaxies and the higher ρ_SFR than predicted by models. At the highest redshifts, we may be seeing a rise in the SFE (Harikane et al. 2023a), or that we are potentially viewing the upper envelope of the |$M_{1500}-\mathcal {M}_{\mathrm{halo}}$| relation (Mason, Trenti & Treu 2023). Although there have been suggestions that the early JWST studies may be subject to significant contamination by low-redshift interlopers, particularly as the overlap in sources between studies has hitherto been relatively low [see discussion by Bouwens et al. (2023)], we can be encouraged by the number of spectroscopic confirmations of many photometrically selected z > 8 candidates (e.g. Curtis-Lake et al. 2023; Arrabal Haro et al. 2023a, b).

The prospects for building upon this LF determination with further JWST Cycle 1 surveys are excellent. In the bright-to-intermediate regime, PRIMER will be instrumental in providing homogeneous coverage of COSMOS and UDS at even greater areas than those probed in this study, and to comparable depths. At the very brightest luminosities, COSMOS-Web (Casey et al. 2023) will provide >0.5 deg² of imaging, which will be crucial for uncovering the brightest (and rarest) high-redshift galaxies. At the faint-end, exceptionally deep surveys such as NGDEEP (Bagley et al. 2023b) and JADES (Robertson et al. 2023) will help to better constrain the evolution of the faint-end slope α.

ACKNOWLEDGEMENTS

We would like thank the referee for very useful comments and a constructive report, which have helped improve the manuscript. We would like to express our appreciation to all of those who have worked tremendously hard to make all of these JWST surveys happen. We are also thankful that all of these data sets were made publicly available immediately.

We would also like to thank all those involved in the creation of the ancillary HST data that was publicly available, and to all of those who have provided gravitational lensing models publicly for the community.

We would like to thank Adi Zitrin for feedback on the lensing models utilized, and Marco Lam for useful discussions. We would also like to thank Andrea Ferrara for useful discussion, and for providing additional data from the Ferrara, Pallottini & Dayal (2023) model.

This work is based in part on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with programs 1063, 1324, 1345, 1433, 1963, 2282, 2561, 2727, 2732, 2736, 2738, 2756 and 2767. The authors acknowledge the associated teams for developing their observing programs with a zero-exclusive-access period.

This work is based in part on observations obtained with the NASA/ESA Hubble Space Telescope, retrieved from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute (STScI). STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555.

This work utilizes gravitational lensing models produced by PIs Bradac, Natarajan & Kneib (CATS), Merten & Zitrin, Sharon, Williams, Keeton, Bernstein and Diego, and the glafic group. This lens modeling was partially funded by the HST Frontier Fields program conducted by STScI. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. The lens models were obtained from the Mikulski Archive for Space Telescopes (MAST).

This research has made use of the SVO Filter Profile Service (http://svo2.cab.inta-csic.es/theory/fps/) supported from the Spanish MINECO through grant AYA2017-84089 (Rodrigo, Solano & Bayo 2012; Rodrigo & Solano 2020)

This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency.

This work is based on observations taken by the RELICS Treasury Program (GO 14096) with the NASA/ESA HST, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.

This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

DJM, CD, RJM, JSD, RB, and MLH acknowledge the support of the Science and Technology Facilities Council. ACC acknowledges the Leverhulme Trust for their support via a Leverhulme Early Career Fellowship. FC and TMS acknowledge support from a UKRI Frontier Research Guarantee Grant (grant reference EP/X021025/1). RSE acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 669253).

For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

DATA AVAILABILITY

At the time of writing, the data sets used in this manuscript were all publicly available from MAST and on the CADC.

References

APPENDIX A: DISCUSSION OF HIGH-REDSHIFT CANDIDATES

In the following sections, we discuss the results of our search for high-redshift galaxies field-by-field, along with a (non-exhaustive) comparison to some previous searches in the literature. The variation between JWST high-redshift candidate samples has been a theme in the literature [see Bouwens et al. (2023) for a comprehensive discussion]. We stress that differences in the candidate lists can result from a number of factors such as differing reductions, photometry measurement and selection criteria such as non-detection and Δχ² thresholds.

The nomenclature adopted for our final candidate IDs is Field-redshift-catalogueID (e.g. GLASS-z11-1481). When referring to candidates from other studies we use their naming convention.

A1 SMACS 0723

We select one high-redshift candidate from SMACS 0723 in our final sample, SMACS0723-z12-7442.

There have been numerous studies of this cluster field to date, yielding contrasting samples. We only compare candidates that have passed the additional criteria of F150W − F200W > 0.75 or F150W < 2σ. Adams et al. (2023b) reported a sample of four 9 < z < 12 galaxies, but only one of which can be compared to our final sample. Their candidate with ID 10234 was initially selected in our sample, where we found |$z_{\mathrm{phot}} = 11.8^{+0.2}_{-0.1}$|⁠. However, there appeared to be a marginal F090W detection upon visual inspection of the imaging, and so we discarded the candidate.

Atek et al. (2023) reported a sample of high-redshift candidates, including two potential z > 15 candidates and two potential z = 12 candidates. We do not include their ID z16a due to a >3σ detection in F150W, which rules out a potential z > 13 solution. While we do find a |$z_{\mathrm{phot}} = 15.6^{+0.5}_{-0.9}$| solution for their ID z16b, we exclude it from our final sample as it has Δχ² ≃ 2. Both their objects z12a and z12b have Δχ² ≃ 1 between the high-redshift and low-redshift solutions, so are also excluded.

Most of the SMACS 0723 sample from Donnan et al. (2023a) was not included in our final sample due to our more stringent selection and brightness criteria. Their candidate with ID 8347 is in our final sample as SMACS0723-z12-7442. While we find |$z_{\mathrm{phot}} = 12.2^{+0.5}_{-0.4}$| for their object with ID 1566, it just misses out on being in our final sample through having Δχ² = 3.4. Finally, their candidate with ID 10566 had been excluded from our final sample due to a >2σ detection in the HST ACS imaging in our present analysis, potentially due to a noise spike that contaminated the aperture.

A2 GLASS

Using the maximum-depth combination of the epoch1 + 2 GLASS survey, we report a sample of 12 z ≥ 9.3 galaxy candidates, including five in the range 10.5 < z < 11.5 and three at z > 12. The high density of bright z ∼ 10 galaxies seen in Castellano et al. (2023) is confirmed by this study, as we re-select all five of their candidates at similar redshifts.

With the first epoch of GLASS imaging, Naidu et al. (2022) and Castellano et al. (2022) unveiled two particularly bright z = 10 − 12 galaxies. We recover both of these objects in our final sample as GLASS-z10-3283 (GLASS-z10/GHZ1) and GLASS-z12-28072 (GLASS-z12/GHZ2).

A comprehensive study by Bouwens et al. (2023) also covered the GLASS data set. Their final sample includes the two aforementioned candidates in common with our study, Naidu et al. (2022) and Castellano et al. (2022). We re-select their z = 13.7 candidate, GLASSP2H-3576218534 in our final sample as GLASS-z14-33570. Although we also recover a high-redshift solution for their candidate ID 4002721259, the Δχ² was insufficient to include in our final sample. Their alternative sample includes an additional z = 13.3 candidate, GLASSP1H-4015021230, but we cannot confirm its redshift owing to a flat p(z) distribution. Finally, GLASSP1YJ-4016420474 and GLASSP1YJ-4043920397 are z = 9.6 and z = 9.7 candidates that we do not include as they appear to have a marginal detection in F115W. This is likely a result of the increased depth when adding epoch 2, as the epoch 1-only version of the data appeared consistent with non-detections for both objects.

A3 CEERS

The first epoch of observations was split into two fields, NE and SW. In CEERS-NE we selected six candidates, among which are CEERS-93316 (Donnan et al. 2023a) and Maisie’s galaxy (Finkelstein et al. 2022b). While the former has recently been revealed by spectroscopy to be lower redshift, Maisie’s galaxy has been spectroscopically confirmed to be at high-redshift (DDT 2750, PI Arrabal Haro; Arrabal Haro et al. 2023a). We also note that CEERS-NE-z11-2543 has also been recently spectroscopically confirmed at high-redshift with the same program. Given its new spectroscopic redshift z_spec = 4.9, we exclude CEERS-93316 from the subsequent analysis. Four of the five remaining candidates were included in Finkelstein et al. (2023).

There are five galaxies from Finkelstein et al. (2023) that we also find to be z > 9.0, but do not include in our final sample due to insufficient Δχ² and/or <8σ detection. A further six candidates are not included due to >2σ detection in F115W, but five of those are consistent with their z_phot = 8–9 solutions. The only higher redshift candidate excluded by our study for a F115W detection is their candidate ID CEERS-7227, for which we find the F115W to be marginally detected at 2.6σ. We note that Finkelstein et al. (2023) find a 2.2σ flux for this object in F115W. Only CEERS-NE-z10-3069854 is unique to this study.

In CEERS-SW, we find six more high-redshift candidates, although only two of the final sample are in common with Finkelstein et al. (2023; our IDs CEERS-SW-z11-26592 and CEERS-SW-z11-3961). Our candidates CEERS-SW-z11-20687 and CEERS-SW-z14-37967 were also found in Bouwens et al. (2023) and Donnan et al. (2023a) respectively.

Four objects in the Finkelstein et al. (2023) sample (their IDs 7603, 1748, 2324, and 4012) were also found to have preferred high-redshift solutions in our catalogues, but were detected at <8σ so were excluded. Finally, four z = 8–9 candidates were F115W>2σ and so were excluded.

Epoch two of the CEERS observations were split into six NIRCam pointing pairs, and we adopt this naming convention into our candidate IDs. We report 12 high-redshift candidates across the ∼60 arcmin². Two of these objects (CEERS-2-5-z10-9617 and CEERS-2-6-z10-3530) have very recently been independently confirmed with spectroscopy to lie at z ≃ 10 (Arrabal Haro et al. 2023b).

A4 Stephan’s Quintet

Our search over the Quintet field yielded one high-redshift candidate over 11.5 < z < 12.5. Harikane et al. (2023b) also reported two candidates at z > 12 from the Quintet field. Their z ∼ 12 candidate was in our initial catalogue, and with our photometry appears to be a robust candidate at |$z_{\mathrm{phot}} = 11.7^{+0.4}_{-0.5}$| and Δχ² = 4.8. It also passes our additional F150W-F200W colour criteria, but we do not include it in our final sample as it is only detected at ∼6σ in F200W and F277W. For the z = 16 candidate, we obtain a low-z preferred solution of z = 0.6, with other solutions at z = 16.4 and z = 4.9 being of similar probability (within 1σ). We cannot rule out that it is indeed a z = 16.4 galaxy, but based on the current fits we do not include it in our final sample. Interestingly, the high- and low-redshift solutions that we derive are very similar to CEERS-93316. The low-redshift solution with emission lines was ultimately the spectroscopic redshift, suggesting that candidates around z = 16.4 are vulnerable to this class of contaminant, especially in the absence of photometric constraints from key medium bands such as F250M and F300M.

A5 UNCOVER

From our analysis of the UNCOVER field, we amass a sample of sixteen candidates. We find UHZ1 from Castellano et al. (2023) as UNCOVER-z10-32757 with a similar |$z_{\mathrm{phot}} = 10.4^{+0.1}_{-0.2}$| and Δχ² > 100. One additional candidate was a duplicate of GLASS-z10-3283 and so was excluded.

We also find two components from the triply lensed z ≃ 10 galaxy from Zitrin et al. (2014), which has now been spectroscopically confirmed to lie at z = 9.8 (Roberts-Borsani et al. 2023). We find JD1B as UNCOVER-z11-15452 which was observed to be exceptionally bright M₁₅₀₀ = −20.9, but is magnified ≃ 11 ×. JD1A was in our initial catalogue as UNCOVER-z10-15137 with a robust |$z_{\mathrm{phot}} = 10.3^{+0.3}_{-0.9}$|⁠, Δχ² = 17 solution. This component had a magnification of 13.4 and M₁₅₀₀ = −17.3. As it was only 7σ_local in each of our detection catalogues, it was not propagated to the final sample. The remaining component, JD1C, was of insufficient signal-to-noise to make it into our catalogues at >5σ, due to enhanced background and crowding.

A6 DDT 2756

In the DDT 2756 parallel field to Abell2744, we found more examples of bright and robust z ∼ 10 candidates, providing more evidence that the region around Abell2744 is over-dense with such galaxies. We report five 9.5 < z < 11.7 candidates, including a particularly robust (Δχ² ≃ 43) candidate with |$z_{\mathrm{phot}} = 10.6^{+0.1}_{-0.3}$|⁠, DDT2756-z11-1001979.

The DDT 2756 survey was also explored in Castellano et al. (2023), who found one remarkably bright z ∼ 10 candidate, DHZ1. Due to the high-resolution of the F200W imaging, our F200W-detected catalogue had de-blended DHZ1 into two components. The NE component was selected for our final sample as DDT2756-z10-1010612, and yielded a robust high-redshift solution of |$z_{\mathrm{phot}} = 9.5^{+0.6}_{-0.3}$|⁠. Although the SW component was excluded from our final sample due to a 2.2σ detection in F115W, we performed SED fitting and found an ultra robust z_phot = 9.5 ± 0.1 (Δχ² > 150) solution, consistent with the NE component. The system has since been spectroscopically confirmed to be z = 9.31 (Boyett et al. 2023).

A7 WHL 0137

We report four high-redshift candidates in the WHL 0137 cluster + parallel data set, three of which are z ≃ 9 candidates and the other with z = 11.1. There has been one reported search for z > 9 galaxies over the WHL 0137 data set to date, with Bradley et al. (2023) uncovering a sample of four galaxies over 8.3 < z < 10.2. We include in our final catalogue one of the four candidates, WHL0137-z11-22312 (their WHL0137-09319). While we found a similar high-redshift solution for their candidate WHL0137-10060, the Δχ² was insufficient to include this candidate in our final sample. Finally, their WHL0137-08004 and WHL0137-1332 candidates were excluded from our sample due to detections in F115W.

A8 MACS 0647

We report three high-redshift candidates in the MACS 0647 cluster + parallel data set. There has been a detailed analysis of MACS 0647 JD (Coe et al. 2013) by Hsiao et al. (2023). We recover two of the three multiply imaged components of JD: MACS0647-z11-20148 is JD1 and MACS0647-z11-26400 is JD3. Both have extremely robust solutions Δχ² > 100, and in excellent agreement with the reported value of z_phot = 10.6 ± 0.3 in Hsiao et al. (2023). This object has since been spectroscopically confirmed to lie at z = 10.17 by Harikane et al. (2023a).

We do not include JD2 in our final sample as it has a 2.8σ detection in F435W, which aligns with Hsiao et al. (2023)’s F435W photometry also being >2σ. For completeness, we note that we recover a z = 10.8 ± 0.2 solution for JD2 using our photometry, with a Δχ² = 63 showing it is an extremely robust object despite the F435W SNR>2 measurement. For our subsequent LF analysis in Section 5, we only propagate JD1 so as not to double-count multiple images. In addition to MACS0647 JD, we find a notably robust z = 9.5 galaxy in the cluster module, MACS0647-z9-20158.

A9 RXJ 2129

We report two z ≃ 9 galaxy candidates in this field, however we do not uncover any z ≥ 9.5 candidates. To date, there has been one reported z ≥ 9.5 galaxy in the RXJ 2129 field with the new JWST data (Williams et al. 2023b). This galaxy is reported to be intrinsically very faint with M₁₅₀₀ = −17.4, but high magnification factor (μ ≃ 20). We do not recover this spectroscopically confirmed galaxy due to its position near the RXJ 2129 cluster. The enhanced background in the region around the object had caused the object to fail our signal-to-noise cuts in all of the catalogues.

A10 NEP TDF

Our study of the NEP TDF data set uncovered one z ∼ 10 candidate and one z ∼ 11 candidate. To date, there have been no published candidates in the literature. Adams et al. (2023a) have recently presented 7.5 ≤ z ≤ 13.5 LFs based on a search including candidates across this field, but the associated candidate list is yet to be published (Conselice et al., in preparation).

A11 Cartwheel

Our final sample does not include any galaxy candidates at z > 9 from the Cartwheel galaxy field.

A12 J1235

From our study of the J1235 field we report three high-redshift galaxy candidates over 9.7 ≤ z ≤ 11.8. To date, there have been no other published high-redshift candidates in the literature.