Insight from JWST /NIRCam into galaxy overdensities around bright Ly 𝛼 emitters during reionization: implications for ionized bubbles at 𝑧 ∼ 9

Several studies have detected Lyman-alpha (Ly 𝛼 ) from bright ( 𝑀 uv ≲ − 21 . 5) galaxies during the early stages of reionization despite the significantly neutral intergalactic medium. To explain these detections, it has been suggested that 𝑧 > 7 Ly 𝛼 emitters (LAEs) inhabit large, physical Mpc (pMpc)-scale ionized regions. However, systematic searches for the overdensities of faint galaxies expected to be powering these ionized bubbles around LAEs have been challenging. Here, we use CEERS JWST /NIRCam imaging to investigate the possibility of large-scale galaxy overdensities associated with two very UV-bright LAEs at 𝑧 = 8 . 7 in the EGS field. We design a colour selection to identify objects at 𝑧 = 8 . 4 − 9 . 1, selecting 28 candidates (including the one LAE in the footprint, EGSY8p7). We model the SEDs of these objects and infer that all are moderately faint ( − 21 . 2 ≲ 𝑀 uv ≲ − 19 . 1) with stellar masses of 𝑀 ∗ ≈ 10 7 . 5 − 8 . 8 M ⊙ . All are efficient ionizing agents ( 𝜉 ∗ ion ≈ 10 25 . 5 − 26 . 0 Hz erg − 1 ) and are generally morphologically simple with only one compact ( 𝑟 𝑒 ≲ 140 to ∼ 650 pc) star-forming component. Of the 27 objects besides EGSY8p7, 13 lie within 5 arcmin of EGSY8p7, leading to a factor of four galaxy overdensity in projection at separations < 5 arcmin ( ∼ 1 . 4 pMpc in projection at 𝑧 = 8 . 7). Meanwhile, separations of 10 − 15 arcmin ( ∼ 2 . 7 − 4 . 1 projected pMpc) are consistent with an average field. The spatial distribution of our sample may qualitatively suggest a large ( 𝑅 ≥ 2 pMpc) ionized bubble encompassing both LAEs in the field, which is theoretically unexpected, but may be possible for a galaxy population four times more numerous than the average to create given moderate escape fractions ( 𝑓 esc ≳ 0 . 15) over long times ( ≳ 200 Myr). Upcoming spectroscopic follow up will enable characterization of the size of any ionized bubble that may exist and the properties of the galaxies powering such a bubble


INTRODUCTION
During cosmic reionization, nearly all hydrogen atoms in the intergalactic medium (IGM) were ionized by radiation emerging from the earliest galaxies.Reionization is thus inextricably linked to the formation and evolution of these early galaxies, and significant effort over the past two decades has been dedicated to understanding this connection (e.g.Madau et al. 1999;Fan et al. 2006;Robertson et al. 2015).
Recently, attention has turned to understanding the cause of Ly  emission observed from galaxies at  > 7 (e.g.Vanzella et al. 2011;Oesch et al. 2015;Zitrin et al. 2015;Castellano et al. 2016Castellano et al. , 2018;;Endsley et al. 2021b;Tilvi et al. 2020;Jung et al. 2020Jung et al. , 2022Jung et al. , 2023;;Larson et al. 2022;Tang et al. 2023;Bunker et al. 2023;Saxena et al. 2023).It is expected that Ly  emission from galaxies at these epochs is highly attenuated by the partially neutral IGM, so the detection of Ly  is somewhat unexpected.One possible explanation for strong observed Ly  is that these LAEs are embedded in large ionized bubbles; if a galaxy is surrounded by a large volume of ionized gas, the Ly  it emits can cosmologically redshift far into the damping wing before encountering intergalactic neutral hydrogen, thereby preserving much of the intrinsic emission escaping the galaxy (Wyithe & Loeb 2005;Furlanetto & Pritchard 2006;Weinberger et al. 2018).Such large ionized bubbles could result from significant overdensities of galaxies, where copious amounts of ionizing photons are produced by an unusually abundant galaxy population (e.g.Barkana & Loeb 2004;Furlanetto et al. 2004;Iliev et al. 2006;Dayal & Ferrara 2018;Qin et al. 2022).However, it is also possible that Ly  may be detectable from galaxies in relatively small ionized regions if, for example, the majority of Ly  photons that emerge from the galaxy are already significantly redshifted from systemic (as has frequently been observed for luminous sources; e.g.Erb et al. 2014;Willott et al. 2015;Mason et al. 2018b;Hashimoto et al. 2019;Matthee et al. 2020;Endsley et al. 2022) or the galaxy has a very intense radiation field that facilitates the production and visibility of Ly  (Stark et al. 2017;Endsley et al. 2021b;Simmonds et al. 2023;Tang et al. 2023;Roberts-Borsani et al. 2023).Thus, it is unclear if galaxies observed to be emitting strong Ly  reside in small bubbles primarily powered by their own ionizing emission, or large bubbles generated by the combined ionizing flux of all neighbouring galaxies contributing to the reionization process.
Observations with the Hubble Space Telescope (HST) and groundbased facilities delivered the first insights into the nature of overdensities around UV-luminous LAEs at  > 7. Using dedicated HST/Wide Field Camera 3 (WFC3) imaging, Leonova et al. (2022) found an excess of moderately UV-bright ( uv ≲ −20), photometrically selected galaxies surrounding three extremely UV-luminous ( uv ≈ −22) LAEs in the Extended Groth Strip (EGS) field at  = 7.5 (Roberts-Borsani et al. 2016;Stark et al. 2017),  = 7.7 (Oesch et al. 2015), and  = 8.7 (Zitrin et al. 2015) on projected scales of  ∼ 0.3 physical Mpc (pMpc) (i.e. the physical scale of the WFC3 field of view at these redshifts, ∼ 1 arcmin).Additional spectroscopic efforts with Keck/MOSFIRE have revealed Ly  emission from several moderately bright (−20.5 <  uv < −20.0) galaxies < 0.7 pMpc away from the extremely UV-luminous LAE at  = 7.7 (Tilvi et al. 2020; see also Jung et al. 2022 and JWST/Near Infrared Spectrograph observations by Jung et al. 2023), and at  = 8.7, a second bright LAE only 4 pMpc away from the already known LAE in the field (Zitrin et al. 2015) was confirmed with MOS-FIRE (Larson et al. 2022).Both of these discoveries lent additional evidence that these rare systems occupied significant galaxy overdensities associated with large ionized bubble.Furthermore, the two  = 8.7 LAEs are also accompanied by an overdensity of bright ( uv ≲ −20.8)HST-selected photometric candidates over the entire EGS field (Finkelstein et al. 2022a;Larson et al. 2022) -some of which are now spectroscopically confirmed to lie at  = 8.7 (Tang et al. 2023;Nakajima et al. 2023;Arrabal Haro et al. 2023a) -which may further suggest that these systems exist in a large ionized bubble.
HST could identify signatures of overdensities of moderately bright galaxies relatively close to bright LAEs and larger scale overdensities of brighter galaxies.However, a key question moving forward is whether overdensities of faint galaxies regularly extend to large,  ≳ 1 pMpc spatial scales.Significant large-scale overdensities associated with LAEs may be suggestive of large ionized bubbles being powered by many faint neighbours of LAEs, whereas a lack of such structures may be indicative of smaller ionized bubbles and physics internal to the LAEs facilitating Ly  escape (e.g.outflows leading to large Ly  velocity offsets).To help distinguish these scenarios, it will thus be important to (1) map faint galaxies in regions surrounding extremely UV-luminous LAEs to separations of a few pMpc and (2) obtain constraints on the density and ionizing properties of these faint ( uv ≳ −20) neighbouring galaxies.Fortunately, JWST enables both of these observational studies.First, JWST allows for the photometric identification of faint neighbouring galaxies of LAEs, placing the first constraints on their spatial distribution and ionizing properties.Then, once faint neighbours have been identified, spectroscopic follow up targeting Ly  can begin to characterize the sizes and morphologies of ionized regions in the vicinity of the bright LAEs.
In this work, we characterize the overdensity potentially associated with two of the highest redshift known, extremely UV-bright LAEs, EGSY8p7 (Zitrin et al. 2015) and EGS_z910_44164 (Larson et al. 2022), which lie in the EGS field at  = 8.7.We use deep JWST/Near Infrared Camera (NIRCam; Rieke et al. 2023) imaging from the Cosmic Evolution Early Release Science (CEERS) program1 (Finkelstein et al. 2022b, Finkelstein et al. in prep), which extends up to ∼ 15 arcmin away from EGSY8p7 and ∼ 27 arcmin away from EGS_z910_44164, allowing us to search for faint neighbouring systems at large projected separations.Ultimately, this enables us to place constraints on the spatial extent and amplitude of overdensities around the LAEs using a large sample of galaxies, then discuss implications for any ionized bubble that may be present.
This paper is organized as follows.In Section 2, we describe the JWST and ancillary HST imaging used in this work, as well as our method for measuring photometry.We present our selection and sample in Section 3, then discuss implications for a galaxy overdensity and ionized bubble in Section 4. Finally, we summarize our conclusions in Section 5. Throughout this work, we adopt a flat ΛCDM cosmology with ℎ = 0.7, Ω m = 0.3, and Ω Λ = 0.7.All reported uncertainties correspond to the marginalized 68 per cent credible interval and all logarithms are base-10 unless stated otherwise.Finally, all magnitudes are given the AB system (Oke & Gunn 1983) and all distances are physical unless stated otherwise.

