A CEERS Discovery of an Accreting Supermassive Black Hole 570 Myr after the Big Bang: Identifying a Progenitor of Massive z

Introduction

One of the most consequential periods in cosmic history is the Epoch of Reionization (EoR), where the material between galaxies underwent a significant transition from neutral to ionized hydrogen. Satisfactory explanations for the sources of radiation that contributed to this process, when it started, and how long it lasted are all presently wanting, primarily from the lack of necessary observatories and instruments tuned to the early Universe-until now. With the launch of JWST, we are on the cusp of an exciting new era in astronomy, as detailed studies of galaxies in the first billion years are finally possible.

One of the key questions JWST was designed to answer was when and how the first black holes (BHs) formed. Defining this epoch will help constrain the role these sources played in reionization alongside ionizing photons from massive stars. Supermassive black holes (SMBHs) exist at the centers of massive galaxies in the present-day Universe and exhibit a wellstudied correlation with the velocity dispersions (and stellar masses) of galaxy bulge components (see the review by Kormendy & Ho 2013). These massive objects have >13 Gyr of cosmic time to grow to their present-day mass, which is possible via standard accretion scenarios with stellar-mass BH seeds (∼1-10 M e ; possibly up to 100 M e if a Population III star remnant). However, the surprising discovery of z > 6 quasars with the Sloan Digital Sky Survey (SDSS) with BH masses of log (M BH /M e ) > 9 (e.g., Fan et al. 2006;Bañados et al. 2018;Wang et al. 2021;Farina et al. 2022) challenges models of BH growth. Such objects require very early (z ∼ 25-30) stellar-mass seeds with near-unity duty cycles of Eddington-limited accretion and/or super-Eddington accretion to grow to such a mass from a stellarmass seed in  1 Gyr (e.g., Volonteri et al. 2021 and references therein). None of these scenarios are adequately predicted by simulations (e.g., Volonteri et al. 2021;Fontanot et al. 2023). This tension is further exacerbated by the similarly massive quasars now being discovered at z > 7 (e.g., Mortlock et al. 2011;Bañados et al. 2018; see Fan et al. 2022 for a recent review), with the current highest-redshift known active galactic nucleus (AGN) being a bright quasar at z = 7.64 (Wang et al. 2021).

These surprisingly massive BHs in the first Gyr of cosmic history have led to an alternative seeding theory: direct collapse black holes (DCBHs;Bromm & Loeb 2003). In this model, minihalos irradiated by Lyman-Werner photons (11.2-13.6 eV) cannot form molecular gas and thus do not form Population III stars. As these halos grow, they eventually cross the atomic cooling regime, at which point the gas will begin to cool rapidly via H I cooling (primarily Lyα). Further fragmentation-inducing cooling is avoided due to the H 2 -suppressing UV background, leading to collapse into massive BHs in the range of 10 4-6 M e , preceded by a brief phase as a supermassive star (see Smith & Bromm 2019;Woods et al. 2019 and references therein for a detailed explanation of this process). A similarly massive BH seed could be formed as the remnant of a supermassive star powered by WIMP-like dark matter annihilation (so-called "dark stars"; Ilie et al. 2012;Freese et al. 2016). Such massive BH seeds, forming at z ∼ 10-15 (by necessity after the first generation of stars), could more feasibly grow into the observed z ∼ 6-7 quasar population (e.g., Madau et al. 2014;Natarajan et al. 2017;Regan et al. 2019;Latif et al. 2021;Pacucci & Loeb 2022;Massonneau et al. 2023;Trinca et al. 2023). These massive quasars must represent only the extreme cases-there likely exists a much larger population of lower-mass BHs and/or obscured BHs waiting to be discovered.

While DCBHs could alleviate the tension between observed BH masses and our theories of BH growth, such objects have yet to be observed. One clear pathway to understanding BH growth is to observe more SMBHs in the EoR. Identifying modest-sized BHs at earlier cosmic times could provide further evidence as to whether DCBHs are a necessary pathway. Additionally, the discovery of such a population would both better explain how the observed z ∼ 6 quasar population originally formed and inform on the potential contribution of AGNs to reionization both through X-ray heating (e.g., Jeon et al. 2014) and through ionizing photon contributions (e.g., Finkelstein et al. 2019;Giallongo et al. 2019;Grazian et al. 2020Grazian et al. , 2022;;Yung et al. 2021).

Prior to JWST, only the most massive SMBHs at high redshifts could be identified. However, the spectroscopic capabilities of JWST now enable the search for signs of AGN activity from less luminous sources, particularly those embedded in galaxies whose stellar emission dominates the total galaxy luminosity (e.g., Endsley et al. 2022;Kocevski et al. 2023) and/or where the bulk of the accretion emission is obscured (e.g., Fujimoto et al. 2022; see also Furtak et al. 2023).

Here we report the discovery of the first known AGN at z > 8. The galaxy harboring this AGN was first identified as a candidate z ∼ 8 galaxy by Roberts-Borsani et al. (2016a;named EGSY-2008532660). Its spectroscopic redshift was measured via Lyα emission via Keck/MOSFIRE to be z Ly = a 8.683 0.004 0.001 -+ from Zitrin et al. (2015;as EGSY8p7); at the time, and for several years, it was the farthest known Lyα emitter. It later gained a potential N V detection at z sys = 8.667, from Mainali et al. (2018), from additional Keck/MOSFIRE data, and it was also the brightest z > 8.5 galaxy found in the CANDELS survey (Finkelstein et al. 2022a;EGS_z910_6811). Here we report the results from JWST with the Cosmic Evolution Early Release Science (CEERS) Survey (Finkelstein et al. 2022b) data set. The NIRSpec spectroscopy (Arrabal Haro et al. 2023, in preparation) of this source has an ID from the micro shutter array (MSA) of 1019; thus, and hereafter, we refer to it as CEERS_1019.

In Section 2, we present the data from four different JWST observational modes from the CEERS Survey. In Section 3, we describe our method of emission-line fitting, while we explore the detected emission lines in detail in Section 4. In Section 5, we analyze the properties of this galaxy from the available imaging data and from our spectroscopic data in Section 6. We discuss the implications of these results in Section 7 and present our conclusions in Section 8. Throughout this paper we assume a flat Planck cosmology, with H 0 = 67.36 km s -1 Mpc -1 , Ω m = 0.3153, and Ω Λ = 0.6847 (Planck Collaboration et al. 2020). All magnitudes are in the AB system, and all restframe wavelengths are vacuum.

NIRSpec Data Reduction

The details of the CEERS NIRSpec data processing are presented in Arrabal Haro et al. (2023). We summarize the main steps of the reduction here. The NIRSpec data are processed with the Space Telescope Science Institute Calibration Pipelinefoot_1 (version 1.8.5;Bushouse et al. 2022a) and the Calibration Reference Data System (CRDS) context jwst_1027.pmap. We use the calwebb_detector1 pipeline module to subtract the bias and the dark current, correct the 1/f noise, and generate count-rate maps (CRMs) from the uncalibrated images. At this stage, the parameters of the jump step are modified for an improved correction of the "snowball" eventsfoot_2 associated with high-energy cosmic rays.

Table 1 Source Information for CEERS_1019 R.A. 215.0353914 [deg] (1) Decl. 52.8906618 [deg] (2) m F160W 25.2 ± 0.03 [AB mag] (3) m F356W 24.8 ± 0.01 [AB mag] (4) z phot 8.84 0.25 0.12 -+ (HST+IRAC) ( 5) z phot 8.72 0.06 0.04 -+ (+JWST) ( 6) z spec(Lyα) 8.683 0.004 0.001 -+ (Keck) ( 7) z spec(Lyα) 8.6854 ± 0.0045 (NIRSpec) ( 8) z spec( [O III]) 8.6788 ± 0.0002 (NIRSpec) ( 9) log(M å ) 9.5 ± 0.3 [M e ] (10) log(sSFR) -7.9 ± 0.3 [M e yr -1 ] (11) log(M BH ) 6.95 ± 0.37 [M e ] (12) T e ([O III]) 18,630.76 ± 3.68 [K] (13) 12+log(O/H) 7.66 ± 0.51 (14) Z 0.095 0.06 0.21 -+ [Z e ] (15) n e 1.9 ± 0.2 × 10 3 [cm -3 ] (16) UV Slope β 1.76 0.13 0.12 --+ (17) A v 0.4 ± 0.2 [mag] (18) Note. Column (1): right ascension. Column (2): declination. Column (3): observed AB magnitude in the F160W filter from Finkelstein et al. (2022a). Column (4): observed AB magnitude in the F356W filter (Section 2.2). Column (5): photometric redshift measured with HST (prior to JWST photometry; Finkelstein et al. 2022a). Column (6): photometric redshift measured including JWST photometry from CEERS (Section 5). Column (7): spectroscopic redshift measured from Keck/MOSFIRE spectroscopy via Lyα emission-line detection from Zitrin et al. (2015). Column (8): spectroscopic redshift measured from JWST/NIRSpec spectroscopy via Lyα. Column (9): spectroscopic redshift measured from JWST/NIRSpec spectroscopy via [O III] emission-line detection (Section 4.1). Column (10): stellar mass from Prospector SED fit (Section 5). Column (11): specific star formation rate (Section 5). Column (12): BH mass (Section 6.1). Column (13): electron temperature (Section 6.2). Column (14): T e -based metallicity (Section 6.2). Column (15): metallicity (Section 6). Column (16): electron density (Section 6.2). Column (17): UV spectral slope (Section 5.1). Column (18): dust attenuation (Sections 5 and 6.4).

