A gas-enshrouded and gas-reddened black hole at cosmic dawn

Observations and data reduction

MoM-BH*-1 has been observed with JWST by three programmes. It was imaged in cycle 1 by the PRIMER survey¹⁴ (JWST-GO-1837; principal investigator: J. Dunlop) using the MIRI (5 January 2023 and 16 January 2023) and NIRCam (7 August 2023 and 9 August 2023) instruments. In cycle 2 (19 December 2023), the EXCELS survey¹⁵ (JWST-GO-3543) obtained 1.5 h of NIRSpec G395M spectroscopy. In cycle 3 (15 December 2024), we targeted MoM-BH*-1 as part of the ‘Mirage or Miracle’ NIRSpec prism survey (JWST-GO-5224). We included MoM-BH*-1 as a high-priority target second in importance only to luminous z > 10 sources in our UDS masks, because it appeared in several priority target lists—AGN/LRD candidates selected based on compact morphology and template fitting with EAZY⁴⁷, extremely massive galaxy candidates, sources with peculiar red colours and the literature LRD candidates⁴⁸. 

We use the v.7.2 images of the PRIMER field released on the DAWN JWST archive (DJA) reduced using the grizli software⁴⁹. PSF-matched photometric catalogues based on these images were produced in ref. ⁵⁰. We use the public v.3 NIRSpec reductions of the EXCELS grating data from the DJA derived using the msaexp software^37,51,52. The MoM data are reduced with the same pipeline following the same choices.

Although we found no relevant radio or ALMA archival data, MoM-BH*-1 has been observed with Chandra⁵³. Similar to virtually all LRDs^9,54,55, it remains undetected in the X-rays (L_x < 44.5 erg s⁻¹ (1σ) at rest frame 5–90 keV). Key empirical properties of the source are summarized in Extended Data Table 1.

Emission line fitting

We use a custom NIRSpec emission line fitting package⁵⁶ to simultaneously fit emission lines in the grating and prism spectra. The advantage of this approach is that despite the low signal-to-noise ratio (SNR) in either mode, features may be robustly recovered because of their occurrence at the same wavelength across both dispersers.

We first fit the Hβ line and [Oiii] doublet, and then use the redshift as a prior to fit Hγ (Fig. 1). We model Hβ as a single emission line with three absorbers constrained to have negative flux (one at line-centre and two on either side of zero velocity) motivated by the symmetric absorption troughs on either side of the central double-peak (Extended Data Fig. 2). The systemic redshift is tied to [Oiii] and the broad Hβ component, with the absorbers allowed to range freely. The number of absorbers is decided based on the maxima reached in the reduced χ², which is similar for three and four absorbers, but we opt for parsimony. This large number of absorbers may be merited to account for secondary peaks at ±2,400 km s⁻¹ (Extended Data Fig. 2). Including narrow Hβ at the systemic redshift leads to completely unconstrained flux degenerate with absorption and no improvement, so we neglect this component. Resulting fits, in which we sample the posterior with the NUTS sampler implemented in numpyro⁵⁷, are shown in Extended Data Fig. 1 and reported in Extended Data Table 2.

Although the formal errors on the derived fluxes and line widths of many of the components is significant, the key features relevant to this analysis are robustly recovered—extremely broad Hβ emission, and deep, broad absorption wiping out about 25% of the emission flux, including approximately 100% of the flux at line-centre. Another notable aspect of these fits is the location of the two absorbers—they are recovered at very similar velocities, but on either side of line-centre (approximately ±1,500 km s⁻¹), albeit with significant uncertainties ((-1,53{2}_{-113}^{+345},mathrm{km},{{rm{s}}}^{-1}) and (+1,55{6}_{-1,378}^{+232},mathrm{km},{{rm{s}}}^{-1})). We explore this symmetry, which extends not only to the location of the absorbers but also to the detailed structure of the entire emission line profile over a few 1,000 km s⁻¹ in Extended Data Fig. 2.

Morphology

We use pysersic⁵⁸ to fit a Sersic profile to the imaging. We focus on the F356W and F444W imaging, in which the source is well-detected. We follow the procedure described in ref. ⁵⁹—we build an empirical PSF using stars in the field, and then use this PSF with pysersic to sample the posterior Sersic parameters with numpyro. The source is unresolved, and we are able to place a 99% upper limit on the effective radius of <117 pc consistent with the BH-dominated interpretation.