IMAGES AND PHOTOMETRY
In this work, we wish to investigate the large-scale neighbourhoods of two extremely UV-luminous ( uv ≈ −22) Ly  emitters at  = 8.7 (Zitrin et al. 2015;Larson et al. 2022) in the EGS field.
To achieve this, we combine deep multi-band JWST/NIRCam and HST/Advanced Camera for Surveys (ACS) imaging over EGS to perform a dropout selection designed to identify sources at redshifts similar to the LAEs (see Section 3 for details).In this section, we describe the imaging data and photometric measurements used for our analysis.
The JWST/NIRCam imaging utilized in this work were taken as part of the CEERS survey in June and December 2022.These NIRCam data consist of imaging in six broadband filters (F115W, F150W, F200W, F277W, F356W, and F444W) and one medium band filter (F410M) over ten independent pointings for a total area of ∼ 85 arcmin2 .At the redshifts of interest for this work (around  ∼ 8.7), the NIRCam filters cover a rest-frame wavelength range of  rest ≈ 1000 − 5100 Å, probing both the rest-UV (largely redward of the Ly  break) and the rest-optical (including the strong emission lines H and [O iii]4959,5007 up to  ∼ 8.9).
We produce co-added mosaics for each NIRCam band following the methods described by Endsley et al. (2023b).We begin from the individual uncalibrated NIRCam exposures (*_rate.fits)and process these exposures with the JWST Science Calibration Pipeline, 2 utilizing the most recent NIRCam photometric calibration reference files (first released as jwst_0942.pmapin early October 2022).We align the astrometry of the co-added images for each filter, module, and pointing combination to the Gaia astrometric reference frame using the tweakreg package, adopting the Gaia-matched mosaic from the HST/WFC3 F160W Complete Hubble Archive for Galaxy Evolution (CHArGE) project (Kokorev et al. 2022) as a reference image.This alignment achieves an rms offset of ∼ 6 − 15 mas relative to the HST mosaics; see Chen et al. (2023) for details on the astrometric calibration.During this process, we also implement a 1/  noise mitigation algorithm and global background subtraction using the sep package (Barbary 2016) as described by Endsley et al. (2023b).Finally, all co-added mosaics for each filter are sampled to the same world coordinate system at a resolution of 30 mas pixel −1 .
To assist in identifying  ∼ 8.7 Lyman break galaxies across the CEERS NIRCam footprint, we also utilize optical imaging taken with HST/ACS.The EGS field has been observed in the F435W, F606W, and F814W ACS bands as part of the All-Wavelength Extended Groth Strip International Survey (AEGIS; Davis et al. 2007), the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011;Koekemoer et al. 2011), and the Ultraviolet Imaging of the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey Fields (UVCANDELS3 ; PI: Teplitz) program.The ACS mosaics used in this work were produced using the grizli software (Brammer et al. 2022) as part of the CHArGE project.All HST mosaics are registered to the Gaia astrometric frame and the same world coordinate system with pixel scale of 40 mas pixel −1 .
The HST and JWST imaging introduced above span more than a factor of three in angular resolution, ranging from a point spread function (PSF) full width at half maximum (FWHM) of ≈ 0.045 arcsec in ACS/F435W to ≈ 0.15 arcsec in NIRCam/F444W.To ensure robust colour calculations, we must therefore account for variations in the PSF across different images.However, we simultaneously wish to preserve as much signal-to-noise as possible, particularly in the shorter wavelength filters that probe the Ly  break.We thus homogenize the PSFs of all the ACS and NIRCam short wavelength (SW; F115W, F150W, F200W) mosaics to that of ACS/F814W, while the NIRCam long wavelength (LW; F277W, F356W, F410M, F444W) mosaics are homogenized to the PSF of WFC3/F160W to enable future studies using HST/WFC3.
We identify sources across the CEERS NIRCam footprint by running Source Extractor on an inverse variance weighted stack of the PSF-homogenized F150W and F200W mosaics.The variance images are computed by running sep on each mosaic and thus incorporate both background and detector noise.We calculate photometry following procedures previously used for HST-based analyses of  ≳ 6 galaxies (e.g.Bouwens et al. 2015Bouwens et al. , 2021;;Finkelstein et al. 2015Finkelstein et al. , 2022a;;Endsley et al. 2021a).Our methods are described in full by Endsley et al. (2023b), but we provide a brief overview here.First, we use the background-subtracted mosaics to compute photometry in elliptical Kron (1980) apertures following procedures similar to those used for previous analyses of  ≳ 6 galaxies (e.g.Bouwens et al. 2015Bouwens et al. , 2021Bouwens et al. , 2023;;Finkelstein et al. 2015Finkelstein et al. , 2022aFinkelstein et al. , 2023;;Endsley et al. 2021a), adopting aperture sizes set by a Kron parameter of  = 1.2 to maximize signal-to-noise (Finkelstein et al. 2022a).Photometric uncertainties for each object in the catalog are calculated as the standard deviation of flux measured in apertures of the same size as the aperture used for the object, which are randomly placed in nearby object-free regions as determined from a Source Extractor segmentation map specific to each band.Finally, we perform aperture corrections in two stages to obtain total fluxes (see section 2.3 of Finkelstein et al. 2022a for details).First, we multiply all ACS and SW fluxes and their errors by the ratio of F200W flux measured in apertures with  = 2.5 versus  = 1.2 from the PSF-homogenized F200W mosaic for ACS and the NIRCam SW filters, and an inverse variance weighted stack of the LW mosaics for the LW photometry.Second, the photometry is corrected for flux outside the  = 2.5 aperture using the measured PSFs of ACS/F814W and WFC3/F160W from the mosaics, normalizing the encircled energy distribution at large radii (> 1 arcsec) based on reported values in the instrument manuals.To account for the different PSFs of the ACS+SW and LW images, we multiply the aperture size of the ACS+SW bands by a factor of 1.5 to obtain the aperture for the LW photometry, where the factor of 1.5 reflects the typical size ratio of twelve UV-bright  ∼ 6 − 8 EGS galaxy candidates (Chen et al. 2023) measured from the PSF-homogenized F200W and F444W mosaics.
Though the NIRCam images have a relatively high angular resolution, it is nevertheless possible for flux from neighbouring objects to contaminate our photometric apertures.To mitigate the possibility of such contaminated photometry, we have developed a neighbour subtraction algorithm heavily based on techniques developed for Spitzer/Infrared Array Camera (IRAC) deconfusion (e.g.Labbé et al. 2010Labbé et al. , 2013;;Bouwens et al. 2015;Endsley et al. 2021a) and described in detail by Endsley et al. (2023b).In brief, given an object in the catalog, we fit Sérsic profiles to neighbours identified in the Source Extractor catalog using the non-PSF matched mosaics, accounting for the PSF of each filter during the fit.Using the bestfitting Sérsic profile parameters, we then subtract the neighbouring objects in the PSF-homogenized images, accounting for the changes in the PSF.We apply this algorithm to sources in our colour selected sample (selected as described in Section 3.1) that are contaminated by flux from a neighbouring object within the  = 2.5 Kron aperture as determined by a visual inspection.