The resulting CRMs are then processed with the cal-webb_spec2 pipeline module, which creates two-dimensional (2D) cutouts of the slitlets, performs the background subtraction making use of the three-nod pattern, corrects the flat fields, implements the wavelength and photometric calibrations, and resamples the 2D spectra to correct the distortion of the spectral trace. The pathloss step accounting for the slitloss correction is turned off at this stage of the reduction process (see Figure 2). Instead, we introduce slit-loss corrections based on the morphology of the sources in the NIRCam bands and the location of the slitlet hosting them.

The one-dimensional (1D) spectra of the sources were obtained via an optimal extraction (Horne 1986) with a spatial weight profile from the trace of the source, such that the pixels near the peak of the trace are maximally weighted. To create the extraction profile, the 2D signal-to-noise ratio (S/N) spectrum was collapsed in the spectral direction for each grating independently, taking the median value at each spectral pixel and fitting a Gaussian to the positive trace. This source is effectively unresolved at all wavelengths (the intrinsic size FWHM is less than the JWST point-spread function, or PSF, FWHM), so the central trace is unaffected by upper and lower negative traces. A significant trace was measured in both the G140M and G395M gratings, but was not as apparent in the G235M grating such that it could be fitted with a similar Gaussian profile. We thus used the profile from the G140M grating for the extraction from 2D to 1D in the G235M grating.

We evaluate the scaling factor for the slit-loss correction compared to the NIRCam photometry for this source, again using the reduction without the path-loss correction applied. In all NIRCam filters, we repeatedly calculate the fluxes enclosed within the 0 2 × 0 46 rectangle aperture, using the source and MSA shutter positions, and estimate the scaling factors to match them with the total flux measurements. The scaling factors range from approximately 2.0-2.5 among the NIRCam filters. We correct the spectrum with the scaling factor of the filter whose central wavelength is closest to that of the observed wavelength. This self-consistently applies an aperture correction accounting for the variable PSF across the observed wavelength range. We assess the uncertainty in these slit losses by comparing them to another method based solely on the F444W filter for multiple sources in the CEERS observations. We find differences in the correction factor on the order of 10%-40% (Fujimoto et al. 2023) on average for all sources. We note that the slit-loss correction is multiplicative, such that any systematic uncertainty in this correction does not affect the significance measurements for emission lines.

We test the accuracy of the error spectrum computed by the pipeline by comparing the normalized median absolute deviation of the science spectrum (masking out emission lines and removing a smoothed continuum) to the median of the error spectrum in each grating individually. The real data show fluctuations ∼1.5-2 times larger than the typical error value. Thus, we measure and scale the error spectrum up by this scale factor in each grating.

Once the errors are corrected, and the spectra are scaled to the NIRCam photometry, the three M gratings are combined into a single spectrum, resampling to a common wavelength array at the overlapping wavelengths and adopting the mean flux at each pixel, weighted by the flux errors. This combined spectrum is then used for the remainder of the analysis in this paper.

NIRCam Imaging

The galaxy discussed in this work was observed in the CEERS JWST/NIRCam (Rieke et al. 2003(Rieke et al. , 2005;;Beichman et al. 2012;Rieke et al. 2023) imaging taken in 2022 December (in CEERS NIRCam Field 8); these imaging data, including the detailed reduction steps, are described in Bagley et al. (2023). 44Photometry of this source was measured using Source Extractor (SE; Bertin & Arnouts 1996) in dual-image mode. The photometry procedure is broadly similar to that described in full by Finkelstein et al. (2022b), though with a few differences we detail here. First, to improve color accuracy, especially in the Hubble Space Telescope (HST)/WFC3 bands, we do a two-step PSF correction. For any filters with a PSF FWHM smaller than the NIRCam F277W (which includes HST/Hubble's Advanced Camera for Surveys (ACS), F606W and F814W and NIRCam F115W, F150W, and F200W), we convolve the images with a kernel designed to match these images' PSFs to that in F277W. For images with larger PSFs (WFC3 F105W, F125W, F140W, and F160W and NIRCam F356W, F410M, and F444W), we derive a correction factor. This is done by convolving the F277W image with a kernel designed to match the PSF to the larger PSF image, with a correction factor derived as the ratio of the flux in the broadened image to that in the original F277W image. In this way, we correct for light not captured in our default aperture in the bands with larger PSFs, without needing to smooth all bands to that with the largest PSF (which would be WFC3 F160W).

We tested this process using source injection simulations to confirm that accurate colors were recovered. Finally, following Finkelstein et al. (2022b), we derived accurate total fluxes, first by correcting the fluxes in all bands by an aperture correction derived in F277W (as the ratio between the flux in our small Kron aperture used to measure colors and in the default larger "FLUX_AUTO" Kron aperture), and then correcting this by a factor of 8%-10%, to account for missing light from the wings of the PSF, with this correction factor derived via the source injection simulations (see Finkelstein et al. 2022b for more details).

This final photometry catalog includes measurements over the full CEERS NIRCam wavelength range in the F115W, F150W, F200W, F277W, F356W, F410M, and F444W filters, which have exposures of ∼3000 s per filter (∼6000 s for F115W), as well as in the existing HST/CANDELS ACS and WFC3 F606W, F814W, F105W, F125W, F140W, and F160W bands. We include the F098M (nondetection) data for this source from Finkelstein et al. (2022a).

The measured photometric redshift for this galaxy before JWST was ( ) z HST 8.84 phot 0.25 0.12 = -+ (Finkelstein et al. 2022a). With the addition of the JWST/NIRCam imaging (and by removing the blended Spitzer/IRAC data), we measure a photometric redshift with EAZY (using the same process as Finkelstein et al. 2022b) nearly equal to the spectroscopic redshift (z spec = 8.679) of ( ) z HST JWST 8.68 phot 0.15 0.09 + = -+

Methods: Emission-line Search

To search for emission-line features in our 1D spectrum, we utilize an automated line-finding code first published in Larson et al. (2018) and outlined here. This code uses a Monte Carlo Markov Chain (MCMC) routine to fit a model that consists of a Gaussian line plus a continuum constant to a given wavelength range, with four free parameters: the continuum level, line central wavelength, line FWHM, and integrated line flux. To run the MCMC, we use an IDL implementation of the affineinvariant sampler (Goodman & Weare 2010) to sample the posterior similar to that used in Finkelstein et al. (2019), which is similar to the Python emcee package (Foreman-Mackey et al. 2013).

As a first pass-through for emission-line features, we search the entire spectrum with our automated code. At each wavelength pixel, we do an initial S/N check (flux/error). If this is greater than unity, we run our fitting routine (this "precheck" is not required, but helps with efficiency, as no detectable emission feature would have S/N < 1 at the line center). We use a fitting range of 100 times the pixel scale in the G140M grating on either side of the center pixel, which equates to 630 Å across the full spectrum. This ensures we fit the same wavelength range for every line, regardless of pixel scale. We fit a single Gaussian feature to the spectrum:

We also impose noninformative priors on each free parameter, designed to increase the efficiency of the line detections while removing the chance for the MCMC chain to exit a range of realistic parameters. For the continuum constant ( f c ), we let it vary between -3 to 10 times the average flux across the fitting range. These values are larger than the 1σ noise level and broadly encompass the typical continuum values for our source. We restrict the peak wavelength (λ 0 ) to be the wavelength at that pixel ± one pixel, such that we fit a Gaussian within each pixel. The actual value in Å varies, as the pixel scale (Δλ) for NIRSpec is wavelength-dependent, with Δλ G140M = 6.4 Å, Δλ G235M = 10.7 Å, and Δλ G395M = 17.9 Å. We limit the FWHM (FWHM ( ) 2 2 ln 2 2.355 s s = = to ± 30 km s -1 of the observed FWHM of the [O III] 5008 line, the brightest feature in the spectrum (see Figure 4 and Table 3). Tying the FWHM to [O III] 5008 in velocity units rather than pixels is critical, as the wavelength dispersion changes across the observed wavelength range. The line flux prior requires the line flux ( f 0 ) to be greater than the average of the flux over the plotting range after a >1σ clipping (estimate of the continuum at that location) and less than 100 times the average flux, well outside the expected range of line fluxes.

To define the starting point for our MCMC routine, we used the IDL routine mpfit, a Levenberg-Marquardt least-squares fitting routine that fits a Gaussian with the above parameters (Markwardt 2009). We then run our MCMC fitting code with 10,000 iterations and 100 walkers on each pixel and determine the best line fit results as outlined in the following sections. We use the median of the last 10,000 steps of our MCMC chain for our fit parameters. To measure the MCMC error on our parameters, we use the robust sigma calculation: using the median absolute deviation as the initial estimate, then weighting points using Tukey's biweight (Equation (9) from Beers et al. 1990), assuming a Gaussian posterior distribution. As our MCMC fits marginalize over the line flux and the continuum flux, the errors on both are included in our final line fit error.

We place several constraints on the resulting Gaussian fits when determining a successful line detection, including the wavelength and comparison to neighboring pixels. First, we mask out the edges of the wavelength range for the three M Gratings (within 0.97-5.1 μm), where the grating transmission curve falls below 70%. To determine the correct peak pixel for the line, we select the fit at the pixel that has the highest peak flux (the flux at the peak of the Gaussian fit) within the FWHM of the line. We calculate an integrated S/N from the MCMC Gaussian fit as the median line flux divided by the robust sigma of the line flux as described above. We also measure a peak S/N value, calculated by taking the maximum flux of the Gaussian fit and calculating the ratio of this to the noise per pixel within our fitting range (similar to Larson et al. 2022). For our initial search through the spectrum for emission lines, we only consider lines found at >5σ in both integrated and peak S/N. This initial pass is run for masking out significant features in the spectrum to fit and remove the detected continuum, as described in Section 3.1, and for identifying potentially unexpected emission-line detections in sources at this redshift. In Section 3.2, we discuss how these initial S/N measurements overestimated the significance of our emission-line detections and how we determine the final line flux errors and S/N values for our detected emission lines.

Determining the Significance of the Emission Lines

Upon inspection of the initially identified lines described above, it was found that this methodology tended to overestimate the S/N, especially for faint lines. We thus empirically derive line flux errors and appropriate S/N measurements for each line that passes the above selection cuts in the following way. Using the line-masked spectrum as above, we find the 20 closest "blank" pixels, free of any other nearby detected emission lines, on either side of the emission line (line center + FWHM, such that we are not overlapping the emission line at all). At each of these blank pixels, we insert a fake emission line with the same parameters as our detected emission line (F(λ), FHWM, and f c ) and a line center (λ 0 ) at the wavelength of that pixel. We then run the same line-fitting routine at each pixel and record the recovered line flux from the MCMC (F(λ) out ) for each of our 40 simulated lines. Our reported line flux error (ΔF(λ)) is then the median absolute deviation of these recovered line fluxes, as reported Figure 3. The spectrum of CEERS_1019 obtained with JWST/NIRCam's wide-field slitless spectroscopy mode. This consists of spectra with both the column ("C"; blue) and row ("R"; purple) grisms, both taken with the F356W blocking filter, with a combined spectrum shown in black. While these data have a higher noise level than the NIRSpec spectra, we observe detections of the same [O II] 3727+3729 , [Ne III] 3870 , and bluer Balmer lines that we see in NIRSpec (Figure 5). While we do not use the grism measurements in this paper, these data highlight the utility of this mode for obtaining slitless measurements of modestly faint emission lines out to high redshifts. . Plotted in green is the fit to the continuum in the spectrum (see Section 3.1) plus the detected emission lines, as discussed in Section 4. We require an S/N >2.3 (with our simulation-estimated noise as described in Section 3.2) for a line to be considered detected. [10

Lyα 1215 11775.10 ± 6.58 14.4 5.41 ± 0.38 209.77 ± 121.21 L L 1196.74 ± 470.98 10.67 ± 0.93 218.61 ± 332.94 N V 1238+1242 11969.40 ± 3.39 2.0 1.74 ± 0.88 358.93 ± 19.66 0.03 ± 0.03 1.70 ± 0.72 L 3.90 ± 2.30 -522.45 ± 85.20 N IV] 1483+1486 14357.60 ± 8.77 7.6 3.36 ± 0.44 363.99 ± 20.77 1.23 ± 0.85 2.13 ± 0.69 L 7.96 ± 1.47 16.15 ± 183.26 [O II] 3727+3729 36074.80 ± 4.09 11.9 2.94 ± 0.25 357.55 ± 20.83 1.80 ± 0.35 1.15 ± 0.35 L 34.79 ± 5.80 11.43 ± 34.71 [Ne III] 3869 37456.30 ± 2.29 12.84 3.5 ± 0.27 358.44 ± 17.43 L L L 44.27 ± 7.07 8.51 ± 19.64 He I 3889 +Hη 3890 37656.80 ± 9.73 2.3 0.66 ± 0.29 357.13 ± 18.12 L L L 8.43 ± 3.87 71.95 ± 77.78 [Ne III] 3968 +Hò 3971 38407.90 ± 3.75 5.5 2.04 ± 0.37 347.65 ± 21.01 0.48 ± 0.43 1.56 ± 0.40 L 26.53 ± 6.07 -23.34 ± 30.11 Hδ 4102 39717.90 ± 6.07 3.1 1.14 ± 0.37 358.56 ± 18.29 L L L 15.64 ± 5.50 53.32 ± 46.36 Hγ 4341 42018.20 ± 5.25 6.6 2.60 ± 0.39 360.56 ± 18.78 L L L 39.36 ± 8.48 -27.63 ± 38.12 [O III] 4364 42244.90 ± 9.91 3.0 1.16 ± 0.38 361.11 ± 17.79 L L L 17.78 ± 6.45 19.06 ± 70.68 Hβ 4862 47059.80 ± 4.91 10.2 7.32 ± 0.72 354.80 ± 19.61 3.79 ± 0.69 3.53 ± 1.51 1196.26 ± 349.07 128.94 ± 18.96 -31.55 ± 32.06 [O III] 4960 48009.60 ± 1.87 26.3 13.08 ± 0.50 358.52 ± 17.37 L L L 235.35 ± 27.51 1.00 ± 13.64 [O III] 5008 48473.50 ± 1.14 47.2 39.58 ± 0.84 358.66 ± 12.23 L L L 725.59 ± 82.2 -0.59 ± 124.41 C III] 1908 * 18493.30 ± 9.22 1.63 1.70 ± 1.05 642.90 ± 156.84 L L L 6.04 ± 3.85 309.51 ± 149.63 Mg II 2796+2803 27077.90 ± 6.67 1.18 1.62 ± 1.37 358.76 ± 20.82 0.67 ± 1.01 0.96 ± 1.20 L 11.56 ± 9.90 144.74 ± 74.18 Note. Identified emission-line values in CEERS_1019; see Figure 5 for individual plots. Column (1): line name and rest-frame vacuum wavelength in [Å]. Column (2): observed wavelength of the line in NIRSpec in [Å].

Column (3): S/N of the line (Section 3.2). Column (4): total integrated line flux in F λ ( units of 10 -18 erg s -1 cm -2 Å -1 ; Section 3). Column (5): FHWM of the line fit (in km s -1 ) for single-Gaussian fits, the blue line in the doublet (Section 3.4.2) and asymmetric (Section 3.4.1) fits, or the narrow component in the dual-component (Section 3.4.3) fits. Column (6): flux of the line fit (in 10 -18 erg s -1 cm -2 Å -1 ) for single-Gaussian fits, the blue line in the doublet (Section 3.4.2) and asymmetric (Section 3.4.1) fits, or the narrow component in the dual-component (Section 3.4.3) fits. Column (7): flux of the red line in the doublet (Section 3.4.2) fit or the broad component in the dual-component (Section 3.4.3) fit (in 10 -18 erg s -1 cm -2 Å -1 ). Column (8): FHWM of the red line in the asymmetric (Section 3.4.1) fit or the broad component in the dual-component (Section 3.4.3) fit (in km s -1 ). Column (9): rest-frame EW (Rest EW) of the line (in Å; Section 3.3). Column (10): velocity offset of the emission line from the systemic redshift of CEERS_1019 (z [O III] =8.679; Section 4.1; in km s -1 ). Lines detected at a low significance are reported below the horizontal line and are included for completeness. * While the C III] line is a doublet with peaks at 1906 and 1908 Å in the rest frame, we do not resolve this doublet and show the best fit, which uses a single-Gaussian profile with a broad FWHWM (Section 4).

