Two Remarkably Luminous Galaxy Candidates at z ≈ 10–12 Revealed by JWST

Our analysis is based on some of the first JWST/NIRCam data sets that have been observed and released over extragalactic fields. In particular, we analyze the two Early Release Science (ERS) programs: GLASS and CEERS.

The first is the NIRCam parallel field of the ERS program GLASS (PID: 1324; Treu et al. 2022). This program obtained a single NIRCam pointing in seven wide filters F090W, F115W, F150W, F200W, F277W, F356W, and F444W, observed for 3.3, 3.3, 1.7, 1.5, 1.5, 1.7, and 6.6 hr, respectively.

We also analyzed the first four NIRCam pointings from the ERS program CEERS (PID: 1345; S. Finkelstein et al. 2022, in preparation) that have been observed to date. The CEERS images include F115W+F277W, F115W+F356W, F150W+F410M, and F200W+F444W short- (SW) and long-wavelength (LW) exposure sets, for a typical integration time of 0.8 hr per filter, except for F115W that obtained double this integration time.

The combined area of these five NIRCam fields used in our analysis amounts to 49 arcmin² (≈10 arcmin² in GLASS, and ≈40 arcmin² in CEERS) reaching an unprecedented 5σ depth that ranges between 28.6 and 29.6 AB mag at 4 μm as measured in 032 diameter apertures (see Table 1).

The NIRCam images were reduced with the standard JWST pipeline up to stage two using the reference files from jwst_0942.pmap, as well as our own sky flats. Additional, chip-dependent zero-point offsets were applied based on our own reductions of the flux standard star J1743045 as well as observations of the Large Magellanic Cloud. For some filters, the individual corrections in different modules of the NIRCam detector can reach up to 30%. ²⁶ The resulting images were then processed with the grizli ²⁷ pipeline for proper alignment onto a common world coordinate system that was matched to the Gaia DR3 catalogs. To do this, grizli recomputes the traditional simple imaging polynomial distortion headers that were common for HST data. This allows us to use drizzlepac as with HST images to combine the individual frames and produce distortion-corrected mosaics. Additionally, grizli mitigates the 1/f noise and masks "snowballs" that are most prominent in short-wavelength filters.

The pipeline was used to drizzle images at 20 mas/40 mas pixels for the SW and LW data, respectively. In the following, our analysis is based on pixel-matched 40 mas images, however. For more information on the grizli processing see G. Brammer et al. (2022, in preparation).

After processing the new NIRCam images, we produced photometric catalogs in all fields using the SExtractor software (Bertin & Arnouts 1996). Sources were detected in dual mode with two different detection images: F200W or a weighted combination of all the LW filters. Fluxes were measured in small circular apertures of 032 diameter and were corrected to total using the AUTO flux measurement from the detection image. An additional correction of typically a few percent only was applied for remaining flux outside of this aperture based on the predicted encircled energy for the JWST point-spread functions (PSFs). Flux uncertainties were estimated based on sigma-clipped histograms of circular apertures placed throughout the images in random sky positions. These were then used to rescale the drizzled rms maps. Thus, our uncertainties are as close to the data as possible. The 5σ depths derived in this way are listed in Table 1.

Given that JWST is a completely new facility for which calibration is still ongoing, it is important to test the resulting images for any issues. Indeed, the NIRCam data revealed several features that are not accounted for in the standard pipeline. In particular, the SW data suffer from significant scattered light, if there are bright stars in the vicinity. This is particularly pronounced in the GLASS parallel field where a 10th mag star just outside of the field seems to cause artificial images across the field. This is particularly pronounced in the F090W filter in the B4 detector. However, also other detectors and filters seem to be affected, albeit to a lower extent. While this issue could be overcome with improved processing, or with additional observations at a different roll angle, this does not significantly affect the current analysis. Since we are searching for very high-redshift galaxies that disappear at shorter wavelengths, these data issues only introduce some level of incompleteness. Most importantly, the areas of the two candidate sources presented later in this Paper are not affected.

Where available, we compared our JWST photometry in the shorter wavelength filters with existing HST data to check both for issues with residual distortion or with magnitude zero-point offsets. In particular, for the GLASS parallel field, we used HST images made available by Kokorev et al. (2022) from the BUFFALO survey (Steinhardt et al. 2020). However, only a small portion of the field is covered by both HST and JWST. For the CEERS data, we make use of rereductions from the imaging taken by the CANDELS survey (e.g., Koekemoer et al. 2013). No significant offsets or issues were detected. Zero-point corrections remained small (<15%; see the next section).