SAMPLE SELECTION AND PROPERTIES
We now wish to build a sample designed to assess whether there is evidence of any large-scale galaxy overdensities surrounding the two extremely UV-luminous ( uv ≈ −22) Ly  emitters at  = 8.7 in the EGS field.These LAEs lie at nearly the same redshift (EGSY8p7 at  sys = 8.678 and EGS_z910_44164 at  sys = 8.610; Tang et al. 2023;Larson et al. 2023) and are separated by only 4 pMpc, suggesting they may be part of the same overdensity powering a very large ( ≳ 2 pMpc) ionized bubble (Larson et al. 2022).
To investigate the possibility of a large-scale  ∼ 8.7 overdensity within EGS, we identify Lyman break galaxies at similar photometric redshifts over the entire CEERS NIRCam imaging footprint.This footprint extends to approximately ∼ 15 (∼ 27) arcmin away from EGSY8p7 (EGS_z910_44164), corresponding to maximum separations of ∼ 4 (7.5)pMpc in projection at  = 8.7.As the number of  ∼ 9 galaxy candidates in EGS that have thus far been targeted with sensitive rest-optical spectroscopy is limited, in this work we aim to quantify the overdensity using photometrically selected Lyman break galaxies (as in, e.g.Castellano et al. 2016;Endsley et al. 2021b;Leonova et al. 2022).Consequently, we focus on angular separations, acknowledging that the full 3D separations can be larger given the redshift uncertainties associated with broadband dropout selections (e.g. a redshift separation of Δ = 0.1 corresponds to a line-of-sight separation of ∼ 2.6 pMpc at  = 8.7).
We first require our sample to have apparent magnitudes of F200W < 28, S / N > 7 in F200W, and S / N > 3 in at least two LW NIRCam filters to ensure the selected sources are real (criteria (i) and (ii)).Next, criteria (iii), (iv), and (v) are designed to ensure the presence of a complete Ly  break between ACS/F814W and NIRCam/F150W and effectively no flux blueward of the break, as expected at  > 8 from the very high Lyman series and Lyman continuum opacity through the partially neutral IGM (Inoue et al. 2014).We note that fluxes in the F814W and F115W dropout filters are set to their 1 upper limits if detected at < 1 significance.Criterion (vi) is designed to identify partial F115W dropouts while selecting against a complete Ly  break in F115W to place an upper limit on our redshift range.Finally, criteria (vii) and (viii) are designed to remove candidates with extremely red spectral energy distributions (SEDs) in order to allow us to distinguish  > 8 sources, which are typically blue (e.g.Bouwens et al. 2014;Bhatawdekar & Conselice 2021;Topping et al. 2022b;Cullen et al. 2023), from intrinsically red galaxies at lower redshifts.We note that there is no ACS/F435W imaging for ∼ 15 arcmin 2 of the total ∼ 85 arcmin 2 NIRCam imaging footprint.For these areas, we adjust criterion (iii) to require S/N < 2 only in F606W and F814W and adjust criterion (iv) to  2 opt < 3.5, as we only calculate  2 opt using two filters.In general, we expect this selection to efficiently identify moderately faint objects at redshifts of  ∼ 8.4 − 9.1.We estimate our selection completeness by generating mock power law SEDs (   ∝   ) with rest-UV continuum slope of  = −2.We normalize these SEDs to a range of F200W magnitudes (F200W = 24 − 30, d = 0.05) at a variety of redshifts ( = 5 − 12, d = 0.02), applying the IGM attenuation model of Inoue et al. (2014).For our fiducial completeness models, we assume a Ly  EW of 0 Å, but explore the impact of EW Ly  > 0 Å below.We then convolve the normalized SEDs with the transmission functions of our ten ACS and NIRCam filters to obtain mock observations, perturb the 'observed' fluxes appropriately for the depth of each filter as measured from real sources, and apply the selection criteria described above.Using these models, we estimate that our colour selection is ≥ 50 per cent (≥ 80 per cent) complete from  ∼ 8.4 − 9.1 ( ∼ 8.5 − 9.1) at F200W = 27 with a maximum efficiency of ∼ 90 per cent due to upscattering in the optical ACS bands and our requirements of ACS nondetections and  2 opt < 5.At fainter magnitudes, F200W ∼ 27.5, our completeness function broadens in redshift as objects at slightly lower ( ∼ 8.1 − 8.4) and higher ( ∼ 9.1 − 9.4) redshifts scatter into our selection with ≲ 10 per cent completeness.
If the galaxies we wish to select reside in an overdensity, and associated large ionized bubble, we may expect that they have reasonably strong Ly  emission.Thus, we wish to ensure that our selection would identify such objects at redshifts similar to those of the two known LAEs.At  ∼ 8.7, Ly  transmits through F115W, resulting in bluer F115W − F150W colours at fixed redshift if a galaxy has strong Ly  emission, ultimately shifting our selection to higher redshifts for LAEs.To quantify this effect, we model our selection completeness adopting Ly  EWs ranging from 5 Å to 50 Å.For Ly  EW Ly  = 10 Å (approximately the median value for  uv ≲ −20.5 galaxies at  ∼ 7 found by Endsley et al. 2021b), our selection identifies moderately faint (F200W < 28) objects at redshifts of  ∼ 8.5 − 9.3.For stronger Ly  (EW Ly  = 30 Å), the selection identifies candidates at redshifts of  ∼ 8.6 − 9.4, still near the redshift of the LAEs.In particular, we highlight that our selection recovers EGSY8p7 itself (EGS-69787 in our catalog), which is at a systemic redshift of  sys = 8.678 (Tang et 2. The observed and model SEDs of a subset of our sample: three objects with spectroscopic redshifts (top row), and the six objects closest to EGSY8p7 (all of which are within 3 arcmin of EGSY8p7; center and bottom rows).The observed photometry is shown as teal circles, with open circles denoting 2 upper limits, and the model photometry is shown as purple diamonds.For the objects without a spectroscopic redshift, we also show the photometric redshift posteriors that result from the models restricted to photometric redshifts of  phot = 4 − 12 described in Section 3.3.We highlight that most of these objects have red observed F410M − F444W colors (F410M − F444W ≳ 0.5), which is likely due to strong [O iii]+H  emission boosting the F444W flux; a 0.5 mag F444W excess approximately corresponds to an EW of ∼ 650 Å at  = 8.7.observations, Tang et al. 2023;Larson et al. 2023; upper bound from ground-based Keck/MOSFIRE observations, which is impacted by an OH sky line, Zitrin et al. 2015).We would also expect our selection to identify EGS_z910_44164 at  sys = 8.610 (Tang et al. 2023) with EW Ly  ≈ 5 Å (Larson et al. 2022), though this source falls outside the CEERS NIRCam footprint.
To further ensure that our conclusions do not depend strongly on the possible lack of the strongest LAEs (EW Ly  ≳ 75 Å) at  ∼ 8.7 identified by our primary selection, we design a secondary selection to identify these objects, described in Appendix A. Using this selection, we identify 13 additional candidates over the field, one of which lies 0.7 arcmin away from EGSY8p7.We highlight that, in general, our conclusions do not significantly change if we supplement our fiducial sample with these potential strong LAEs at  ∼ 8.7.
After applying our fiducial selection criteria, we conduct a visual inspection of the objects that are identified.We remove sources that (1) have photometry that is likely contaminated by a nearby bright source or artifact (primarily diffraction spikes or 'snowballs' caused by a significant cosmic ray event), (2) appear to have by-eye detections in any of the three optical ACS bands (which are expected to be essentially completely attenuated by the IGM at our redshifts of interest; Inoue et al. 2014), or (3) appear to have unreliable noise measurements (often due to being near a detector edge).We additionally remove one object (EGS-25508) with a pattern of photometric excesses in F150W, F200W, and F277W, followed by a drop in flux in the longer wavelength filters that is highly suggestive of emission lines at  ∼ 2 boosting the 1.5 − 3 m filters, and for which the SED models described in Section 3.2 find no acceptable solutions at  ≥ 6.
After visual inspection, we obtain a final colour-selected sample of 27 dropout candidates (one of which is EGSY8p7), which we use as our fiducial sample for the remainder of this work.However, for the interested reader, we summarize the objects that we removed during our visual inspection in Appendix B. For each of the 27 objects in our fiducial sample, we perform a second visual inspection to determine whether the  = 2.5 Kron apertures used for Table 1.Observational and physical properties of the candidates identified with our colour and S/N selection criteria (Section 3.1) ordered by increasing angular separation from EGSY8p7.We report the observed F200W magnitudes and observed rest-UV continuum slopes.We also report properties inferred from our beagle models: the photometric redshifts (or systemic spectroscopic redshift without errors, if available), stellar masses, sSFRs, and intrinsic ionizing photon production efficiencies before processing through dust and gas.For these SED models, we have assumed that these objects lie at  ≥ 4 to ensure that our inferences are representative of the physical properties these systems would have at high redshift rather than the properties these objects would have at lower redshifts; see Section 3.3.Before presenting our sample in detail, we investigate whether  ∼ 9 candidates previously identified with HST imaging in EGS satisfy our selection criteria.Bouwens et al. (2021) reported four objects with photometric redshifts of  phot = 8.2 − 9.2 (EGSY-9597051195, EGS910-3, EGS910-10, and EGS910-15) that fall within the CEERS imaging footprint, three of which we recover (one of which is EGSY8p7, and the other two of which are now confirmed to lie at  spec = 8.7; Tang et al. 2023;Arrabal Haro et al. 2023a).We find that the fourth, EGSY-9597051195, has only a weak break between F115W and F150W (F115W − F150W = 0.1) and is consequently best fit by a photometric redshift of  phot = 7.2, lower than we consider in this work.Finkelstein et al. (2022a) also reported four HST-selected objects at  phot = 8.2 − 9.2 that lie within the CEERS imaging (EGS_z910_6811, EGS_z910_68560, EGS_z910_26890, and EGS_z910_40898), all but one of which we select (noting that EGS_z910_6811 is EGSY8p7 and EGS_z910_26890 is the same object as EGS910-15 identified by Bouwens et al. 2021).EGS_z910_68560 is not detected within 1 arcsec of its expected position in the NIRCam imaging; given its reported F160W magnitude of  ab = 25.8 reported by Finkelstein et al. (2022a) and the depth of the JWST imaging (5 depth of ∼ 28.3 mag in F150W), this suggests that this object may be a transient source.
Additionally, Leonova et al. (2022) conducted a targeted HST/WFC3 search for galaxy candidates that lie very close to extremely luminous, high redshift galaxies in the EGS field.They identified eight objects, including EGSY8p7 and EGS910-3, and our selection recovers four (EGSY8p7, EGS910-3, EGSY-66990, and EGSY-69674).One of the candidates that we do not identify falls in a detector gap in the SW NIRCam imaging of CEERS (EGSY-27727 in Leonova et al. 2022).Another (EGSY-56680, with reported F160W magnitude of  ab = 27.2) is not detected in the CEERS imaging which, as for EGS_z910_68560, may imply that this object could be a transient source.Finally, the remaining two objects we do not select (EGSY-6887 and EGSY-12580) both have bluer F115W − F150W colours than we allow during our colour selection (F115W − F150W = 0.4 and −0.2, respectively) as well as evidence of F410M excesses consistent with strong rest-optical emission lines transmitting through F410M at slightly lower redshifts than we target in this work ( ∼ 7 − 7.5).
We emphasize that the identification of objects by HST that are likely at slightly lower redshifts than we target is unsurprising.Photometric redshift uncertainties with HST photometry can be large (e.g.Leonova et al. 2022 reported median upper uncertainties of Δ ≈ 0.2 and median lower uncertainties of Δ ≈ 0.5).Furthermore, while objects down to  ∼ 7.5 drop out at least partially in WFC3/F105W, objects at  ∼ 7.5 − 8 do not significantly drop out in NIRCam/F115W, allowing us to cleanly remove these slightly lower redshift objects from our sample.Overall, the availability of the NIR-Cam/F115W filter significantly improves our ability to distinguish between objects at  ∼ 7.5 − 8 and  ∼ 9.Moreover, overlapping JWST/NIRCam broad and medium bands then enable us to constrain the photometric redshifts of candidates via strong emission lines even more precisely.