in column 4 of Table 3. The S/N of our emission lines is taken as the measured line flux from our line fit divided by this error, S/N = F(λ)/Δf (λ), and is reported in column 3 of Table 3. To derive a consistent threshold in S/N where we consider a line to be detected, we examine the results of the automated line-finding routine, searching for lines with restframe wavelengths that do not match a known line (using a large line list). The point below which we detect no "unknown" (spurious) emission lines is S/N = 2.3, which we use as our detection threshold for emission lines in CEERS_1019, unless otherwise noted below.

Asymmetric Gaussian Fit

Given the detection of Lyα in our NIRSpec spectra, we want to fit an appropriate asymmetric Gaussian profile to the data to account for intergalactic medium (IGM) absorption on the blue side of the line. For this fit, we add an extra FWHM parameter, such that we now have "blue" and "red" FWHMs of the line. The blue FWHM can vary between 0, as the blue side of the line may be fully attenuated, and FWHM [O III] + 30 km s -1 , which is the same maximum as the single-Gaussian fit. We start the MCMC chain for this parameter at FWHM

[O III] . The red FWHM is less restricted, as it can vary between FWHM [O III] and 2000 km s -1 , and the MCMC chain is started at FWHM [O III] + 30 km s -1 , the maximum value for the blue FWHM. Following the same method as Jung et al. (2020), the equation for the asymmetric Gaussian is:

This returns the continuum flux ( f c ), peak line flux value ( f 0 ), a peak wavelength (λ 0 ), and the blue-and red-side line widths (σ blue and σ red ). To measure our line flux, we integrate this emission-line profile and report the integrated line flux in column 4 of Table 3 3.4.2. Doublet Gaussian Fits Several of the UV and optical emission lines are not a single emission-line feature, but rather a doublet. To accurately fit these lines, we implement a doublet Gaussian line profile to our MCMC routine:

For this fit, we employ the continuum flux ( f c ) plus two Gaussians with the same FWHM (σ) and individual line fluxes for each line ( f 1 and f 2 , respectively), with the blue peak (λ 1 ) fixed to ±6 pixels from the expected location and the red peak fixed to the redshifted separation between the two doublet lines (δλ = δλ rest (1 + z)). This returns a total combined line flux for the doublet ( f (λ)) as reported in column 4 in Table 3, as well as individual line fluxes for each component, as reported in columns 6 and 7 of that same table.

Broad Hβ Feature

The observed profile of the Hβ emission line shows a strong narrow component, similar to [O III], but also a weaker broad component. To explore the significance of this broad component, we perform a dual Gaussian fit to this line, as described in Section 3.4.3. Uncertainties on both components were derived via similar simulations, as described in Section 3.2. From this fit, we find a distinct broad component to the Hβ emission line, with an FWHM = 1196 ± 349 km s -1 . We show this fit (green) in Figure 6. The flux in this broad component is comparable to the flux in the narrow component, with a broadto-narrow flux ratio of 0.93.