Cloudy modelling

That extremely dense gas might blanket the LRDs was already inferred in the work that defined this class of sources⁹—about 10% of these LRDs showed clear signs of Balmer absorption, with this being a lower limit due to resolution and SNR¹¹ as opposed to the ≪0.1% in pre-JWST AGN⁶⁰, with only a handful examples⁹. Local AGN⁶¹ do not display Balmer breaks anywhere as strong as MoM-BH*-1 (about 8), with the strongest breaks (≲2.5) occurring in quiescent galaxy AGN that reach the stellar population maximum shown in Fig. 2b. Inspired by Balmer absorption in LRDs, and using Cloudy models for continuum absorption⁵, demonstrated that the very gas producing strong Balmer absorption likely also produces Balmer breaks. However, note the breaks in these models are weaker than what we require to explain MoM-BH*-1 and are abrupt instead of the smooth rollover seen in this source. Here we build on this work.

The grid of parameters we explore is summarized in Extended Data Table 3. This grid spans a wide range to capture the extreme spectral shape at hand. We start with an intrinsic AGN continuum SED that is parametrized using a series of power laws and with a ‘big bump’ temperature³². This AGN continuum is passed through gas surrounding the central source that is defined in terms of its gas density, column density, metallicity and turbulent velocity. The irradiation of the cloud is modulated by an ionization parameter (log(U)). A large turbulent velocity of the absorbing gas 500 km s⁻¹ is motivated by the width of the absorption lines we observe (see Fig. 1). This implies a high Mach number, which is consistent with recent models of AGN disks⁶²—whether such high turbulence can be sustained across a large envelope remains to be seen. The turbulence is important in producing a smooth rollover instead of a sharp break⁶. Finally, we apply a uniform dust screen⁶³ with A_V = 0−3 in the post-processing—this is the only post-processing step we apply.

To select models that are consistent with the data, we require: (1) Hβ in net emission with 30 < EW [Å] < 45; (2) Hγ in absorption with −5 < EW [Å] < 0 Å; (3) a strong Balmer break and high optical-to-UV ratio such that ({f}_{{rm{4.5mu m}}}^{lambda }/{f}_{{rm{1.8mu m}}}^{lambda } > 6) and ({f}_{{rm{4.5mu m}}}^{lambda }/{f}_{{rm{2.8mu m}}}^{lambda } > 3); and (4) fluxes in MIRI bands within 2σ of the observations. The few thousand models that satisfy these constraints are re-simulated at higher resolution, retaining only hydrogen that is relevant to the key features for simplicity and speed. Of these, we select the one that closely follows the detailed shape of the continuum. Furthermore, we experiment with the covering factor and distance between the gas and central source, and find these quantities to be degenerate with the AGN SED and ionization parameter. We find the ‘net transmitted’ flux (that is, the sum of the attenuated incident continuum and diffuse continua or lines) produces a good match to the data, not the ‘total’ flux (which also includes reflected continua or lines).

Key features of fiducial Cloudy model: blanketed by gas, not by dust

First, we emphasize that this model (Fig. 3 and Extended Data Fig. 3) is a model that matches various features of interest based on our limited grid. The relevant parameter space that we have done our best to sample is high-dimensional, degenerate and much remains unknown (for example, the intrinsic SEDs of early AGN with effects such as photon trapping). Therefore, we caution against detailed inferences beyond the feasibility of the broad physical picture that we present here (an accretion disk embedded in dense gas).

Notably, the selected model features an extreme column density (about 10^25.8 cm⁻²) comparable to the most enshrouded systems observed⁶⁴ and a high gas density (10¹¹ cm⁻³) conducive to Balmer absorption. The turbulent velocity of 500 km s⁻¹ is commensurate with the width of the central Hβ absorber. Metal-poor gas expected of a dwarf galaxy at z ≈ 8 is preferred. The AGN slope parameters are within the range of literature SEDs^6,65.

An important feature of the fiducial model is that it is virtually dust-free (A_V = 0.15 mag). The UV-weakness arises entirely because of the extreme hydrogen opacity. This is a crucial constraint on LRD models that typically invoke significant amounts of dust to suppress strong UV emission from classical AGN^21,24,66 or to explain the weakness of Hβ relative to Hα¹². Although significant A_V ≈2–3 helps dense gas AGN models produce a smooth Balmer break in the rest-optical⁶, this is ruled out by longer wavelength constraints in our source (Extended Data Fig. 3) and more generally by stringent infrared constraints on the LRDs^34,35. We note that given the very low A_V ≈ 0, the exact details of the dust geometry (for example, screen compared with clumps) do not have a bearing on our results.