Photometric redshifts
We next wish to assess whether the candidates we have identified may lie at redshifts similar to that of the two known bright LAEs.We therefore infer the photometric redshifts of the objects in our sample with the BayEsian Analysis of GaLaxy sEds tool (beagle; Chevallard & Charlot 2016).beagle computes the stellar and nebular emission of galaxies by combining the latest version of the Bruzual & Charlot (2003) stellar population synthesis models with the cloudy photoionization code (Ferland et al. 2013); see description by Gutkin et al. (2016).The updated Bruzual & Charlot (2003) models are underpinned by single star stellar isochrones computed by the PAdova and TRieste Stellar Evolution Code (parsec; Bressan et al. 2012;Chen et al. 2015).For these models, we adopt a Chabrier (2003) initial mass function (IMF) with mass range 0.1 − 300 M ⊙ , an SMC dust prescription (Pei 1992), and the IGM attenuation model of Inoue et al. (2014).We place a uniform prior on redshift ranging from  phot = 0 − 15 and assume a constant star formation history (CSFH) with a log-uniform prior on age ranging from 1 Myr to the age of the Universe at the redshift under consideration.We adopt uniform priors on logarithmic stellar mass (5 ≤ log( * /M ⊙ ) ≤ 12), band optical depth (−3 ≤ log( v ) ≤ 0.7), stellar metallicity (−2.2 ≤ log( * /Z ⊙ ) ≤ 0.24), and ionization parameter (−4 ≤ log() ≤ −1).Finally, we assume that the total interstellar gas-and dust-phase metallicity is equal to the stellar metallicity with dust-to-metal mass ratio to   = 0.3, noting that beagle self-consistently models the effects of depletion.
We find compelling  ∼ 8 − 10 solutions for all objects but two in our colour-selected sample.EGS-25508 has no acceptable solutions at  phot ≥ 6, so we remove this object from our sample.EGS-64770 is best fit at  phot = 7.7, and while we retain this object in our sample, we caution that it may be at a slightly lower redshift than we target in this work.For the remainder of the sample, the median  phot = 8 − 10 SEDs inferred by our models generally fit the data well within uncertainties.If the SED model fits the observed photometry exactly within uncertainties, we would expect  2 = 9 or  2 = 10 depending on the filters available at the location of the object, 5 and we find that the  phot = 8 − 10 SEDs have  2 values of 1.3 − 19.4,indicating that we find reasonable fits at  = 8 − 10.
Our goal for these models is to assess the feasibility of our colour selected candidates lying at redshifts near to those of the known LAEs in the EGS field.We find that the large majority of our sample (18/23, or 80 per cent of the objects without spectroscopic redshifts) have very large probabilities (≥ 95 per cent) of lying at  > 6, and twelve of those also have large, ≥ 90 per cent probabilities of lying in the more narrow redshift range of  = 8.2 − 9.2.However, we acknowledge that in some cases, there are non-negligible probabilities of lower redshift ( ∼ 2) solutions that are at times degenerate with the  ∼ 8 − 10 solutions.These lower redshift solutions are often extremely dusty with extreme emission lines, with typical stellar masses of  * ≈ 10 7−8 M ⊙ , typical ages of a few hundred Myr, stellar metallicities of  * ∼ 0.1 − 0.2 Z ⊙ , and strong dust attenuation ( v ∼ 0.7, where  v = 1.086 v ).These degeneracies underline the need for spectroscopic confirmation of photometric candidates to rigorously assess the neighbourhoods of bright LAEs during reionization.Nevertheless, we emphasize that, though these lower redshift solutions exist, we still find credible  ∼ 8 − 10 SEDs for all candidates in the sample.Thus, we proceed with our entire colour selected sample to evaluate any potential overdensity associated with the LAEs.