We consider multiple measures of the significance of this broad feature. The first is the S/N of the broadline flux from this fit, which is 2.5. However, this does not consider how good the fit is without the broadline component. We next try a narrowline-only fit, with the FWHM allowed to be free (i.e., not tied to [O III]), shown as the blue line in Figure 6. We find a single-line FWHM of 500 ± 52 km s -1 , significantly broader than the fiducial narrowline-only FWHM of 367 ± 17 km s -1 (the purple line in Figure 6), making it clear that the broad component affects even a single-line fit.

Finally, we calculate theBayesian Information Criterion (BIC) between the dual-component and narrowline-only fits. The BIC is a method of comparing the goodness-of-fit between two models, accounting for differing degrees of freedom. It is defined as BIC= ( )

, where k is the number of free parameters and N is the number of data points (the number of spectral pixels, in this case; Liddle 2004). Our narrowline-only fit with four free parameters ( f c , f 0 , λ 0 , and σ-see Section 3) has a BIC of 21.4 when the FWHM is tied to [O III] and 21.8 when the FWHM is free. The narrow+broad fit with six free parameters ( f c , λ 0 , f nar , f broad , σ nar , and σ broad -see Section 3.4.3) has a BIC of 18.7. The ΔBIC=2.7-3.1 gives "positive" evidence (using the scale of Jeffreys 1961, where ΔBIC > 6 is "strong" evidence) that a broadline component is necessary to fit the Hβ line successfully.

This significant broad component of Hβ thus indicates the presence of high-velocity gas, which we interpret as emitting from the broadline region of an AGN. Large-scale outflows can similarly produce broad velocity widths, but would be apparent in all emission lines (e.g., Amorín et al. 2012;Hogarth et al. 2020). We do not observe any broad features in lines with much higher S/N, like [O III]λ5008. Instead, [O III] and other forbidden lines have widths that are inconsistent with the outflow scenario and are instead consistent with the narrow Hβ component. Broadline AGNs typically exhibit broad components for permitted lines (like Hβ) and narrow widths of forbidden lines (like [O III]; e.g., Schmidt et al. 2016;Vanden Berk et al. 2001), as observed in our spectrum. We discuss further evidence for the AGN nature of this source in Section 4.4 below.

Lya

The Lyα line was first detected by Zitrin et al. (2015) at z Lyα = 8.683 with Keck/MOSFIRE. We detect an emission line at this same wavelength, which we fit with an asymmetric Gaussian profile, as described in Section 3.4.1, and measure a line flux f Lyα = 5.41 ± 0.38 × 10 -18 erg s -1 cm -2 . Using our measured central wavelength and the vacuum wavelength for Lyα of 1215.67 Å, we measure a Lyα-based redshift of 8.6854 ± 0.0045, 218 ± 333 km s -1 higher than that of [O III]. This is lower, although consistent within ∼1σ, of the measured Lyα velocity offset (relative to N V) of +362 km s -1 from Mainali et al. (2018).

We find a clear asymmetric line profile with an extended redside tail (FWHM red > 1000 km s -1 ). Although it becomes more difficult to detect Lyα at this high redshift, due to the increasingly neutral IGM, one can expect the escape of Lyα that is scattered and redshifted enough to avoid resonant scattering (Dijkstra et al. 2004), as seen in the extended red tail in our Lyα spectra. Additionally, a sharp cutoff at the blue edge of the line profile may indicate a significant contribution of resonant scattering, which can be caused by a proximate optically thick medium (Mason & Gronke 2020). Particularly, reionization simulations predict the sharp blue-side edge of the red peak, due to the infall motion of neutral gas around a galaxy, as seen in our Lyα spectrum. The cutoff location is found at ∼200 km s -1 in our spectra, which is comparable to the predicted velocity of the infalling gas around an M UV = -22 galaxy at this redshift (Park et al. 2021).

Our Lyα line flux is ∼seven times fainter than the line flux measured from MOSFIRE. This could be due, in part, to the MOSFIRE slits being wider than the NIRSpec micro shutters (0 7 versus 0 2). Our NIRSpec slit-loss correction would not correct for extended Lyα emission in this source. The MOSFIRE observations of this source were taken ∼8-9 yr prior to the JWST/NIRSpec data, ∼1 yr in the rest frame, and it is possible that the variability common in AGNs might contribute to the difference in the measured line fluxes. AGN variability is typically ∼20% (e.g., MacLeod et al. 2012) and is unlikely to fully explain this large difference.

The N V 1243 line was detected from this source with Keck/ MOSFIRE by Mainali et al. (2018) at 12019.5 Å, with the other line in the doublet, N V 1238 , obscured by a skyline at 11981 Å.

In our JWST/NIRSpec data, we do not detect a line at either of these wavelengths, but there is a 1.8σ peak at 12007 Å. This would correspond to a rest-frame wavelength of 1240 Å between the peaks of the N V doublet. The velocity offset of the N V line from the systemic redshift z [O III] = 8.679, if indeed the λ1238 line is ∼ + 442 ± 85 km s -1 . This is consistent (within the error) with the offset we measure for the Figure 5. Fits to the emission lines identified in the JWST/NIRSpec spectrum of CEERS_1019. Each panel shows an individual emission-line fit, with the type of line profile as described in Section 3. The panels are plotted in F λ (10 -19 erg s -1 cm -2 Å -1 ) vs. observed wavelength (μm), are presented in order of increasing wavelength, and are scaled vertically to show the extent of the highlighted emission line. Emission-line values are tabulated in Table 3 and discussions of specific lines are given in Section 4.

C III] of ∼ + 310 ± 150 km s -1 . If it is the λ1243 line, the velocity offset is ∼ -521 ± 85 km s -1 .

Other Potential Broad Emission Lines

As noted above, broadline (Type 1) AGNs typically exhibit broad permitted lines and narrow forbidden lines, interpreted as high-velocity dense gas in a "broadline region" near the BH (and dominated by its gravitational influence) and lowervelocity gas in a "narrowline region" that is more distant from the BH and dominated by the galaxy kinematics (e.g., Netzer 2015). Unfortunately, the Hα line is redshifted beyond the wavelength range of NIRSpec for our target (though it is observable by MIRI spectroscopy). The other lines in the Balmer series (Hγ, Hδ, etc.) are too weak for strong constraints on the presence of broad components (see Figure 6). In AGNs with a broad Hβ line, broad features are also typically observed in the permitted UV lines Lyα, C IV, C III], and Mg II with similar or stronger fluxes to the broad Hβ line (e.g., Vanden Berk et al. 2001).

We detect a significant broad width in the observed Lyα line of our target. The blue-side absorption of the Lyα forest makes it difficult to measure the line width precisely, but we note that our asymmetric Gaussian fit (described in Section 3.3) finds a red-side width of FWHM = 1360.6 ± 479.9 km s -1 : comparable to the best-fit width of the broad Hβ component.

The best-fit Gaussian to the C III] line also indicates a broadline width of FWHM = 643 ± 157 km s -1 . The observed C III] feature is not well fit by two narrow Gaussians, with widths constrained to be the same as the [O III]λ5008 line for each of the doublet lines. This could indicate that the doublet is unresolved in a noisy part of our data, leading to a single line with a broader width. Alternatively, the broader FHWM of this line compared to [O III] may have some physical implications for this galaxy. The line width of the best-fit broad C III] profile is somewhat narrower than would be expected for a stratified broadline region, in which the higher-ionization C III] gas orbits closer to the BH and consequently has a broader velocity width than Hβ. The line center of the best-fit broad C III] profile is redshifted by +310 km s -1 with respect to the systemic redshift (determined from the narrow [O III]λ5008 line as described in Section 4.1). The shift in line center and narrower-thanexpected width may indicate a blueshifted absorption component in the C III] line, as commonly observed in broad absorption-line quasars (e.g., Gibson et al. 2009).

We additionally see a significant broad component in the best-fit profile for the N IV]λ1486 line. This feature is not typically observed in luminous quasars at z  4, but it is detected as a feature in the NIRSpec prism spectrum of the z = 10.6 galaxy GN-z11 (Bunker et al. 2023). Its presence in GN-z11 has been interpreted as evidence for high-density and nitrogen-enhanced gas (Cameron et al. 2023;Senchyna et al. 2023), both properties consistent with expectations for broadline regions around rapidly accreting and low-mass AGNs (Hamann & Ferland 1999;Matsuoka et al. 2017).

On the other hand, we do not observe broad components in the C IV and Mg II emission lines. The absence of C IV can be explained by the low S/N in its region of the NIRSpec spectrum, and the large flux uncertainty at that wavelength range allows for an undetected broad C IV line that is of comparable brightness to the broad Hβ feature. The Mg II feature is more puzzling: the observed Mg II profile is best fit by narrow Gaussians for each of the doublet lines, consistent with galaxy (or narrowline region) emission, rather than the expected broad feature of an AGN. It is possible that this line is intrinsically broad but is affected by intervening absorption by the Mg II doublet: this requires a particular blueshift velocity of the absorption component, but is not uncommon in quasars. Alternatively, the AGN may be attenuated by dust, such that its UV emission lines are weaker than the Balmer lines, as implied by the best-fit SED model, although this conflicts with the detection of broad Lyα, N IV], and C III]. Deeper spectra of this source would more effectively test for the presence of the broad UV lines expected for AGNs, especially in investigating the hypothesis for blueshifted absorption affecting the broad C III] line and the Mg II line. Alternative single-Gaussian line profile fits are shown where the FWHM is restricted to FWHM [O III] ± 30 km s -1 (blue) or left unrestricted (purple). In both cases, the fit to the Hβ emission line is worse than the dual-component fits, as indicated by the higher Bayesian Information Criterion (BIC) measurement (Section 4.2). Right: scaled-down fits of the same dual-component fit from Hβ to the peak of the Hδ (top) and Hγ (bottom) emission lines. In both cases, the broad component is not distinguishable from the noise, indicating that a broad component could be present in these lines below our noise threshold.