Resonant Balmer scattering and insights from Lyα-like shell models

The Hβ profile of MoM-BH*-1 bears a remarkable resemblance to double-peaked Lyα^67,68,69. This motivates us to explore whether Hβ is behaving like a resonant line under the high densities in MoM-BH*-1. It may be the case that this is a common phenomena in the LRDs, but that a symmetric double-peaked line with large separation is easier to observe in MoM-BH*-1 because of negligible Hβ emission from the sub-dominant host galaxy.

Although Hβ normally has an easy cascade escape route, saturating the 2p state can effectively trap Hβ photons, which may be plausible through Lyα pumping in these optically thick environments⁷⁰. We test this idea with the COLT⁷¹ radiative transfer code with which we implement a simple shell model (directly analogous to Lyα shell models) for Hβ. The dominant parameters in these calculations are the thermal velocities of the inner and outer shells, the relative velocities between the shells and the optical depth encountered by Hβ. Strong turbulence is not accounted for in these toy models.

In Extended Data Fig. 4, we show how a relatively narrow Gaussian Hβ line manifests as a double-peaked profile. We confirm the basic trends seen with the Lyα line hold up for Balmer resonant scattering and can serve as a powerful guide to interpret MoM-BH*-1. With this simple toy model, we are able to generate a match for the velocity separation and intensity of the strongest peaks seen in the MoM-BH*-1 Hβ profile.

If scattering is underway, this has some important implications for the physics of the situation. The observed width of Balmer lines (especially when fitted as an absorption system to a broad Gaussian line; that is, the most common baseline model^9,21,72) may not trace the kinematics of the broad-line region but are a consequence of resonant scattering. Therefore, SMBH masses in these systems that are based on line widths of the Balmer lines may be severely overestimated by up to 2 dex (Extended Data Table 4). More sophisticated treatments of resonant scattering effects in LRDs⁷³ differ in assumptions and details, but make the same basic point.

SMBH properties and ramifications for LRD parameters

Given the unique conditions (for example, unusual gas density) in this source, we must exercise care in applying the local scaling relations typically used to derive SMBH properties. Here we explore implications for a variety of approaches to back out the SMBH properties that are summarized in Extended Data Table 4. Given the significant systematic uncertainties, these calculations must be seen as order-of-magnitude estimates to bracket the range of possibilities.

First, we simply apply Hβ-based local scaling relations³⁹ to the absorption-corrected and dust-corrected Hβ line. The absorption correction is typically done by emission line fitting assuming an underlying Gaussian or Lorentzian^9,21 (as in Extended Data Fig. 1). For the dust correction, note that the Hα and Hβ line ratio¹³ is challenging to use given that their observed fluxes reflect radiative transfer in the dense gas. The standard approaches used in the literature at the moment infer a significant A_V of a few magnitudes either from continuum slope fitting or through SED modeling^6,10,12,13. Going by the optical continuum slope¹⁰ would imply an A_V ≈ 2 for MoM-BH*-1 implying a BH mass of about 10^8.3M_⊙. This is comparable to the stellar mass of the host galaxy (<10^8.5M_⊙ at 95% confidence), that is, this is an apparently ‘overmassive’ BH relative to local scaling relations between host galaxy mass and BH mass^13,41.

However, a key insight from the SED of MoM-BH*-1 is that the A_V required to explain the LRDs is likely negligible as the SED is intrinsically UV-weak below the Balmer break. That is, the LRDs are red not because of obscuration from dust, but because of opacity from gas. Accounting for this, with A_V = 0, the SMBH mass derived now is about 5 × 10⁷M_⊙ with a luminosity of approximately 15% the Eddington limit.

Next, we note that the observed line width (after correcting for absorbers) may not be faithfully tracing the broad-line region, violating the basic ansatz for scaling relations (Extended Data Figs. 4 and 5). The underlying BH emission line profile in LRDs may be as complex as the one shown in Extended Data Fig. 2, but is perhaps difficult to disentangle from strong emission from the host. The complex structure may signify radiative processes that are modifying the emission arising from the broad-line region. As an example of one such process, if Hβ is undergoing resonant scattering through the atmosphere around the BH, then the intrinsic BLR width before scattering could be as low as around 600 km s⁻¹ (Extended Data Fig. 5). This would yield an SMBH mass of approximately 10⁶M_⊙.