Sample properties
The ionizing, star-forming, and morphological properties of the  ∼ 8.7 galaxies in our sample grant insight into not only their potential contribution to reionization, but also how these early galaxies first assembled.We therefore wish to infer these physical properties from the photometric SEDs we have measured (see Figure 2 for example SEDs from our sample).To do this, we use beagle to re-fit the photometric data with a very similar model set up as our models designed to infer photometric redshifts for selection described in Section 3.2, but restrict the redshift to  phot = 4 − 12 (uniform prior) and logarithmic stellar metallicity to −2.2 ≤ log(/ ⊙ ) ≤ −0.3 (uniform prior).We have set an upper limit on stellar metallicity of  ∼ 0.5 Z ⊙ to avoid solutions with unreasonably high, near-solar metallicities for our  ∼ 8.7 galaxies.For objects with known redshifts (EGSY8p7/EGS-69787, EGS-83141, EGS-87873, EGS-109389, and EGS-110201), we fix the redshift to the systemic redshift.We report selected physical properties inferred from these SED models, as well as selected observational properties, for each object in our sample in Table 1.
The objects in our sample, excluding EGSY8p7, have observed F200W magnitudes ranging from F200W = 26.2− 28.0 (median F200W = 27.3),implying that these objects are inherently moderately faint if they lie at our redshifts of interest,  ∼ 8.7.We also find that most objects have moderately blue observed rest-UV continuum .The inferred stellar masses, sSFRs, and intrinsic stellar ionizing photon production efficiencies (before stellar photons are processed through dust and gas) of our sample as a function of absolute UV magnitude with error bars showing formal uncertainties from the SED models, which may be a factor of a few larger when considering systematics introduced by e.g.altering the stellar IMF or model SFH.We also infer the properties of the two bright LAEs in the field, shown as colored points.The red square is EGSY8p7 (Zitrin et al. 2015) and the blue-green diamond denotes EGS_z910_44164 (Larson et al. 2022).We note that we adopt the HST and Spitzer photometry reported by Finkelstein et al. (2022a) for EGS_z910_44164, as it is not covered by the CEERS NIRCam imaging.
At  ∼ 8.7, the ≳ 3 m imaging from JWST probes rest-frame optical emission (see SEDs in Figure 2), allowing us to gain detailed insights into the star-forming and ionizing properties of the sample.For example, if F356W, F410M, and F444W are all brighter than the bluer filters, this may empirically suggest the presence of a Balmer break built up by a mature stellar population.Indeed, one object in our sample (EGS-91290) is ∼ 0.4 − 0.8 mag brighter in F356W, F410M, and F444W than in the bands that probe the rest-UV continuum, and is best fit with an old CSFH age of ∼ 200 Myr (implying a formation redshift of  ∼ 12 assuming the median photometric redshift of this object,  phot = 8.71) and a moderately large stellar mass of ∼ 10 8.8 M ⊙ (assuming a CSFH).
Moreover, at  ≲ 8.9, the strong rest-optical emission lines H  and [O iii]4959,5007 transmit through the F444W filter.Thus, a red F410M − F444W color implies both the presence of strong [O iii]+H  emission and a redshift of  ≲ 8.9.Specifically, at  = 8.7, F410M − F444W = 0.5 suggests EW [O iii]+H  ≈ 640 Å, slightly less than the median EW of ∼ 770 Å expected for  ∼ 6.5 − 8 galaxies at similar luminosities found by Endsley et al. (2023b).In our 6 Calculated by fitting   ∝   to the observed F150W, F200W, and F277W fluxes, i.e. the filters expected to probe the rest-UV that are not attenuated by the IGM at our redshifts of interest. 7Calculated by integrating the beagle model spectra over a value unity tophat from rest-frame wavelengths of  rest = 1450 − 1550 Å. sample, we find that nine of the 27 objects have F410M − F444W ≥ 0.5, consistent with the very large [O iii]+H  EWs expected for this population (e.g.De Barros et al. 2019;Endsley et al. 2021aEndsley et al. , 2023b)).Of the remaining 18 objects, seven have F410M We next investigate whether this distribution of F410M − F444W colors is consistent with expectations for [O iii]+H  emission in this population.Following the methods of Gehrels (1986) assuming binomial statistics, we estimate the fraction of the sample that we would expect to observe with F410M − F444W < 0.5 given the [O iii]+H  EW distribution for  uv ≲ −19 galaxies at  ∼ 6.5 − 8 found by Endsley et al. (2023b).To a 1 confidence level, we find that we would expect to observe up to 14 of 27 objects with F410M − F444W ≲ 0.5 if all are at  ≲ 8.9, while we observe 18 in our sample.However, we note that we expect to observe more objects with F410M − F444W ≲ 0.5 colours if any objects in our sample are at redshifts of  ≳ 8.9, as the stronger [O iii]5007 component of the [O iii]4959,5007 doublet quickly redshifts out of F444W at  ≳ 8.9.Thus, the large number of objects lacking strong F444W excesses could suggest that some of the objects in our sample may lie at slightly higher redshifts of  ∼ 8.9 − 9.1, as would potentially be expected from our selection completeness (Figure 1).Alternatively, it may be possible that more objects than expected are members of the population of weak [O iii]+H  emitters observed in photometric samples of galaxies with similar luminosities at  ∼ 6−8 (e.g.Endsley et al. 2023b,a), as well as at least one spectroscopically confirmed at  = 8.807 with an [O iii]5007 EW of ∼ 370 Å (Fujimoto et al. 2023).Such weak line emission could potentially suggest several physical situations: (1) a high occurrence rate of low metallicities ( ≲ 0.01 − 0.05 Z ⊙ , which have been observed in large numbers in photometric  ≳ 6 galaxy samples from JWST; e.g.Endsley et al. 2023b,a), leading to weaker [O iii] emission than is found in more chemically evolved systems; (2) a large fraction of galaxies being observed after a rapid downturn in star formation rate (SFR) over the most recent ∼ 10 Myr, leading to a relatively small population of OB stars to produce strong nebular line emission; or (3) large Lyman continuum escape fractions, leading to a decrease in the interactions between hydrogen-ionizing photons and the interstellar medium that produce nebular lines (e.g.Zackrisson et al. 2013).For these weak line emitters, our beagle models with a CSFH and fixed  esc = 0 generally tend to find solutions at  ∼ 8.5 − 8.9 (driven primarily by the strength of the observed Ly  break) with low metallicities rather than solutions at higher redshifts.More model flexibility in quantities such as the escape fraction or the SFH may be justified to better fit the large population of weak line emitters and could ameliorate the models' need to fit these systems with low metallicity solutions, but ultimately, upcoming future spectroscopy will be necessary to confirm the redshifts and emission line properties of these objects to precisely characterize their properties and spatial distribution.
We now turn to investigating insights provided by the NIRCam photometry into the physical properties of the objects in our sample.In general, we find properties consistent with expectations for  ∼ 9 galaxies previously studied with HST and JWST (see Figure 3 for a summary of the properties we infer).Under the assumptions of a CSFH, we find moderate stellar masses of  * ≈ 10 7.5−8.8M ⊙ (median  * ≈ 10 8.0 M ⊙ ).As expected for the fainter luminosities we can probe with JWST, these masses are smaller than the masses inferred for the two LAEs and other HST-selected objects at the same redshifts (e.g.Tacchella et al. 2022), but are broadly consistent with inferences for  ≳ 6 galaxies at similar luminosities studied with JWST (e.g.Endsley et al. 2023b;Bradley et al. 2022;Leethochawalit et al. 2023;Santini et al. 2023;Furtak et al. 2023;Whitler et al. 2023a;Robertson et al. 2023).We also find that the objects in our sample generally appear to be dominated by recent bursts of star formation, with large CSFH-based specific star formation rates (sSFR; sSFR ≈ 5 − 600 Gyr −1 with median 42 Gyr −1 ), and are very efficient ionizing agents ( * ion ≈ 10 25.5−26.0Hz erg −1 , median 10 25.6 Hz erg −1 , where  * ion is the intrinsic stellar ionizing photon production efficiency before processing through dust and gas).These values are similar to the large SED-inferred sSFRs and ionizing photon production efficiencies for galaxies at  ∼ 6 − 8 (e.g.Endsley et al. 2021a;Stefanon et al. 2022;Endsley et al. 2023b).
We note that throughout our SED modelling, we have assumed a Chabrier ( 2003) stellar IMF and a constant SFH.Adopting an alternative IMF could change the inferred stellar masses by ∼ 0.2−0.4dex (Wang et al. 2023;Woodrum et al. 2023) and increase uncertainties by a factor of ∼ 2 − 3 over formal errors from SED models.An alternative form of the SFH model (e.g. a non-parametric model that allows significant stellar mass formation at early times) may lead to larger stellar masses than the CSFH model, typically by factors of up to ∼ five for objects such as those in our sample with CSFH ages of a few tens of Myr (e.g.Whitler et al. 2023b;Endsley et al. 2023a).The form of the SFH may also alter the ionizing photon production efficiencies inferred; for example, a two-component parametric SFH (a delayed- SFH at early times and a constant component at recent times that is decoupled from the delayed- model; see Endsley et al. 2023b,a) infers smaller ionizing photon production efficiencies by less than a factor of three for the majority our sample.Ultimately, these effects may combine to increase uncertainties by factors of a few over formal SED model uncertainties, though many of these model choices may have competing effects on the inferred galaxy properties and thus quantifying systematic effects will require a better understanding of e.g. the true IMF and the typical galaxy SFH during reionization.
In addition to the SED-based properties we infer, the high angular resolution of the JWST/NIRCam imaging allows us to examine the morphological properties of the sample.As shown in Figure 4, the majority of the faint objects we have identified in this work are less morphologically complex than the very bright EGSY8p7.They often have only one star-forming component while EGSY8p7 has three, consistent with findings from previous studies that clumpy morphologies become more common at brighter luminosities (e.g.Bowler et al. 2017;Chen et al. 2023).Following the methods of Chen et al. ( 2023), we derive the half-light radii,   , of these objects by fitting their F150W surface brightness profiles at original resolution (i.e.before PSF homogenization) with fixed  = 1 Sérsic profiles (i .e. exponential profiles) using the lmfit package (Erwin 2015).For the three components of EGSY8p7, we assume an exponential profile for each component and fit them simultaneously, and find that each of its three star forming complexes are extremely compact (  ≲ 140 pc in the rest-UV, upper limit corresponding to the F150W pixel scale of 30 mas) and are each separated by ∼ 900 pc.We also find that the single components of the faint galaxies in our sample are often relatively compact.Specifically, we measure half-light radii in the rest-UV of the five brightest single-component candidates in our sample range from   ≲ 140 pc to ∼ 650 pc.This implies star formation rate surface densities (Σ sfr = (SFR/2)/( 2  )) ranging from Σ sfr ∼ 6 M ⊙ yr −1 kpc −2 up to Σ sfr ≳ 80 M ⊙ yr −1 kpc −2 , which are larger than observed in  ∼ 0 galaxies (e.g.Shibuya et al. 2015) but similar to those found for other systems at  ∼ 9 (Ono et al. 2022).
Using the resolved NIRCam imaging, we can also investigate the possibility for a spatially extended or offset component of old stars in EGSY8p7 itself, whose light would be outshined by any young stars present in the integrated photometry (for a more detailed discussion of this effect in the general UV-bright population, see Chen et al. 2023).If such a component were to exist, this would imply that EGSY8p7 may have been contributing to the beginning stages of the growth of an early ionized bubble.However, we find no clear evidence that any component of EGSY8p7 is dominated by a mature (≳ 100 Myr) stellar population.All three of the star-forming complexes labelled in Figure 4  Figure 5.The observed surface density of the objects in our sample as a function of angular separation from EGSY8p7.On the top -axis, we convert this angular separation to a projected physical separation from EGSY8p7 at its systemic redshift,  sys = 8.678 (Tang et al. 2023;Larson et al. 2023).At the closest separations we consider, < 5 arcmin (corresponding to ∼ 1.4 pMpc at  = 8.678), we find a surface density of (0.45 ± 0.12) arcmin −2 , a factor of ∼ four overdense relative to expectations derived from the evolution of the UV LF (Bouwens et al. 2021) after accounting for the completeness of our selection.This overdensity decreases with increasing distance away from EGSY8p7 to a factor of ∼ 1.9 overdense at separations of 5 − 10 arcmin (∼ 1.4 − 2.7 pMpc) until the densities are consistent within Poisson uncertainties with an average field at separations of 10 − 15 arcmin (∼ 2.7 − 4.1 pMpc).and 7.7.This suggests that EGSY8p7 is dominated by stars formed in a recent burst of star formation, and if there is an older population of stars, they are outshined by the recently formed population.
However, though the rest-UV and optical SED is likely dominated by a young stellar population formed in a burst of star formation, it may be possible that the SED at longer wavelengths is more sensitive to the presence of a more mature stellar population.It has been demonstrated by Papovich et al. (2022) that, relative to pre-JWST constraints from Spitzer/IRAC data, the addition of MIRI data covering the rest-optical SED of EGSY8p7 leads to significantly more precise constraints on the stellar mass.Here, we investigate whether the addition of longer wavelength MIRI data can significantly influence quantities such as the SFH and stellar mass when starting from the NIRCam SED for EGSY8p7.We note that the F410M NIRCam filter already provides a constraint on the rest-frame optical continuum ( rest ≈ 3980 − 4430 Å), though MIRI probes much longer rest-frame wavelengths (up to  rest ∼ 8900 Å with F770W).
To conduct this investigation, we repeat our beagle SED modelling analysis assuming a CSFH, adding the MIRI F560W and F770W photometry reported by Papovich et al. (2022) to the NIRCam SED.In general, we find similar properties with the MIRI data as without: CSFH-based stellar mass of log( * /M ⊙ ) ∼ (8.8−9.0)±0.1 and age of ∼ 4−6 Myr, suggesting that the NIRCam imaging probing the rest-optical continuum may be sufficient to constrain the mass of the stellar population dominating the SED.However, we also investigate whether the addition of MIRI photometry enables us to place stronger upper limits on the early SFH than is possible with NIRCam alone.If so, it may be possible to more precisely quantify the contribution of ionizing photons from early star formation in EGSY8p7 towards the early growth of an ionized bubble.
To this end, we model the NIRCam and MIRI SED of EGSY8p7 with a non-parametric SFH model using a prior that not only allows, but tends to prefer an extended early period of star formation.We follow the methods of Whitler et al. (2023b), assuming the nonparametric 'continuity' prior of the Bayesian SED modelling tool Prospector (Leja et al. 2019;Johnson et al. 2021b) and adopting eight age bins for the SFH, with the two most recent age bins fixed to 0 − 3 Myr and 3 − 10 Myr and the rest spaced evenly in logarithmic time.With this model, we infer a factor-of-six larger stellar mass (log( * /M ⊙ ) ∼ 9.7) than the CSFH model, both with and without MIRI.Ultimately, this suggests that for objects similar to EGSY8p7, with very young light-weighted ages that suggest a recent burst of star formation, the SED continues to be dominated by nebular continuum emission well into the rest-frame optical.Thus, even with photometry probing wavelengths as red as  rest ∼ 8900 Å, significant systematic uncertainties in the stellar masses of burst-dominated objects can still exist, which may ultimately require rest-near-IR observations or dynamical masses to resolve (e.g.Tang et al. 2022).This stellar mass difference is caused by a long, early period of early star formation (SFR ∼ 10 M ⊙ yr −1 , leading to a half-mass age of ∼ 160 Myr) in the non-parametric model, which suggests that EGSY8p7 could have feasibly been contributing to the reionization process for several hundred Myr.Moreover, there is evidence that EGSY8p7 may host an active galactic nucleus (AGN) (Mainali et al. 2018;Larson et al. 2023), which could have undergone past periods of rapid growth in the past despite being a sub-dominant component of the rest-UV and optical emission at the time of observation.In a similar manner as past star formation activity, past AGN activity could also contribute copious amounts of ionizing photons towards the growth of a large ionized bubble at early times (Larson et al. 2023), though such activity may also be short-lived or episodic.