Emission-line Flux Ratio Measurements

Despite potential lingering issues in the absolute flux calibration of the NIRSpec instrument, we find that the relative flux calibration is consistent, at least for pairs of emission lines near one another in wavelength. We measure [O III] 5008 / [O III] 4960 = 3.03 ± 0.13, which is consistent with the atomic physics calculation within errors (Storey & Zeippen 2000). We are confident in using the ratios of emission lines close in wavelength. Still, we acknowledge that additional flux calibration and improvements of the NIRSpec instrument are needed to trust widely separated line ratios fully. We also note that our wavelength-dependent aperture correction using the differing spatial profile of CEERS_1019 across the three NIRSpec filters may add additional systematic errors in distant line ratios. We thus only report ratios for lines that fall within the same NIRSpec filter: G395M. Since we also use the NIRCam filters to perform a wavelength-dependent flux correction (see Section 2.1.1), the errors on the line ratios that span across the two redder NIRCam filters (F356W and F444W) may be underrepresented. We report the following line ratios in this paper given the above caveats.

We list the measured line ratios for this source in Table 4. All reported line ratios are measured without using the broadline fit to the Hβ line, and only the narrow component from the dual-component fit (See Figure 6). For any of the line ratios reported where one line is not detected at a significant level, we report the 1σ upper limit for the line and mark it with an asterisk in the tables of line fluxes and ratios (Tables 3 and 4).

Constraints from Continuum Emission

We explore fitting the photometric SED of this object with stellar and AGN emissions simultaneously to explore which dominates the observed continuum emission. We first use the FAST v1.1 (Kriek et al. 2009;Aird et al. 2018) SED fitting code, including a star-forming galaxy and an AGN component, as described in Kocevski et al. (2023). With this code, we try fits using two sets of AGN templates. The first uses eight empirically determined AGN templates. These include five AGN-dominated templates from the Polletta et al. (2007) SWIRE template library (namely, the Torus, TQSO1, BQSO1, QSO1, and QSO2 templates) and three composite SEDs of X-ray-selected AGNs with absorption column densities of N H = 10 22-23 , 10 23-24 , and 10 24-25 cm -2 from Silva et al. (2004). See Appendix A of Aird et al. (2018) for additional details. We also try a second fit with a low-metallicity (Z 0.4 Z e ) AGN component from CLOUDY modeling, chosen to mimic the SED of a radio-quiet AGN (see Section 6.5 for more details). We measure the ratio of light from the best-fitting AGN model to the total stellar+AGN model in two wavelength windows: a rest-UV window at 0.15-0.25 μm and a restoptical window at 0.51-0.6 μm (both designed to avoid strong emission lines). In both of these fits, the photometry constrains the AGN to be subdominant to the stellar emission. With the fiducial AGN template, the stellar emission comprises 99.5% of the UV flux and 85% of the optical flux. With the lowmetallicity AGN model, stellar emission comprises 82%-83% of both the UV and optical flux.

We try an alternative fit using the Cigale code (Boquien et al. 2019;Yang et al. 2020Yang et al. , 2022)). In the fit, we adopt a modified Schartmann et al. (2005) AGN accretion disk model. We allow the deviation from the default optical spectral slope (the δ AGN parameter) to vary from -1 to 1 and the polar dust extinction (the E(B -V ) PD parameter) to vary from 0 to 1 (see Section 4 of Yang et al. 2022 for details). Measuring the results in the same two wavelength windows as above shows that the stellar emission comprises 70% and 59% of the continuum emission in the rest-UV and rest-optical, respectively.

These results imply that the observed continuum emission is dominated by stellar light. We thus proceed to perform Bayesian SED modeling with stellar-only models to explore what constraints can be placed on the stellar population of the host galaxy. Our fiducial results come from the Prospector SED fitting code (Johnson et al. 2021), following Tacchella et al. (2022) by using both a continuity and bursty prior on the star formation history (SFH). Both SFH priors lead to similar posterior distributions, implying that the data are informative and the impact of the prior is minimal. We find that this galaxy has a stellar mass of M å ≈ 10 9.5±0.3 M e and is actively forming stars (specific star formation rate, or sSFR ≈ 10 -7.9±0.3 yr -1 ), doubling its mass every ∼100 Myr. We find a mass-weighted age of t 34 Myr

50 29 119 = - +

. We also note that using the conversion from the star formation rate (SFR) to Balmer line flux from Kennicutt & Evans (2012), the Hβ flux estimated from this SFR (1.3 -5.3 × 10 -18 erg s -1 cm -2 ) is fully consistent with our measured narrowline Hβ flux (3.8 ± 0.7 × 10 -18 erg s -1 cm -2 ) within errors. We list the physical properties of this galaxy from our SED fitting in Table 5 and show the model fits in Figure 7.

We also fit these data with BAGPIPES (Carnall et al. 2018). We generally followed the procedures in Papovich et al. (2022), but we have now used Binary Population and Spectral

Table 4 Emission-line Ratios for CEERS_1019 Emission-line Ratio Measured Value Convention Oxygen [O III] 5008 / [O III] 4960 3.03 ± 0.13 [O III] 4363 / [O III] 5008+4960 0.02 ± 0.01 RO3 [O III] 5008 / [O II] 3727+3729 13.46 ± 1.18 O32 O/H [O III] 5008+4960 + [O II] 3727+3729 / Hβ 14.67 ± 2.68 R23 [O II] 3727+3729 / Hβ 0.78 ± 0.16 [O II] 3727+3729 / Hγ 1.13 ± 1.14 [O III] 4363 / Hγ 0.45 ± 0.16 [O III] 5008 / Hβ 10.44 ± 1.91 R3 Ne [Ne III] 3869 / [O II] 3727+3729 1.18 ± 0.14 [Ne III] 3869 / Hγ 1.34 ± 0.22 Balmer Series Hγ / Hβ 0.69 ± 0.16 0.36 ± 0.06 Total Hβ Hβ / Hγ 1.46 ± 0.34 2.82 ± 0.49 Total Hβ Hβ / Hδ 3.33 ± 1.24 6.42 ± 2.18 Total Hβ Hβ / Hò 2.43 ± 0.76 4.69 ± 1.29 Total Hβ Note. Common emission-line flux ratios as measured from JWST/NIRSpec for CEERS_1019. All of the ratios that include Hβ use the narrow component flux from the dual-component Gaussian fit (Section 3.4.3), as reported in column 6 of Table 4, unless otherwise noted as using the Total Hβ, which is the value in column 4. Synthesis (BPASS) v2.2.1 stellar population models (Eldridge et al. 2017), an SFH represented as a Gaussian mixture model (Iyer et al. 2019), and fit over a range of ionization parameters log U = [ -4, -1] to model the nebular emission. From these fits, the galaxy has a stellar mass of  M M log 9.3 0.1 =  *

, consistent with the fits from Prospector, though a slightly higher log(sSFR) of -7.5 ± 0.2 yr -1 . The systematic uncertainties on stellar mass and SFR here are ≈0.2 dex. The inferred properties, such as stellar mass, SFR, and dust attenuation, are generally consistent with and within the systematic uncertainties of the Prospector fits, but require a slightly higher SFR to account for differences in modeling assumptions. We show the best-fitting Prospector and BAG-PIPES models in Figure 7, alongside images of this galaxy in all observed filters and photometric redshift results from EAZY.

Morphology

We investigate the morphology of this source using the fitting codes Galfit (Peng et al. 2002(Peng et al. , 2010) ) and statmorph (Rodriguez-Gomez et al. 2019). We use the Galfit least-squares fitting algorithm to fit the galaxyʼs surface brightness profile. We use a 100 × 100 pixel cutout of the F200W science image as input, the corresponding error array (the "ERR" extension) as the input sigma image, and the empirically derived PSFs. As an initial guess, we use the source location, magnitude, size, position angle, and axis ratios from the Source Extractor (SE) catalog. The Sérsic index is allowed to vary between 0.01 and 8, the magnitude of the galaxy between 0 and 45, the effective radius (R e ) between 0.3 and 100 pixels, and the axis ratio between 0.0001 and 1. We also allow Galfit to oversample the PSF by a factor 9. We then visually inspect the best-fit model and image residual for each source to ensure that the fits are reasonable and that minimal flux remained in the residual.

We find that a single Sérsic component poorly fits the source. Visually, the source is extended asymmetrically, with three distinct components (see Figure 8). The optimum fit is obtained when the central region is fit with both a point source and a Sérsic component, and the two other regions to the west and the northeast are each fitted with their own Sérsic components. The requirement of a point-source component for the fit is consistent with the source having an AGN. We similarly performed fits on the F356W and F444W images and again found that single-Sérsic fits do not work because of the multiple components, even with the lower resolution at these longer wavelengths. The asymmetric nature of the source and the presence of the three separate components are consistent with the galaxy being involved in a major merger. We use statmorph to measure the concentration parameter and size as a function of wavelength on PSF-matched cutouts and find that C = 2.82 and R e = 6.49 pixels (0.91 kpc) for F356W and C = 2.93 and R e = 5.64 pixels (0.79 kpc) for F444W, suggesting that the F444W light is more concentrated than the F356W or F200W emission, supporting the visual impression in Figure 2. At the redshift of this galaxy, the F444W emission is dominated by the [O III] emission line (see Figure 7).

Black HoleMass

In this section, we estimate the virial mass of the central BH, assuming that the observed broad Hβ emission traces the kinematics of gas in the broadline region. The measurement of BH masses with single-epoch spectra can be made using the width of the broad Hβ emission line and the rest-frame 5100 Å continuum luminosity, L 5100 , which has been shown to correlate with the distance to the broadline region using reverberation mapping (e.g., Kaspi et al. 2000;Cackett et al. 2021). However, this assumes the rest-frame 5100 Å continuum luminosity is dominated by light from the AGN, which is likely not the case here (see Section 6.1). As a result, we instead use Equation (10) from Greene & Ho (2005), which