We confirm the photometric redshifts for the two GLASS candidates and derive stellar population properties using the Prospector SED fitting code (Leja et al. 2017, 2019; Johnson et al. 2021). The SED parameter space explored by Prospector is more expansive than EAZY's linear template combinations, and therefore it acts as an important check on our derived redshifts. We use FSPS (Conroy et al. 2009, 2010; Conroy & Gunn 2010a) with the MIST stellar models (Choi et al. 2017). We adopt the 19 parameter physical model and parameter choices described in Tacchella et al. (2022) that fits for the redshift, stellar and gas-phase metallicities, stellar mass, star formation history, dust properties, active galactic nuclei (AGN) emission, and scaling of the intergalactic Medium (IGM) attenuation curve. We make slight modifications to their setup—in particular we explore a broader redshift range of z = 0.1–20 and keep two bins fixed at lookback times of 0–5 and 5–10 Myr in the star formation history following Whitler et al. (2022) to capture recent bursts that may be powering extreme nebular emission expected to occur generically at the redshifts of interest (e.g., Labbe et al. 2013; De Barros et al. 2019; Endsley et al. 2020). We adopt a "continuity" prior on the star formation history, which limits the amount of variance across consecutive time bins resulting in smooth histories (Leja et al. 2019; Tacchella et al. 2022). For further details we direct readers to Table 1 and Section 3.4 of Tacchella et al. (2022).

The redshift fits from Prospector are in excellent agreement with EAZY—we find ${z}_{\mathrm{Prospector}}={10.4}_{-0.5}^{+0.4}$ for GL-z10 and ${z}_{\mathrm{Prospector}}={12.4}_{-0.3}^{+0.1}$ for GL-z12. The photometry and redshift inference for these sources are summarized in Figures 1 and 2, with fluxes listed in Table 2. We confirm that no significant data quality issues affect the z > 10 candidacy of the sources in the imaging. We derive p(z > 10) ≈ 100% for GL-z12 and p(z > 9.4) ≈ 100% for GL-z10, with their dramatic ≳1.8 mag breaks explained by total absorption of photons bluewards of Lyα by neutral hydrogen in the intergalactic medium.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Both galaxies are detected at very high significance in all filters longward of their break, by virtue of our selection. While they appear very luminous in the JWST data, these sources have UV absolute magnitudes (M_UV ≈ −21) that correspond to ${L}_{\mathrm{UV}}^{* }$ at z ∼ 8–10 (see, e.g., Bouwens et al. 2021). This also makes them 1 mag fainter than GN-z11 and even 2.5 mag fainter than the possible z ∼ 13 galaxy candidate HD1 (Harikane et al. 2022). Hence, these sources are not really extreme outliers (see also Figure 3). Nevertheless, it is interesting that the first few images with JWST already reveal two such bright sources. We will discuss their implications on the UV luminosity function (LF) in a later section.

Zoom In Zoom Out Reset image size

The nondetections of both sources in deep, shorter wavelength images essentially rule out a lower redshift solution. Nevertheless, it is interesting to explore the nature of possible contaminants. We thus rerun our photometric redshift codes and force them to find lower redshift fits. The best z < 6 solutions in our low-z EAZY runs for these sources are ≈10⁸⁻⁹ M_⊙ quiescent galaxies at z ≈ 3.4 (z ≈ 2.5) with Balmer breaks straddling the dropout filter (silver SEDs in Figures 1 and 2). However note that Balmer breaks, even in the most pathological cases (e.g., 4000 Å falls just redward of the dropout filter in a supersolar metallicity galaxy as old as the age of the universe at z ≈ 2–3.5), can only produce drops of ≲1.5 mag (assuming no attenuation). The best-fit low-z solutions have A_V ≈ 0.1—stronger attenuation that deepens the break is disfavored by the blue continuum slope at wavelengths longer than the break. In other words, the best-fit low-z solutions predict >5σ detections in bands where we find no flux, and continuum slopes redder than we observe.