A GALAXY OVERDENSITY NEAR UV LUMINOUS Ly 𝛼
EMITTERS AT  = 8.7 With our sample of galaxies selected to lie near  = 8.7, we now turn to evaluating the potential for a pMpc-scale overdensity and ionized bubble associated with the two  = 8.7 LAEs in the EGS field.Using our NIRCam-based sample and SEDs, we can extend previous HST-based investigations of the environments of the LAEs (e.g.Finkelstein et al. 2022a;Leonova et al. 2022) to both fainter luminosities and larger physical scales.The depth of the NIRCam imaging allows us to identify objects as faint as  uv ∼ −19 (∼ 0.1 * uv ), probing the faint galaxy population thought to provide the bulk of ionizing photons needed to create an ionized bubble.Furthermore, the imaging extends up to ∼ 15 arcmin away from EGSY8p7 and ∼ 27 arcmin away from EGS_z910_44164 (physical separations up to ∼ 4 pMpc and ∼ 7 pMpc in projection at  = 8.7, respectively), allowing us to search for faint galaxies over pMpc scales that may trace correspondingly large ionized bubbles.We concentrate primarily on the galaxy population near EGSY8p7, as EGS_z910_44164 does not fall in the CEERS NIRCam imaging footprint, but also comment briefly on the distribution of galaxies between the two LAEs.We focus on the sample colour selected as described in Section 3.1, as the inclusion of the objects identified by our LAE-targeted selection (Appendix A) does not significantly impact our conclusions.
To quantify the amplitude and spatial extent of any overdensity that may exist, we first must determine the number of objects expected in an average field given our selection function and the UV luminosity function (UV LF).To achieve this, we assume a Schechter parametrization of the UV LF and the Bouwens et al. ( 2021) redshift evolution of the Schechter parameters, then convolve with our  et al. 2023;Fujimoto et al. 2023;Arrabal Haro et al. 2023b).Objects from the literature are shown as grey diamonds (Bouwens et al. 2021), squares (Finkelstein et al. 2022a), and triangles (Leonova et al. 2022).We note that we do not show objects from the literature if the NIRCam photometry indicates they do not lie in our targeted redshift range; see the discussion in Section 3.1.The CEERS NIRCam imaging footprint is outlined in grey.Finally, the candidates we identify in this work are shown as filled, coloured points -stars if they are spectroscopically confirmed (Tang et al. 2023;Fujimoto et al. 2023;Arrabal Haro et al. 2023a) and circles if they are not -and the two bright  = 8.7 LAEs (Zitrin et al. 2015;Larson et al. 2022) are shown as large stars.In the left panel, we colour our sample by the redshift, either the median photometric redshift inferred by our SED models (see Section 3.3) or the systemic spectroscopic redshift if known.In the right panel, the colour bar shows the object's absolute UV magnitude inferred from our SED models.Besides EGSY8p7 itself, we identify 26 candidates over the entire field (ten of which are within the ∼ 20 arcmin 2 of the imaging that lies within 4 arcmin of EGSY8p7), many more than the ten expected for an average field.Furthermore, most of the sample lies in the northeastern region of the footprint, relatively close to EGSY8p7 and EGS_z910_44164.
analytic selection completeness (described in Section 3.1 and shown in Figure 1).We find an average surface density of 0.11 arcmin −2 in the redshift and magnitude range probed by our selection. 8Thus, we would expect to identify ten objects over the entire ∼ 85 arcmin 2 field.However, we have identified ten candidates within only a four arcminute radius of EGSY8p7 (∼ 20 arcmin 2 of the imaging area), or  ∼ 1 pMpc in projection at  = 8.7, potentially suggesting at least a mild overdensity in the close vicinity of EGSY8p7.
We can now compare this expected surface density to the observed surface density implied by our sample in the field of the  = 8.7 Ly emitters.As we are interested in assessing overdensities that may be associated with EGSY8p7, we calculate the surface densities in concentric circular annuli centred on EGSY8p7 with radii increasing in steps of Δ = 5 arcmin (corresponding to Δ = 1.4 projected pMpc at  = 8.7).These annuli intersect with 30 arcmin 2 , 41 arcmin 2 , and 16 arcmin 2 of the CEERS imaging area for separations < 5 arcmin, 5 − 10 arcmin, and 10 − 15 arcmin, respectively.In marked contrast to the 0.11 arcmin −2 surface density expected for an average field, we find a surface density of (0.45 ± 0.12) arcmin −2 (uncertainties calculated assuming Poisson statistics) in the 30 arcmin 2 of the imaging 8 Assuming the UV LF parametrizations and  ∼ 9 parameters derived by other studies (e.g.Harikane et al. 2023b;Donnan et al. 2023) implies generally similar surface densities ranging from ∼ 0.07 − 0.12 arcmin −2 .that lies within  = 5 arcmin of EGSY8p7 -a factor of four overdensity in projection.As the distance from EGSY8p7 increases, the surface density decreases to (0.30 ± 0.08) arcmin −2 at separations of  = 5−10 arcmin, then (0.13±0.09) arcmin −2 at  = 10−15 arcmin (consistent with the surface density of an average field). 9We show these observed surface densities in comparison with the expected surface density in Figure 5.In particular, we highlight that, though the region within 5 arcmin of EGSY8p7 comprises only a third of the total imaging area, nearly half of our photometric candidates (including two of the four spectroscopically confirmed objects) lie within that area; see Figure 6.These objects are all significantly fainter and inferred to be much less massive than EGSY8p7.They also all have large sSFRs (∼ 20 − 600 Gyr −1 ) and large intrinsic ionizing photon production efficiencies ( * ion ∼ 10 25.6−26.0Hz erg −1 ) that are typical of the full sample (see Section 3.3).
While the NIRCam imaging does not probe the volume immediately surrounding EGS_z910_44164, we can use the imaging between EGSY8p7 and EGS_z910_44164 to qualitatively assess whether there is any evidence for a single large ionized region containing both LAEs.If the LAEs do inhabit the same ionized bubble, their ∼ 4 pMpc separation would require an ionized region with a radius of  ∼ 2 pMpc, at minimum.We note that, in a spherical geometry, this would place the LAEs very close to the edges of the bubble, near intergalactic H i where Ly  may be easily attenuated.This would suggest that an even larger ionized region ( ≳ 3 pMpc), at least along the line of sight, may be necessary to truly facilitate the transmission of Ly  through the IGM.However, even an  = 2 pMpc ionized bubble would be larger than theoretically expected given the large neutral fraction expected at  ∼ 9 (e.g.Lin et al. 2016;Geil et al. 2016;Lu et al. 2023), so we consider an  = 2 pMpc ionized bubble as the most conservative case.
If there is one large ionized region encompassing both LAEs, we may expect a larger surface density of galaxies between the LAEs than in the areas furthest away from them.Indeed, we find that ∼ 80 per cent of the sample falls in the five most northeastern NIRCam pointings of the CEERS imaging, approximately between EGSY8p7 and EGS_z910_44164 (∼ 6 − 15 arcmin away from EGS_z910_44164, or ∼ 1.7 − 4 pMpc in projection at  = 8.7), while only three are in the most southwestern region of the imaging, ∼ 20 − 25 arcmin (∼ 6 − 7 projected pMpc) away from EGS_z910_44164 (see Figure 6).At face value, this could suggest a single large,  ≳ 2 pMpc-scale ionized region encompassing both EGSY8p7 and EGS_z910_44164.Such a structure could also potentially extend along the line of sight, as our photometric selection extends up to  ∼ 9 and the fraction of objects with weak F444W excesses in our sample may imply that at least some of them could lie at  ∼ 8.9 − 9.1 (∼ 5 − 10 pMpc along the line of sight from  = 8.7); see discussion in Section 3.3.
To investigate the plausibility of the presence of a large-scale ionized bubble, we estimate the size of the H ii region that could be created by the galaxy population we have explored in this work.Adopting the methods described by Endsley & Stark (2022) (see also Shapiro & Giroux 1987;Haiman & Loeb 1997;Cen & Haiman 2000), we approximate the radius of a spherical H ii region,  H ii , that a given population of galaxies within  = 2 pMpc of EGSY8p7 (corresponding to a sphere of volume 33.5 pMpc 3 ) could create using the following equation: The first two terms of Equation (1) describe the growth of ionized bubbles due to photoionizations from the ionizing sources and the Hubble flow, respectively, while the last term accounts for recombinations in the IGM.At the redshift, , corresponding to the time step under consideration,  () is the Hubble parameter and ⟨ H i ()⟩ is the average neutral hydrogen density. b is the case B recombination coefficient, where we have assumed a temperature of 10 4 K, leading to a value of  b = 2.59 × 10 −13 cm 3 s −1 (Osterbrock & Ferland 2006). = ⟨ 2 H i ⟩/⟨ H i ⟩ 2 is the 'clumping factor' that accounts for inhomogeneities in the IGM (Madau et al. 1999).Motivated by theoretical expectations from simulations for the average IGM (e.g.Shull et al. 2012;Finlator et al. 2012;Kaurov & Gnedin 2015;Gorce et al. 2018), we assume a fixed value of  = 3, though we note that if the region is overdense, the clumping factor may be higher (e.g.Mao et al. 2020;Bianco et al. 2021), leading to smaller bubbles due to the increased recombination rate; we briefly present predicted bubbles sizes with clumping factors of  = 10,  = 30, and  = 50 in Appendix C. For  ion , we assume the median ionizing photon production efficiency of our sample inferred from our SED models,  ion = 10 25.7 Hz erg −1 .We note that this quantity is the ionizing photon production efficiency after the processing of stellar photons through dust and gas has occurred, in contrast to the previously re- The colourbar shows the predicted radius of a spherical ionized region,  H ii , as a function of the length of time the galaxies have been producing ionizing photons (i.e. the age) and the escape fraction of hydrogenionizing Lyman continuum photons (  esc ).We assume an ionizing photon production efficiency of  ion = 10 25.7 Hz erg −1 and show the predicted bubble sizes if the number density is that of an average volume in the left panel, and if the volume is a factor of four overdense in the right panel.To ionize an  ≥ 2 pMpc spherical volume, galaxies in an average field would need very high escape fractions that are rarely seen even in individual objects with Lyman continuum detections.However, a galaxy population that is a factor of four overdense may be able to power an  ≥ 2 pMpc bubble with only moderate escape fractions (  esc ∼ 0.15) if they are producing ionizing photons for several hundred million years.ported intrinsic stellar ionizing photon production efficiency,  * ion .To calculate the non-ionizing UV luminosity,  uv , we again assume the same Schechter parametrization of the UV LF as we adopted to estimate the expected surface density in the field, and consider the contributions of all galaxies  uv = −15 and brighter (approximately the faint end turnover of the UV LF tentatively measured in lensed fields, Bouwens et al. 2022).To obtain  uv in a field that is a factor of  overdense, we assume that the overall shape of the UV LF is the same as that in an average field, but with a higher normalization by a factor of  (i.e.we multiply the normalization of the Schechter function by ).Finally, we integrate Equation (1) over time assuming that all galaxies in the region have been producing ionizing photons at the same rate with the same escape fraction for their entire lifetime.
Figure 7 shows the predicted radii of H ii regions for a variety of ages and escape fractions of ionizing photons (  esc ) assuming an average field (left panel) and a field that is a factor of four overdense (right panel).Under the assumptions described above, we conclude that for  uv ≤ −15 galaxies in an average field to power an ionized bubble large enough to encompass both LAEs within the age of the Universe, they would all require escape fractions of  esc ≳ 0.6, significantly larger than expected for the population (and are rare even in individual objects with Lyman continuum detections at  ≲ 3; e.g.Izotov et al. 2018;Pahl et al. 2021;Flury et al. 2022).
However, if the volume hosts a factor of four galaxy overdensity,  uv ≤ −15 galaxies may be able to create an  ≥ 2 pMpc bubble with an escape fraction of  esc ≈ 0.15 over ≳ 200 Myr (approximately a third of the age of the Universe at  = 8.7), or with slightly higher escape fractions over slightly shorter times (e.g. esc ≈ 0.25 for 100 Myr).Such escape fractions, while somewhat high, are found in lower redshift samples and are consistent with inferences from indirect tracers for  esc for galaxies during reionization (e.g.Meyer et al. 2020).Ultimately, this suggests that the galaxy population we have considered in this work, combined with the fainter galaxies under our detection limit, may be able to create a large enough ionized region to contain both LAEs over long times.
However, we highlight that though the escape fractions may be physically plausible, the galaxy population we have considered in this work would require upwards of a few hundred Myr to carve out an  ≥ 2 pMpc bubble.From our CSFH models, we have inferred that the SEDs of our sample generally tend to be dominated by much shorter, recent bursts of star formation.Specifically, we have inferred a median CSFH sSFR of∼ 42 Gyr −1 , which implies star formation on timescales of only a few tens of Myr (Section 3.3) -insufficient time to carve out an  ≥ 2 pMpc bubble under our assumptions for Equation ( 1).If such a bubble does exist, this could imply that the amplitude of the galaxy overdensity around EGSY8p7 is stronger than a factor of four (which may be possible if most or all of the photometric candidates in this work are at redshifts very close to  = 8.7).It could also imply galaxies fainter than  uv = −15 are contributing to the process of creating the bubble, or that there is simply a larger population of faint galaxies (i.e. the faint-end slope of the UV LF is steeper than we have assumed here).Alternatively, we have shown that a much older stellar population formed during an extended early period of star formation can be hidden under the young stellar population dominating the SED for burst-dominated objects such as those in this work (Section 3.3; see also Tacchella et al. 2022;Topping et al. 2022a;Whitler et al. 2023b).Such an extended early SFH, if it lasted for several hundred Myr, could enable the creation of a large ionized bubble, and may in fact be necessary for the objects in our sample to carve out an  ≥ 2 pMpc bubble.We note that the average SFH from a variety of simulations (e.g.Dawoodbhoy et al. 2018;Ocvirk et al. 2020;Legrand et al. 2022) are qualitatively consistent with an early period of increasing star formation activity.Thus, we may indeed expect galaxies to have started creating a large ionized bubble at relatively early times, though a detailed accounting of ionized bubble growth would require careful tracking of  esc ,  ion , and  uv over time, so we defer a comprehensive examination of the impact of the SFH on ionized regions to future work.
If an  ≥ 2 pMpc ionized bubble exists, it may be a rare and individual occurrence.Theoretical predictions by Lu et al. (2023) suggest that even an  = 2 pMpc bubble is very unlikely, even if the region hosts a factor of 3 − 5 photometric galaxy overdensity (consistent with both the findings of Leonova et al. 2022 and the overdensity we measure in the inner ∼ 1.4 projected pMpc around EGSY8p7).Alternatively, it could be indicative that the reionization topology at  = 8.7 is characterized by larger bubbles than expected.We highlight that the small probability of finding a large ionized region at  = 8.7 is, in large part, due to the very large volumeaveraged fraction of intergalactic neutral hydrogen expected at  ∼ 9 ( H i ∼ 0.8−0.9;Mason et al. 2019b;Ocvirk et al. 2021;Rosdahl et al. 2022;Lu et al. 2023;Mitra & Chatterjee 2023;Mutch et al. 2024).If the IGM neutral fraction is lower, then ionized bubbles are expected to be larger given the same dominant source population of ionizing photons.Alternatively, at fixed neutral fraction, ionized bubbles can be larger (though also more rare) if brighter, more massive galaxies are the dominant source of ionizing photons (e.g.Mesinger et al. 2016;Geil et al. 2016;Lu et al. 2023).If an  ≥ 2 pMpc ionized bubble exists at all, it would facilitate significant transmission of Ly  through the IGM for both the known bright LAEs and their faint neighbours (of which there may be > 2 − 4 within a 2 × 2 arcmin 2 area around the central LAE; Qin et al. 2022).Thus, upcoming deep JWST/NIRSpec observations (JWST GO-4287; PI: Mason) aimed at characterizing the ionizing properties, spatial distribution, and Ly  line properties (e.g.velocity offsets) of any faint LAEs that may exist in this volume will be important to distinguish these scenarios.