Table 5 Physical Properties of CEERS_1019 SFH Prior M UV,obs UV Slope β log M å log SFR 50 log sSFR 50

Continuous 22.07 0.05 0.05 --+ 1.77 0.12 0.11 --+ 9.4 0.2 0.3 -+ 1.5 0.2 0.2 -+ 7.9 0.4 0.2 -- + 0.4 0.2 0.2 -+ 0.3 1.0 0.3 -- + Bursty 22.09 0.05 0.06 --+ 1.80 0.10 0.09 --+ 9.5 0.4 0.2 -+ 1.5 0.1 0.4 -+ 7.9 0.4 0.2 -- + 0.4 0.2 0.4 -+ 0.6 1.1 0.6 -- +

Note. Physical properties of CEERS_1019 as derived through SED fitting to the HST and JWST photometry at the spectroscopic redshift. SED fitting with Prospector (Johnson et al. 2021) using both a continuity and bursty SFH prior, as described in Tacchella et al. (2022). Column (1); SFH prior used in the fit. Column (2): absolute UV magnitude at rest-frame 1500 Å. Column (3): UV spectral slope, b , as measured from the model spectra. Column (4): stellar mass. Column (5): SFR averaged over the past 50 Myr. Column (6): specific SFR averaged over the past 50 Myr. Column (7): dust attenuation at 5500 Å, from SED models. Column (8): stellar metallicity from SED models to the HST photometry. The values in columns (4) to (8) are calculated from SED fitting as in Tacchella et al. (2022).

employs only the Hβ line luminosity, L Hβ , and width:

4 10 10 erg s FWHM 10 km s . 5 BH 6 H 42 1 0.59 H 3 1 2 = ´b b --

This equation is based on the formula of Kaspi et al. (2000), but with L Hβ substituted for L 5100 using the empirical correlation between Balmer emission-line luminosities and L 5100 reported in Greene & Ho (2005).

Using the luminosity and width of the broad Hβ component, we derive a BH mass of ( )

. This mass is 1-2 dex lower than existing samples of luminous quasars with BH mass estimates at z > 5 and more comparable to the low-luminosity AGN reported in Kocevski et al. (2023), found in these same CEERS NIRSpec data. Our measured mass implies that the AGN in CEERS_1019 is powered by the least massive BH known in the Universe during the EoR.

To determine the accretion rate onto the BH relative to the Eddington limit, we derive a bolometric luminosity, L Bol , from the broad Hβ line luminosity and compare it to the Eddington luminosity, L Edd , for our measured BH mass. Assuming an intrinsic broadline Hα/Hβ ratio of 3.06 (Dong et al. 2008) and The source is a clear dropout in the optical and is well detected at 1-8 μm. The large inset plot shows the photometric SED from HST and JWST. We also show Spitzer/IRAC measurements from Finkelstein et al. (2022a) with both TPHOT (light gray) and Galfit (dark gray); this source was highly blended in IRAC, and the Galfit measurements appear closer to the NIRCam measurements. The lines are best-fit models from Prospector (gray) and BAGPIPES (purple). The SED is dominated by stellar emission (Section 5), and these stellar models infer a massive (log M/M e ∼ 9.5) and highly star-forming (log sSFR ∼ -7.9 Gyr -1 ) stellar population. The small inset panel shows constraints on the photometric redshift before (green) and after (red and yellow) the inclusion of JWST data. All three are consistent with the spectroscopic redshift, with those including JWST data placing much tighter constraints. a bolometric correction of L Bol = 130 × L Hα (Richards et al. 2006;Stern & Laor 2012), we obtain L Bol = 1.4 ± 0.6 × 10 45 erg s -1 . This results in an Eddington ratio, λ Edd = L Bol /L Edd , of 1.3 ± 0.5. This suggests that the BH is undergoing a rapid growth phase and is accreting at approximately its Eddington limit.

We note that while we observe a broad Hβ line, we do not detect AGN continuum emission at the level that might be expected given the Greene & Ho (2005) relationship between L Hβ and L 5100 . This may indicate the AGN is moderately reddened due to a patchy obscuring medium, similar to what is found in red quasars, which show broad Balmer lines, yet have properties intermediate between type I and II quasars (Greene et al. 2014;Glikman et al. 2023). If this is the case and the Hβ emission is suppressed by dust extinction, then our estimated BH mass should be considered a lower limit to the true virial mass for this AGN. We report our measurements fo the BH in CEERS_1019 in Table 6.

A v from Balmer Decrement

The ratios of Balmer lines can be used to measure nebular attenuation. In particular, gas with Case B recombination and a temperature of T e = 20,000 K (as calculated for our object in Section 6.2) will have intrinsic Balmer ratios of Hβ/Hγ = 2.110 and Hβ/Hδ = 3.817 (Osterbrock 1989). We can estimate the dust attenuation E(B-V) with an assumed attenuation curve (here, we use Calzetti 2001), by comparing our observed line ratios to these intrinsic values. For the observed line ratios, we use the total line flux rather than the narrowline-only flux, because the Hγ and Hδ lines do not have sufficient S/N to separate broad and narrow components. This means that our estimated nebular attenuation may be a mix of attenuation affecting the AGN and the narrow lines.

From the observed line ratios of Hβ/Hγ = 2.82 ± 0.49 and Hβ/Hδ = 6.42 ± 2.18, we find nebular attenuation of:

1. E(B -V ) Hβ/Hγ = 0.64 ± 0.41; 2. E(B -V ) Hβ/Hδ = 0.83 ± 0.54.

The observed line ratios are statistically consistent with the intrinsic values, such that the implied E(B -V ) values are

Table 6 CEERS_1019 BH Properties log M BH 6.95 ± 0.37 [M e ] (1) FWHM Hβ, Broad 1196.26 ± 349.07 similarly consistent with little to no dust. But the large error bars also allow for a broad range of nebular attenuation. In Section 5.1, the Prospector SED fit implied A V = 0.4 ± 0.2, which equates to E(B -V ) = 0.1 ± 0.05. This is consistent within the large error bars of the Balmer decrement estimate for dust attenuation.

Ionization Parameter

Tofoot_8 place this source in the context of other star-forming galaxies and AGNs identified in this early epoch, we investigate the "OHNO" line ratio diagnostic diagram. This diagram, comparing [O III] 5008 /Hβ versus [Ne III] 3870 / [O II] 3727+3729 , has been used at lower redshifts to identify ionization powered purely by star formation from that of an AGN (e.g., Backhaus et al. 2022;Cleri et al. 2023Cleri et al. , 2023)). Recent studies of high-redshift galaxies (z > 5) using JWST spectroscopy have found most if not all of their sources occupying the AGN region of the diagram (i.e., above the theoretical line dividing star formation from AGN, calibrated via low-redshift observations)-results that seem to suggest strong ionization from AGNs and/or very-metal-poor starforming H II regions at these high redshifts (e.g., Trump et al. 2023;Kocevski et al. 2023;Tang et al. 2023;Übler et al. 2023).

Such is the case for CEERS_1019, with high measurements of the ionization indicators [O III] 5008 /Hβ ≡ R3 = 10.44 ± 1.91 and [O III] 4960+5008 /[O II] 3727+3729 ≡ O32 = 13.46 ± 1.18 (see Table 4). Figure 9 shows our galaxy on the "OHNO" diagram (orange stars), using the narrow component of the Hβ flux (the bigger star-also the light green line in Figure 6; the smaller star shows the ratio with the full Hβ flux). Additional sources from recent JWST studies focused on the EoR (Trump et al. 2023;Kocevski et al. 2023;Tang et al. 2023) are included as small gray polygons. The black curve in both panels denotes the boundary between star-forming and AGN regions as defined by Backhaus et al. (2022), with the associated gray shaded region showing the z ∼ 1 "star-forming" region. Underplotted are photoionization models from both the CLOUDY v17 (Ferland et al. 2017; left panel) and MAPPINGS V (Sutherland et al. 2018; right panel) codes, with parameters chosen to showcase a range of possible galaxy properties at this redshift, including elemental abundances, stellar populations, and stellar and gas-phase metallicities. In this figure, each line in the left panel represents CLOUDY models with a specific ionization parameter (increasing log 10 U, shown by the darker colors), where the sizes of the points in each line represent the specific gas-phase metallicity of said model. Similarly, each line in the right panel represents the MAPPINGS models for a specific ionization parameter (similar to the left panel) and isobaric pressure (denoted by the line types).

The CLOUDY models shown in Figure 9 assume a planeparallel geometry with a nebular electron density of 10 3 cm -3 (inferred from the flux ratio of the [O II] 3727+3729 doublet; see Section 6.3), scaled Solar elemental abundances, and cover ionization parameters from log 10 U = [-2.5, -2.1, -1.5]. The star-forming models use BPASS v2.0 (Eldridge et al. 2017) fiducial binary stellar population models with an IMF that extends to 300 M e (described by a slope of α = -1.30 from 0.1 to 0.5 M e and α = -2.35 from 0.5 to 300 M e ) and continuous star formation of 1 M e yr -1 . For these models, we fix the stellar and gas-phase metallicities, covering Z=0.1-0.5 Z e . The AGN models use the table agn model in CLOUDY, which approximates a "typical" radio-quiet AGN, covering the same range of ionization parameters as the star-forming models and spanning nebular metallicities of 0.05-0.5 Z e .