Finally, we are able to constrain the BH properties using our Cloudy model. The bolometric luminosity directly follows from integrating the SED shown in Extended Data Fig. 3 (about 10^44.5 erg s⁻¹). We note that this is an order of magnitude lower than the value expected from local calibrations, further underscoring their inapplicability. We can convert this luminosity into a mass by observing that theoretical models of similar scenarios^74,75,76 predict close to Eddington or super-Eddington accretion to sustain these convective envelopes. Furthermore, the ubiquity of OI emission thought to arise from Lyβ fluorescence in LRDs^21,66,77 is thought to be a hallmark of super-Eddington accretion⁷⁸. Inspired by these lines of evidence, requiring L/L_bol ≈ 1, therefore, allows us to estimate the BH mass as approximately 10^6.3M_⊙. Correcting the potentially overestimated literature BH masses by about 10−100× would bring the typical M_BH/M_⋆ of JWST AGN (about 1–10%; refs. ^79,80,81) closer to the ratio in the local Universe (0.01%; ref. ⁸²).

Host galaxy properties

To constrain the host mass, we can leverage the insight from clustering analyses that the rest-UV light in the typical LRD (FWHM = 1,000−2,000 km s⁻¹), on average, originates almost entirely from the host⁴¹. Of course, this is not true for all LRDs, particularly the most luminous sources with higher FWHM broad lines that display AGN signatures even in the rest-UV²¹, but so far they seem the exception. To construct an empirical M_⋆−M_UV relation, we use the compilation of low-luminosity galaxies at z = 3−7 from the All the Little Things (ALT) survey⁴⁰ in the Abell-2744 field. This is the largest spectroscopic sample of M_UV < −15 galaxies at these redshifts. The ALT stellar masses are derived with the Prospector SED fitting code applied to 27 bands of NIRCam + HST photometry, including all JWST medium and broad bands^83,84. We estimate that a source with M_UV between −17.9 and −18.3 has a median stellar mass of (log ({M}_{star }/{M}_{odot })=7.{4}_{-0.3}^{+0.6}) and a 95% upper limit of (log ({M}_{star }/{M}_{odot }) < 8.5).

Hints of variability

It is of particular interest to test for signs of variability in MoM-BH*-1 as an independent constraint on the physics of the source. In LRDs, the SED is a summation of the host and BH at all wavelengths to varying degrees, but in this case, the stark outshining of the host galaxy means an ‘undiluted’ variability signal may be stronger and more easily detected. One of the handful reports of variability in an LRD to date has been in a source with a large Balmer break implying a high BH fraction³⁰.

Three sets of observations covering 3−5 μm exist for MoM-BH*-1 separated by ≈ 60 days in the rest-frame (Extended Data Fig. 6). It is certainly not ideal to test for variability across three different observing modes with distinct systematics. Nevertheless, it is notable that the source appears to have brightened by 30 ± 7% between the first epoch (NIRCam) and third epoch (NIRSpec prism, no post-processing renormalization to photometry). The prism spectrum suffers from slit losses, and yet this source appears brighter than expected from the NIRCam photometry. It is the only source across the 136 sources with high-SNR photometry (SN > 10) and spectra (50^th percentile of SN > 5) observed in two masks as part of the MoM program that shows this degree of brightening. We also note that typically, the NIRSpec/G395M flux is systematically  ≈ 10–20% lower than the prism flux due to calibration uncertainties⁸⁵ – accounting for this, both the prism and NIRSpec fluxes for MoM-BH*-1 appear to be in agreement and  ≈ 30% brighter than the NIRCam flux. The shape of the SED (and the depth of the Balmer break) is consistent across the NIRCam and prism data (for example, ({f}_{{rm{F410M}}}^{nu }/{f}_{{rm{F356W}}}^{nu }=2.0pm 0.2) in the prism vs. 2.2 ± 0.2 in NIRCam). The difference is in the absolute brightness.

Taken at face value, these measurements imply a ≈4σ detection of variability over a mere two months, marking MoM-BH*-1 as an excellent target for future monitoring campaigns. This may be independent evidence not only for the AGN nature of the source, but also for the volatile environment that exists around it.