SUMMARY
In this work, we have used JWST/NIRCam imaging from the CEERS program (Finkelstein et al. 2023, Finkelstein et al. in prep) in the EGS field to investigate the potential for an overdensity of faint galaxies in the vicinity of two extremely UV-luminous LAEs at  = 8.7.We have built a sample of photometrically selected galaxy candidates, quantified the 2D galaxy overdensity suggested by these objects, and examined the implications of such an overdensity for the existence of large ionized bubbles at  ∼ 9. We summarize our key results below.
(i) We have developed a S/N and colour selection designed to identify moderately faint galaxies at  ∼ 8.4 − 9.1, redshifts similar to those of the two known bright LAEs in the EGS field.Using this selection, we have identified 27 objects, including the one LAE that falls within the CEERS imaging footprint, EGSY8p7.
(iii) The faint objects in our sample typically have simple morphologies.They are often composed of only one star-forming clump, in contrast to the very bright EGSY8p7 ( uv ≈ −22.4), which has three very compact star-forming components.The single components in our sample are also typically compact, with rest-UV half-light radii ranging from   ≲ 140 pc to 650 pc.This implies star formation rate surface densities of at least Σ sfr ∼ 6 M ⊙ yr −1 kpc 2 , and as high as Σ sfr ≳ 80 M ⊙ yr −1 kpc 2 (adopting the SFRs inferred from our SED models).
(iv) Besides EGSY8p7 itself, we have identified ten candidates within only  = 4 arcmin of EGSY8p7 (∼ 20 arcmin 2 of the CEERS imaging area) and 26 candidates over the full ∼ 85 arcmin 2 of the imaging.In contrast, for an average field of the same area, we would expect to identify only ten galaxies over the entire field.Compared to the expected surface density derived from the evolution of the UV luminosity function (Bouwens et al. 2021), we find a factor of four overdensity within 5 arcmin (∼ 1.4 projected pMpc at  = 8.7) of EGSY8p7.At increasing separations from the LAE, the surface density of our sample decreases until it is consistent with an average field at 10 − 15 arcmin (∼ 2.7 − 4.1 projected pMpc).
(v) Qualitatively, the spatial distribution of galaxies that we have identified over the field may suggest the existence of a large,  ≥ 2 pMpc bubble that contains both bright LAEs in the field.However, such an ionized bubble may be larger than theoretically expected.We investigate the properties that would be required for a population of galaxies to power an  ≥ 2 pMpc ionized bubble and find that if the  = 2 pMpc spherical volume around EGSY8p7 is a factor of four overdense, the  uv ≤ −15 galaxy population could create this bubble over ≳ 200 Myr with moderate escape fractions of  esc ≈ 0.15.Alternatively, the same galaxy population could power the same ionized bubble in less time, but would require slightly higher escape fractions (e.g. esc ≳ 0.25 over 100 Myr).
(vi) The SEDs of our sample tend to be dominated by light from stars formed in a recent burst of star formation over timescales of a few tens of Myr, shorter than the timescales required for the creation of a large,  ≥ 2 pMpc bubble.This could imply that extended, early SFHs may be necessary in order for the galaxy population we have considered in this work to power a large, pMpc-scale bubble.
(vii) If a large,  ≥ 2 pMpc bubble at  ∼ 9 exists, it may be a rare occurrence or could suggest that ionized regions at  ∼ 9 may be larger than expected.This could potentially be explained if the IGM neutral fraction at  ∼ 9 is lower than current models assume ( H i ∼ 0.8 − 0.9), or if the reionization process is dominated by the rare population of bright, massive galaxies that create larger ionized bubbles (see also Lu et al. 2023).
(viii) Upcoming follow-up spectroscopy of the faint galaxies discussed in this work will be important to better constrain the properties of any ionized bubble, or bubbles, that may be associated with the LAEs.Spectroscopic observations will first enable us to confirm the existence of a galaxy overdensity associated with the LAEs.Then, detailed measurements of Ly  and other rest-UV and optical emission lines will then allow us to study the topology of the bubble(s) in detail and infer the properties of the ionizing galaxies to gain insights into the mechanisms by which the Universe is reionized.