The MAPPINGS models shown in the same figure are the "Pressure Models" from Kewley et al. (2019a), assuming a (Trump et al. 2023;Kocevski et al. 2023;Tang et al. 2023). Underplotted on both panels are photoionization models, with the left panel showing star-forming (blue) and AGN models (red) using the CLOUDY v17 code (Ferland et al. 2017) and the right panel showing star-forming models (green) from Kewley et al. (2019a) using the MAPPINGS V code (Sutherland et al. 2018). Each line of points represents a specific set of ionization parameters of log 10 U = [-2.5, -2.1, -1.5] for the CLOUDY models and log 10 U = [-2.5, -1.5] for the MAPPINGS models, with increasing ionization toward darker colors (or toward the top right of each panel). For the MAPPINGS models, the lines are further separated by isobaric pressure (denoted by line type). The sizes of the points represent the metallicities of the models.

plane-parallel geometry with a range of pressure (log 10 (P/k)=7-9 cm -3 ), ionization parameter (log 10 U = [-2.5, -1.5]; also referred to as log 10 Q = 8-9), and metallicity (Z = 0.05 -1 Z e ). These star-forming models use the Star-burst99 (Leitherer et al. 2014) stellar population models with a Salpeter IMF (Salpeter 1955) extending to 100 M e . The stellar and gas-phase metallicities follow a prescription such that at lower metallicities, the models have increasingly enhanced α abundances (see Nicholls et al. 2017 for more details). We have included these models in addition to the CLOUDY models (which have Solar abundances scaled to chosen metallicity) to showcase the range of potential galaxy properties for our highredshift galaxy.

Similar to that found in other spectroscopic studies in this epoch (e.g., Trump et al. 2023;Kocevski et al. 2023;Tang et al. 2023), the line ratio when using the single-component fit to the Hβ line for CEERS_1019 (the smaller star) lies in a region on this diagnostic that is difficult to differentiate between AGN and metal-poor high-ionization star-forming H II regions. This is not unexpected, as galaxies at higher redshifts have been shown to generally have higher ionization and lower metallicities in comparison to those at lower redshifts (e.g., Shapley et al. 2003;Erb et al. 2010;Kewley et al. 2019a;Backhaus et al. 2022;Papovich et al. 2022;Sanders et al. 2023). This suggests that at these high redshifts, successfully separating star-forming and AGN sources using these strong line diagnostics can be challenging-evidenced by the lower-metallicity AGN models overlapping with the starforming models in the left panel of Figure 9. Further discussion about the general utility of such emission-line diagnostics at high redshift is discussed in Section 7.3.

However, when focusing on the fit to only the narrow component of the Hβ line for CEERS_1019 (the bigger star), there is a clear distinction from the rest of the high-redshift sources shown in this figure. Of the sources shown, our galaxy is the only one with enough S/N in Hβ to (a) see a broadening of the line profile, and (b) clearly measure both a broad and narrow component for the line. 49 Follow-up spectroscopy of these and other sources that fall in this strong line regime at high redshift may shed light on more "hidden" AGNs in the early Universe (Kocevski et al. 2023).

From Figure 9, CEERS_1019 covers a similar location to the CLOUDY star-forming models with log 10 U ∼ -1.9 to -1.5 (for the single-component fit to Hβ; smaller star) and the AGN models with log 10 U ∼ -2.1 (for both the narrow and singlecomponent fits to Hβ; both stars). Similarly, the MAPPINGS models suggest agreement with log 10 U ∼ -2.5 to -1.5; however, this is also dependent upon the pressure assumed. These results add to the expectation that this galaxy has high ionization powering these strong line ratios. Indeed, this agrees well with the strong O32 line ratio measured for this galaxy, which is often indicative of highly ionized, metal-poor gas (e.g., Schaerer et al. 2022;Williams et al. 2023) and a regime that could suggest Lyman continuum leakage (e.g., Izotov et al. 2018). As a comparison to the ionization parameters gleaned from Figure 9, we estimate this value using two relations from the literature. First, using the theoretical relation between O32, 12+log(O/H), and the ionization parameter described in Kobulnicky & Kewley (2004), we derive log 10 ∼ -2.3.

Finally, using the measured relation between O32 and the ionization parameter explored in Papovich et al. (2022), we derivelog 10 ∼ -2.0. These results are all relatively consistent with one another and further highlight the strong ionizing nature of this source.

Discussion

The "chicken or egg" origin of BH seeds in the first galaxies remains unsolved. Theoretical predictions suggest a mix of "light" (∼10 2 M e ) Population III stellar remnants and/or "heavy" (∼10 5 M e ) seeds formed via the direct collapse of primordial gas, as summarized in recent reviews by Smith & Bromm (2019), Inayoshi et al. (2020), andFan et al. (2022), but the relative mix of each seed type remains unconstrained by observations. Regardless of the seeding mechanisms, our Universe is capable of forming extremely massive SMBHs very early, with the highest-redshift massive quasar at the time of this writing being z = 7.642 (Wang et al. 2021), existing 700 Myr after the Big Bang. Observations with JWST, especially MIRI observations of obscured growth (e.g., G. Yang et al. 2023, in preparation) and NIRSpec observations of broad ( 1000 km s -1 ) and/or high-ionization (e.g., N V, He II, [Ne V], and C IV) emission lines can now allow the first real census of AGNs in low-mass (M * < 10 10 M e ) hosts at high redshift, providing strong constraints on the initial BH seed distribution. Direct observations of early SMBHs will not only constrain mechanisms for early BH growth (e.g., Ricarte & Natarajan 2018a, 2018b), but can also provide further insight into the role ionizing photons from AGNs played in the reionization of the IGM (e.g., Finkelstein et al. 2019;Giallongo et al. 2019;Dayal et al. 2020;Grazian et al. 2020;Yung et al. 2021;Grazian et al. 2022). Here we consider plausible formation mechanisms for this source (e.g., small versus large seeds), as well as discuss what its plausible descendants are. While the inferred mass of this BH is not much larger than the proposed range for DCBH seeds, it is unlikely that we are witnessing a DCBH soon after formation. The constraints we place on both the stellar mass and gas-phase metallicity of the host galaxy are significantly higher than expected (DCBH formation requires near-primordial conditions), implying that this object is a "standard" (albeit very distant) AGN accreting at ∼the Eddington rate, with a continuum SED dominated by stellar emission (in line with the stellar emission dominating the continuum SED; Section 5), similar to the scenarios recently explored by Volonteri et al. (2023).

It is therefore interesting to explore how this AGN was seeded: from a low-mass (∼10-100 M e ) stellar seed or a highmass (∼10 4-6 M e ) DCBH seed, summarized in Figure 10. In between these two scenarios could be a seed from a massive starburst cluster or young ultracompact dwarf galaxy (Kroupa et al. 2020) labeled as "Dense Star Cluster Seed" in Figure 10. The purple curves show idealized BH mass tracks assuming a 100 M e stellar seed. Such a seed could plausibly form at z ∼ 30, but would be unable to accrete for ∼100 Myr due to radiative heating of the gas from the stellar progenitor (e.g., Johnson & Bromm 2007;Jeon et al. 2014). Assuming such a 100 Myr delay, the onset of BH growth would begin by z ∼ 18.5. As shown by the dashed purple line, Eddingtonlimited accretion in this scenario would be unable to reach the inferred BH mass by the redshift measured for this object. The enhanced stellar feedback environment the host galaxy must have experienced in the previous ∼100 Myr further complicates this. Building up to the observed Z  0.01Z e requires a large number of supernova explosions in a short period of time. This would lead to violent, turbulent mixing and heating of the gas, making accretion likely very inefficient (significantly sub-Eddington).

However, there may also be periods of super-Eddington ("catch-up") accretion. These would occur in short episodes, with minimal sub-Eddington accretion between them; such short bursts of BH growth could plausibly build up the ∼10 7 M e SMBH we observe in an otherwise stellar-emissiondominated galaxy (e.g., Volonteri & Rees 2005;Madau et al. 2014;Inayoshi et al. 2016). For example, Pezzulli et al. (2016), when examining the evolution of merger tree simulations, find that super-Eddington accretion modes (e.g., Haiman 2004;Silk 2005;Polletta 2008) can be a significant component of SMBH growth in gas-rich environments, where up to 75% of the SMBH growth can be accounted for by periods of super-Eddington accretion, with intermittent phases of disruption resulting from the rapid depletion/replenishment of the bulge gas reservoir out of which the BHs accrete.

In Figure 10, we highlight one potential version of this scenario, with the purple solid line showing an object that starts with a seed mass of 100 M e , which begins growing at z = 18.5 with episodic 10 Myr periods of super-Eddington (10 times) growth, followed by 50 Myr "breathing" periods with sub-Eddington (0.1 times) growth. Periodic super-/sub-Eddington growth allows this hypothetical object to grow to ∼10 7 M e by z ∼ 8.7 from a small seed, matching our observations. While plausible (in particular, the episodic growth with suppressed, sub-Eddington, accretion periods was predicted by the simulations of Jeon et al. 2012 andMassonneau et al. 2023), this scenario is somewhat contrived, making a stellar seed somewhat unlikely.

A DCBH origin may be more plausible. While the inferred metallicity implies the DCBH event occurred significantly prior to the observed epoch, this could provide the needed time for a massive seed ∼10 6 M e to grow by the needed factor of ∼10 times, without also having to invoke super-Eddington conditions. Given the likely strong stellar feedback environment in the assembling galaxy, this DCBH scenario may be favored. We show two plausible DCBH tracks as the red lines in Figure 10. Both options, of a lower-mass 3 × 10 4 M e DCBH forming at z = 15 or a higher-mass 10 6 M e DCBH forming at z = 10.5, should they grow at the Eddington limit, could reproduce the observed mass of this object by z ∼ 8.7.