In order to allow for possible systematic effects in the new JWST data, we perform further testing. We refit redshifts to multiple versions of photometry for these sources—e.g., by adding PSF corrections using WebbPSF, by increasing the error floor on the photometry, and by extracting photometry using different apertures and detection bands. The only test that produces viable low-z solutions is when we set a 10 nJy error floor on all photometry—this is roughly the level in the SW filters at which the strongest Balmer breaks at z ≈ 2–3 can no longer be ruled out (see open silver squares in bottom-left panels of Figures 1 and 2). This test is a vivid demonstration of why the sensitivity of JWST/NIRCam is required to identify objects like GL-z10 and GL-z12 with confidence.

While the discovery of GN-z11 has already demonstrated that the formation of a billion solar mass galaxies was well underway at ∼400 Myr after the Big Bang, the discovery of these two new sources allows us to derive further constraints on the physical properties of galaxies at this very early epoch of the universe. The Prospector results are summarized in Table 3. In order to efficiently sample the redshift range of interest, we assume a tighter redshift prior (a Gaussian centered on the EAZY p(z) with the width set to the 84th–16th percentile) than in our previous runs when fitting for the redshift.

Table 3. Summary of Properties

	GL-z10	GL-z12
R.A.	+0:14:02.86	+0:13:59.76
Decl.	−30:22:18.7	−30:19:29.1
Redshift zProspector	${10.4}_{-0.5}^{+0.4}$	${12.4}_{-0.3}^{+0.1}$
Redshift zEAZY	${10.4}_{-0.7}^{+0.2}$	${12.2}_{-0.2}^{+0.1}$
Stellar Mass log(M⋆/M⊙)	${9.6}_{-0.4}^{+0.2}$	${9.1}_{-0.4}^{+0.3}$
UV Luminosity (MUV)	$-{20.7}_{-0.2}^{+0.2}$	$-{21.0}_{-0.1}^{+0.1}$
UV Slope (β; fλ ∝ λβ )	$-{1.9}_{-0.1}^{+0.2}$	$-{2.3}_{-0.2}^{+0.2}$
Dust Attenuation (A5500 Å )	${0.3}_{-0.2}^{+0.4}$	${0.1}_{-0.1}^{+0.2}$
Dust Attenuation (A1500 Å )	${0.8}_{-0.5}^{+0.7}$	${0.3}_{-0.2}^{+0.4}$
Age (t50/Myr)	${163}_{-133}^{+20}$	${111}_{-83}^{+26}$
SFR50 Myr (M⊙ yr−1)	${10}_{-5}^{+17}$	${6}_{-2}^{+5}$
reff (kpc)	0.7	0.5
Sérsic Index n	0.8	1.0

Note. SED fitting assumes a continuity prior on the star formation history and a Chabrier (2003) initial mass function (IMF).

Download table as:  ASCIITypeset image

The stellar mass for both objects is constrained to be ≈ 10⁹ M_⊙, comparable to GN-z11 (Oesch et al. 2016; Johnson et al. 2021; Tacchella et al. 2022). We have verified the stellar mass is stable to changes in the star formation history prior by also testing the "bursty" prior from Tacchella et al. (2022), which allows more rapid fluctuations in the SFH from time bin to time bin than in the fiducial model. The star formation rates averaged over the last 50 Myr (SFR₅₀) are typical for galaxies of comparable mass at z ≈ 7–10 (e.g., Stefanon et al. 2022a). The SEDs are consistent with negligible dust attenuation and have blue UV slopes, β ≲ − 2. We note that all these derived properties from the SED fits are collectively consistent with a z > 10 interpretation for these galaxies.

Galaxy–galaxy lensing may be particularly important at the redshift frontier where flux-limited surveys may be preferentially sampling magnified sources (e.g., Wyithe et al. 2011). Here we make a simple estimate of how lensed our sources are by assuming their neighbors are singular isothermal spheres (e.g., Fort & Mellier 1994; Schneider et al. 2006; Treu 2010) following, e.g., McLure et al. (2006), Oesch et al. (2014), and Matthee et al. (2017). For this estimate, galaxy redshifts and stellar masses are based on our EAZY fits. Velocity dispersions that trace the underlying dark matter halos are inferred from the stellar mass by extrapolating the empirical scaling relation in Zahid et al. (2016), which is fit to z < 0.7 quiescent galaxies that span M_⋆ ≈ 10⁹–10¹² M_⊙. We choose a local relation to cover the low masses relevant to the most massive neighbors at < 10'' that are likely to produce significant magnification, while noting that the redshift evolution of such relations at least for M_⋆ ≳ 10¹¹ M_⊙ is expected to be gradual—e.g., ≈ 20% higher dispersion at fixed stellar mass at z ≈ 2, (e.g., Mason et al. 2015b).