Figure 1 .
Figure 1.Simulated completeness for our colour and S/N selection criteria designed to identify bright objects at redshifts close to  ∼ 8.7 (see Section 3.1).The selection generally identifies F200W = 27 objects at redshifts of  ∼ 8.4 − 9.1 with ≳ 50 per cent completeness.Applying this selection to the CEERS NIRCam imaging, we identify 27 galaxy candidates, one of which is EGSY8p7.
Figure3.The inferred stellar masses, sSFRs, and intrinsic stellar ionizing photon production efficiencies (before stellar photons are processed through dust and gas) of our sample as a function of absolute UV magnitude with error bars showing formal uncertainties from the SED models, which may be a factor of a few larger when considering systematics introduced by e.g.altering the stellar IMF or model SFH.We also infer the properties of the two bright LAEs in the field, shown as colored points.The red square is EGSY8p7(Zitrin et al. 2015) and the blue-green diamond denotes EGS_z910_44164(Larson et al. 2022).We note that we adopt the HST and Spitzer photometry reported byFinkelstein et al. (2022a) for EGS_z910_44164, as it is not covered by the CEERS NIRCam imaging.The objects in our sample are much fainter than either of the LAEs and correspondingly, have smaller stellar masses.Our candidates also have relatively large sSFRs (∼ 5 − 600 Gyr −1 with median 42 Gyr −1 ) and are very efficient ionizing agents (  * ion = 10 25.5−26.0Hz erg −1 with median  * ion = 10 25.6 Hz erg −1 ).

Figure 4 .
Figure 4. Rest-frame UV colour images (red = F200W, green = F150W, blue = F115W) of EGSY8p7 (left, large panel) and the six most luminous objects in our sample (remaining six, smaller panels) including three spectroscopically confirmed: EGS-87873 at  = 8.715 (Tang et al. 2023), EGS-110201 at  = 8.876 (Fujimoto et al. 2023), and EGS-83141 at  = 8.763 (Arrabal Haro et al. 2023a).Our sample is generally much fainter than EGSY8p7, with the brightest object being ∼ 1.2 mag fainter.We highlight the complex morphology of EGSY8p7, which has three distinct star-forming clumps, labelled C1, C2, and C3 in the image.In contrast, the majority of the fainter candidates in our sample have much simpler single-component morphologies.

Figure 6 .
Figure 6.The on sky distribution of the candidates identified by our colour selection criteria (Section 3), along with photometric candidates previously identified HST searches over EGS with reported photometric redshifts of  phot = 8.2 − 9.2 (grey points; Bouwens et al. 2021; Finkelstein et al. 2022a; Leonova et al. 2022), and three spectroscopically confirmed galaxies at  spec = 8.881, 8.998, and 8.638 that are not included in our photometric sample (open stars; Tanget al. 2023;Fujimoto et al. 2023;Arrabal Haro et al. 2023b).Objects from the literature are shown as grey diamonds(Bouwens et al. 2021), squares(Finkelstein et al. 2022a), and triangles(Leonova et al. 2022).We note that we do not show objects from the literature if the NIRCam photometry indicates they do not lie in our targeted redshift range; see the discussion in Section 3.1.The CEERS NIRCam imaging footprint is outlined in grey.Finally, the candidates we identify in this work are shown as filled, coloured points -stars if they are spectroscopically confirmed(Tang et al. 2023;Fujimoto et al. 2023;Arrabal Haro et al. 2023a) and circles if they are not -and the two bright  = 8.7 LAEs(Zitrin et al. 2015;Larson et al. 2022) are shown as large stars.In the left panel, we colour our sample by the redshift, either the median photometric redshift inferred by our SED models (see Section 3.3) or the systemic spectroscopic redshift if known.In the right panel, the colour bar shows the object's absolute UV magnitude inferred from our SED models.Besides EGSY8p7 itself, we identify 26 candidates over the entire field (ten of which are within the ∼ 20 arcmin 2 of the imaging that lies within 4 arcmin of EGSY8p7), many more than the ten expected for an average field.Furthermore, most of the sample lies in the northeastern region of the footprint, relatively close to EGSY8p7 and EGS_z910_44164.

Figure 7 .
Figure 7. Predicted sizes of the ionized bubbles that could be created by the  uv ≤ −15 galaxy population within a 33.5 pMpc 3 volume, corresponding to a sphere of radius  = 2 pMpc that would encompass both  = 8.7 LAEs in the EGS field.The colourbar shows the predicted radius of a spherical ionized region,  H ii , as a function of the length of time the galaxies have been producing ionizing photons (i.e. the age) and the escape fraction of hydrogenionizing Lyman continuum photons (  esc ).We assume an ionizing photon production efficiency of  ion = 10 25.7 Hz erg −1 and show the predicted bubble sizes if the number density is that of an average volume in the left panel, and if the volume is a factor of four overdense in the right panel.To ionize an  ≥ 2 pMpc spherical volume, galaxies in an average field would need very high escape fractions that are rarely seen even in individual objects with Lyman continuum detections.However, a galaxy population that is a factor of four overdense may be able to power an  ≥ 2 pMpc bubble with only moderate escape fractions (  esc ∼ 0.15) if they are producing ionizing photons for several hundred million years.

Figure C1 .
Figure C1.Predicted sizes of the ionized bubble that could be created by a four times overdense population of  uv ≤ −15 galaxies within a sphere of radius  = 2 pMpc.We show three different IGM clumping factors (left:  = 10, center:  = 30, right:  = 50).To ionize an  ≥ 2 pMpc spherical volume with larger clumping factors, galaxies require increasingly large Lyman continuum escape fractions over long times due to increased recombination rates.

Table 2 .
Coordinates, angular separations, and physical separations from EGSY8p7 for our sample, ordered by increasing angular separation from EGSY8p7.If a systemic spectroscopic redshift for a given object is available (denoted by a superscript dagger next to the object ID), we adopt the systemic redshift to calculate the physical separation from EGSY8p7.Otherwise, we adopt the median photometric redshift of the object inferred by beagle.Object has a spectroscopic redshift that was used to calculate physical separation.See references in Table1.