Last, extrapolating the inferred BH growth tracks to lower redshifts provides insight into the ultimate descendants of early BHs, such as the one we observe here. As shown in Figure 10, this object is both too low-mass and observed too late to be the plausible progenitor of the z > 7 quasar population. However, should near-Eddington accretion be sustained for significant periods of time, objects like this one could plausibly evolve into the massive z ∼ 6 -7 SDSS quasar population. Further constraining both the progenitors and descendants of highredshift SMBHs is possible, as the population of z > 8 SMBHs is coming into view. Given the discovery of this source in a relatively small data set, it is likely that further JWST spectroscopy-in particular, following up bright sources discovered over wide areas with the upcoming Nancy Grace Roman Space Telescope-will yield such a population (e.g., Yung et al. 2023).

Evolution in the BH-Galaxy Mass Relationship

Figure 11 compares CEERS_1019 with other high-redshift (z ∼ 6) AGNs (Izumi et al. 2021;Kocevski et al. 2023) as well as the z = 0 M BH -M * relationship. Because they are easier to detect with pre-JWST instruments, most high-redshift AGNs are high luminosity. However, AGNs have a well-known luminosity-dependent bias (i.e., the Lauer bias; Lauer et al. 2007), such that higher-luminosity AGNs tend to have overmassive BHs.

JWST has improved our ability to select low-luminosity accreting BHs (e.g., Kocevski et al. 2023), which are expected to lie closer to the intrinsic (overall) M BH -M * relationship. Figure 11 shows an estimate from the empirical TRINITY model (Zhang et al. 2023a) for how different bolometric luminosity thresholds bias the median M BH -M * relationship. With an estimated bolometric luminosity of L bol ∼ 10 45 erg/s and a host stellar mass of M * ∼ 10 9.5 M e , TRINITY would suggest that the present source would be expected to have a ΔM BH ∼ 0.7 dex offset from the median relation. The present source is located around the upper envelope of the M BH -M * relation spanned by the AGN sample from Reines & Volonteri (2015) at z ∼ 0. The present source also has a lower M BH compared to the z = 5.55 AGN from Übler et al. (2023), which is ∼10 times brighter. Qualitatively, this is consistent with the Lauer bias, that brighter AGNs tend to be overmassive BHs compared to their host galaxies. If the present source lies on the median M BH -M * for its luminosity, the Lauer bias estimate from TRINITY would suggest that the overall z = 8.7 median M BH -M * would have M BH = 10 6.3 M e at M * = 10 9.5 M e , i.e., there would be no overall evolution from the z = 0 M BH -M * relationship. If other sources are confirmed to have similarly consistent masses compared to the z = 0 M BH -M * relationship, it would cement a tight relationship between galaxies and BHs that extends into the EoR, and would also provide a strong challenge to theoretical models that have predicted substantial evolution in the M BH -M * relation for high-redshift and low-mass galaxies (see Habouzit et al. 2021 for a survey).

The Efficacy of Line Ratio Diagrams in Identifying AGNs at Early Times

Many of the canonical emission-line ratio diagnostics of AGN versus star formation using rest-frame optical emission lines are calibrated for low to moderate redshifts (z  1; e.g., Baldwin et al. 1981;Veilleux & Osterbrock 1987;Juneau et al. 2011Juneau et al. , 2014;;Trump et al. 2015;Backhaus et al. 2022). As such, the divisions between the star-forming and AGN regions defined by these optical line ratio diagnostics are not necessarily valid at higher redshifts.

Several works offer remedies to these low-redshift diagnostics by quantifying a redshift evolution to the division (e.g., Coil et al. 2015;Cleri et al. 2023), but these new divisions have not been extended to the EoR; in fact, these relations for these diagrams are shown to break down at high z in simulations (Hirschmann et al. 2019(Hirschmann et al. , 2022)). The redshift-evolving diagnostics of these works employ a simple shift from the low-redshift divisions; however, there is insufficient knowledge of high-redshift AGNs to validate this practice for galaxies in the EoR.

For a different avenue to remediate these issues, other works have offered completely new emission-line ratio diagnostics with higher-ionization emission-line ratios, e.g., in the optical with He II/Hβ and [Ne V]/ [Ne III] ( Katz et al. 2023;Cleri et al. 2023) and in the UV with, e.g., C III]/He II, O III]/He II, and C IV/He II (e.g., Feltre et al. 2016;Hirschmann et al. 2019Hirschmann et al. , 2022)). Unfortunately, many of these very high ionization lines are often weak and thus may not always have wellconstrained ratios, as is the case for the object studied in this work. Photoionization modeling (e.g., from Cleri et al. 2023), along with comparisons to data across a broad redshift range (0  z  8.5), suggest that sources of very highly ionizing photons (>54.42 eV; Berg et al. 2021) may easily be confounded with AGNs by traditional BPT-style diagnostics. This result has been shown with well-studied z ∼ 0 extreme metal-poor dwarf star-forming galaxies, which have been used as analogs to galaxies in the EoR (e.g., Berg et al. 2019Berg et al. , 2021;;Olivier et al. 2022).

Many recent studies from early JWST data have shown exactly this: high-ionization, low-metallicity star formation may produce line ratios consistent with AGNs in the lowerredshift diagnostics (e.g., Brinchmann 2023; Katz et al. 2023;Trump et al. 2023;Trussler et al. 2022). From the other side, other work has shown known AGNs that produce line ratios consistent with star-forming galaxies (Übler et al. 2023). This growing body of work suggests that the dichotomous classification of a galaxy as either AGN-dominated or starformation-dominated is inadequate to accurately describe a galaxy's ionizing spectrum in the early Universe.

These works suggest a larger "composite" region of BPTstyle diagrams as a function of redshift. This allows for contributions to the ionizing spectrum from multiple sources (e.g., an accreting BH and extreme metal-poor star formation). As more exotic systems in the early Universe are found, evolved versions of these diagnostics and their use with other information about the galaxy (e.g., the analysis of broad Balmer lines and SED fits as shown in this work) will become increasingly necessary to discriminate between sources of ionization.

Conclusions

We present the discovery of an accreting SMBH at z = 8.679, using spectroscopy from NIRSpec and NIRCam/WFSS and imaging from NIRCam and MIRI from the JWST CEERS Survey (Bagley et al. 2023;Finkelstein et al. 2022b;S. L. Finkelstein et al. 2023,in preparation). This source, denoted here as CEERS_1019, was initially identified as a z ∼ 8 photometric Lyα-break dropout candidate by Roberts-Borsani et al. (2016b), with spectroscopic confirmation via Lyα emission using Keck/MOSFIRE by Zitrin et al. (2015). We detect several strong rest-UV and rest-optical emission lines using medium-resolution JWST/NIRSpec spectroscopy (R ∼ 1000) covering 1-5 μm. From this work, we measure a significant broad component of the Hβ emission line with FWHM ∼1200 km s -1 , which we conclude originates from AGN activity.

Our measurements are based on observations from several JWST instruments, but our key results are derived from the 1-5 μm medium-resolution grating spectra from NIRSpec. We use an automated line-finding code (based on Larson et al. 2018) to identify significant emission features in an unbiased and systematic way, running line-injection simulations to estimate accurate line flux uncertainties. Due to the very high S/N of the [O III] 5008 Å emission line, most of our emissionline fits are tied to the redshift and FWHM of this line. We perform more customized fits when required, including for known doublets, as well as observed broad lines. We also note that a few NIRSpec-detected lines are also observable in the NIRCam WFSS spectra, albeit at lower S/N.

Our key result is that the observed Hβ emission line has a clear broad component, comprising ∼half of the emission-line flux, with an FWHM ∼1200 km s -1 . As this broad component is not seen in the stronger [O III] lines (as would be the case in large-scale outflows), we conclude that its origin is from a broadline region around an AGN. This is supported by weak N V, N IV, and Mg II emission, as well as weak broad C III] emission, and the morphology, which shows a compact point source among three Sérsic-like clumps.

We explore the properties of both this accreting SMBH as well as the host galaxy. Constraints from the continuum SED, including photometry from HST, as well as NIRCam and MIRI, show that the continuum emission is dominated by stellar light, particularly in the rest-UV, and that the stellar population is modestly massive (log M/M e ∼ 9.5) and heavily star-forming (log sSFR ∼ -7.9).

From the width and flux of the broad Hβ emission feature, we estimate the mass of the SMBH to be log (M/M e ) = 6.95 ± 0.37 and that it is accreting at 1.2 ( ± 0.5) times the Eddington limit. From the ratios of narrow emission lines, we find that the gas in this galaxy is modestly metal-poor (∼0.1 Z e ) with little dust attenuation, dense, and highly ionized. Similar to other recent JWST results, CEERS_1019 sits elevated over most star-forming models in diagnostic line ratio diagrams. While this could be interpreted as evidence that the AGN significantly contributes to the narrow lines, the presence of other presumably star-forming-dominated galaxies in a similar line ratio regime means that we cannot rule out stellardominated ionization.

We discuss the implications of the presence of this BH early in cosmic history. We find that it is difficult to explain an SMBH of this mass at z ∼ 8.7 with a stellar seed unless periodic episodes of super-Eddington accretion are possible.

Alternatively, Eddington-limited accretion from a massive (∼10 4-6 M e ) DCBH seed could reach the target mass by the observed epoch. Either scenario is somewhat exotic-the uncovering of a larger population of early SMBHs will place further constraints on their seeding and growth mechanisms.

We conclude by noting that while the broad Hβ component is statistically significant (2.5σ) and is significantly required by the fit (Δ BIC ∼ 3), we can only detect this component in one Balmer line. However, the next few months will see a MIRI spectrum obtained by the MIRI Guaranteed Time Observation (GTO) team (PID 1262), covering Hα. If, as we have concluded, a broadline AGN is present in this source, then the upcoming MIRI Hα spectra should show a well-detected broad line.