For both GL-z10 and GL-z12 we find negligible lensing (μ < 1.1) from all foreground sources within 100. GL-z12 has two relatively massive M_⋆ ≈ 10⁹ M_⊙ neighbors apparent in the bottom-left quadrant of the stamps in Figure 1. Even for these two nearby neighbors (z ≈ 2 at a separation of 08, and z ≈ 3 at 20) the lensing is modest. We further note that GL-z13 has a compact morphology (Section 4.5) that does not show elongation along any particular direction that would hint at strong magnification. Based on these considerations we conclude that the observed luminosities of GL-z10 and GL-z12 are likely to be their intrinsic luminosities.

We fit the sizes of both candidates in the F444W imaging (λ_rest ∼ 3500 Å) using GALFIT (Peng et al. 2010). We create 100 pixel cutouts around each galaxy, then use photutils and astropy to create a segmentation map to identify nearby galaxies. We simultaneously fit any sources with magnitudes (estimated from the segmentation map) up to 2.5 mag fainter than the target galaxy that have centers within 3'' of the galaxy; we mask fainter or more distant galaxies. In our fits, we constrain the center of the target galaxy to be within 10 pixels (04) of the input value, the Sérsic index n to be between 0.01 and 8, the magnitude to be between 0 and 45, and the half-light radius r_e to be between 0.3 and 200 pixels (0012–80). We calculate and subtract off a scalar sky background correction from each cutout, estimated from the masked, sigma-clipped cutout, then fix the sky background component in GALFIT to zero. We use a theoretical PSF model generated from WebbPSF at our 004 pixel scale; we oversample the PSF by a factor of 9 in order to minimize artifacts as we rotate the PSF to the GLASS observation angle calculated from the Astronomer's Proposal Tool file, then convolve with a 9 × 9 pixel square kernel, and downsample to the mosaic resolution.

We find reliable Sérsic fits for both galaxies, with half-light radii of 0.5 and 0.7 kpc, respectively, and disklike profiles (n = 1 and n = 0.8, respectively). The models are shown in Figure 4, and the size, and Sérsic profile estimates are listed in Table 3.

Zoom In Zoom Out Reset image size

The resulting sizes of 0.5 and 0.7 kpc are typical for luminous L^* galaxies at z ∼ 6–9, where measurements have been possible to date (e.g., Holwerda et al. 2015; Shibuya et al. 2015; Bowler et al. 2017; Kawamata et al. 2018; Yang et al. 2022). They are also consistent with expectations from simulations for z > 9 galaxies (e.g., Marshall et al. 2022; Roper et al. 2022). However, at z ∼ 7, the most luminous sources often break up in multiple clumps (Bowler et al. 2017). This is not the case for these two sources, at least down to the resolution limit of order 500 pc for the F444W bandpass. Interestingly, GL-z10 even shows tantalizing evidence for being an ordered disk galaxy at z ∼ 10, based on the exponential light profile and elliptical morphology. If we interpret GL-z10's projected axis ratio of 0.65 using a sample of randomly oriented axisymmetric oblate rotators (following, e.g., Holden et al. 2012; Chang et al. 2013; van der Wel et al. 2014) and adopt c/a ≤ 0.4 as a threshold for disks, we find that the observed axis ratio implies P(disk) ∼ 0.5. Our analysis shows the unparalleled power of JWST to provide accurate profile measurements of early universe galaxies.

The key caveat, as well as the key animating spirit of this work, is that these data are among the first deep extragalactic fields collected by a new great observatory. Systematic uncertainties (e.g., zero-point corrections, treatment of artifacts) can still be significant. We have tested for zero-point offsets by comparing HST and JWST photometry for brighter sources where possible and have not found any major issues at the ≲10% level. Nevertheless, we have attempted to account for remaining uncertainties with conservative choices—e.g., a 10% error floor on all fluxes and focusing on bright galaxies whose > 2 mag breaks are robust to even major uncertainties.

A next caveat applies to the SED models underpinning the stellar population parameters and photo-z fitting. Important details about the nebular emission and nature of massive stars at low metallicities, which dominate the light in these few 100 Myr old star-forming systems, remain unconstrained (e.g., Stanway et al. 2020). These uncertainties directly translate to the parameters we recover from SED fitting and the sources for which we are able to fit high quality redshifts. Fortunately, for our redshift range of interest, from the perspective of continuum fitting for the stellar population analyses, extreme nebular emission from strong rest-frame optical lines is shifted out of all NIRCam bands and the most important feature for the photo-zs in the bright galaxies we study is the Lyman break. Further, there exists a range of plausible, but hitherto unconstrained physical ingredients that are unaccounted for in our models (e.g., primordial AGN, top-heavy IMFs, superluminous Population III stars; Windhorst et al. 2018; Pacucci et al. 2022; Steinhardt et al. 2022). Some of these ingredients may potentially produce large UV luminosities in the absence of substantial stellar mass. Spectroscopic follow-up is therefore essential to confirm the nature of these sources.

These two galaxies enable a first estimate of the number densities of relatively luminous sources at z = 10–13. Given the large uncertainties in redshift, we conservatively estimate a selection volume across this full redshift range, i.e., Δz = 3. Given that the depth of all the fields we studied is far deeper than the bright magnitude of these two sources, these galaxies would have been identified over the full area, without foreground galaxies. This amounts to a volume of 2.2 × 10⁵ Mpc³. The detection of two galaxies with M_UV = − 21 thus results in an estimated UV LF point of $\mathrm{log}\phi [{\mathrm{Mpc}}^{-3}\,{\mathrm{mag}}^{-1}]=-{5.05}_{-0.45}^{+0.37}$. In Section 4.4 we conclude that this UV LF point likely reflects the intrinsic luminosity of the sources, and is unaffected by lensing.

In Figure 5, we compare this estimate with previous UV LF determinations and with extrapolations from lower redshifts. In particular, we show that an extrapolation of the Schechter function trends estimated at z = 3–10 results in an LF at M_UV = − 21 that is a factor >10×lower than our estimate. Interestingly, however, when extrapolating the trends in the double power-law LFs from Bowler et al. (2020) to z ∼ 11.5 (the mean redshift of our sources), we find relatively good agreement. In fact, GN-z11 also lies on this extrapolated LF. However, this would indicate very little evolution in the bright galaxy population at z > 8. Indeed, our estimate is in good agreement with previous z ∼ 10 UV LF determinations and constraints at the bright end from Oesch et al. (2018), Bouwens et al. (2019), Stefanon et al. (2019), Morishita et al. (2020), Finkelstein et al. (2022), Bagley et al. (2022), and Leethochawalit et al. (2022).

Zoom In Zoom Out Reset image size

Finally, we also briefly compare our estimates with simulated LFs from the UniverseMachine model (Behroozi et al. 2019) and from the Delphi model (Dayal et al. 2014, 2022). While our estimate is in rough agreement with these predictions at z ∼ 11, the model LFs evolve very rapidly at these early times, such that the z ∼ 12 LF is already >30× below our estimate. This is a general trend of model predictions: a relatively rapid evolution of the LF at z > 10, driven by the underlying evolution of the dark matter halo mass function (see also Oesch et al. 2018; Tacchella et al. 2018; Bouwens et al. 2021). However, the handful of bright galaxies that have been found at z ∼ 10–13 to date appear to oppose this trend.

It is still unclear what the physical reason for this might be. Combined with the discovery of GN-z11 (Oesch et al. 2016), and taking our SED fits at face value (see Section 5.1 for caveats), evidence is mounting that the star formation efficiency in the early universe may be much higher than expected (e.g., Tacchella et al. 2013, 2018; Mason et al. 2015a) in at least a few sources, thus resulting in the early appearance of UV-luminous galaxies with stellar masses as high as ≈ 10⁹ M_⊙ already a few 100 Myr after the Big Bang. The existence of these massive galaxies at such early times raises interesting questions about just how early such galaxies began forming, potentially earlier than current expectations. Wider area data sets will be required to increase the search volume, for more reliable constraints on the number densities of luminous sources.