© The Authors 2024

1. Introduction

Evidence for supermassive black holes (BHs), with masses ranging from a few million to several billion solar masses, has been found in the nuclei of most galaxies in the local Universe. The tight relation with many of the host galaxy properties, and in particular with the central velocity dispersion, has been regarded as indication of co-evolution between BHs and their host galaxies (Kormendy & Ho 2013; Greene et al. 2020). Models and cosmological simulations envisage different possible co-evolutionary (as well as non-co-evolutionary) scenarios, possibly involving galaxy and BH mergers, as well as mutual self-regulation via feedback processes (e.g. Sijacki et al. 2009; Volonteri 2010; Valiante et al. 2016; Inayoshi et al. 2020; Greene et al. 2020; Trinca et al. 2022; Fan et al. 2023; Volonteri et al. 2023; Bennett et al. 2024; Sassano et al. 2023; Koudmani et al. 2022). Most of these scenarios are degenerate in explaining the scaling relations observed locally. Yet, at high redshift different models and simulations expect different properties for the population of BHs and their relations with their host galaxies (e.g. Visbal & Haiman 2018; Valiante et al. 2018; Schneider et al. 2023; Volonteri et al. 2021; Habouzit et al. 2022; Trinca et al. 2022). Therefore, in order to validate, test, and discriminate between different scenarios, it is crucial to explore the population of (accreting) BHs at high redshift, along with their host galaxies.

The search and characterisation of accreting BHs and their host galaxies at high redshift has made tremendous progress during the past 20 years (e.g. Merloni et al. 2010; Bongiorno et al. 2014; Trakhtenbrot et al. 2017; Mezcua et al. 2018; Lyu et al. 2022). However, at z > 4, until recently, observational constraints limited the identification and characterisation of sources primarily to the very luminous quasar regime (see Fan et al. 2023, for a review). Within this context, the discovery of BHs with masses in excess of several billion solar masses already in place at z > 6–7 has been unexpected (e.g. Bañados et al. 2018; Wang et al. 2020), since models and cosmological simulations found it challenging to reproduce the growth of such massive BHs within the limited amount of time since the Big Bang. Different scenarios have been invoked, such as direct collapse black holes (DCBHs), the merging of stars and BHs in nuclear clusters, and super-Eddington accretion from stellar mass BH seeds, possibly originating from ‘Population III’ remnants (Inayoshi et al. 2020; Greene et al. 2020; Ferrara et al. 2014; Trinca et al. 2022; Volonteri et al. 2023; McKee & Tan 2008; Banik et al. 2019; Sassano et al. 2021; Singh et al. 2023; Haidar et al. 2022; Ni et al. 2022; Weller et al. 2023; Beckmann et al. 2023; DeGraf & Sijacki 2020). Testing and differentiating between these different scenarios requires the detection and characterisation of smaller BHs at high redshift.

Major progress has been made thanks to JWST (Gardner et al. 2023; Rigby et al. 2023). Indeed, while several AGN candidates have been identified through JWST imaging and broad-band, photometric spectral energy distribution (Furtak et al. 2023; Onoue et al. 2023; Barro et al. 2024; Yang et al. 2023a; Bogdán et al. 2024; Juodžbalis et al. 2023), JWST spectroscopy has revealed broad-line AGNs at high redshift with moderate to low luminosities. Specifically, both NIRSpec Multi-Object Spectroscopy (MOS) and Integral Field Spectroscopy (IFS) observations, as well as slitless NIRCam grism spectroscopy, have revealed a population of AGNs at z > 4 and out to z = 10.6 with luminosities (10⁴⁴ − 10⁴⁵ erg s⁻¹; Kocevski et al. 2023; Übler et al. 2023, 2024; Harikane et al. 2023; Matthee et al. 2024; Maiolino et al. 2024a; Juodžbalis et al. 2024a; Scholtz et al. 2023) lower than those of classical quasars (≳10⁴⁶ erg/s). The estimated BH masses are between 10⁶ and 10⁸ M_⊙, significantly lower than those typically inferred for quasars at similar redshifts. Interestingly, based on their narrow line ratios, these systems would not be classified as AGNs in classical diagnostic diagrams, such as the BPT diagram (Baldwin et al. 1981), since they are primarily located in the region populated by star-forming galaxies in the local Universe. This offset compared to their lower redshift counterpart, in the above diagrams, is primarily interpreted as being due to high-redshift AGNs being hosted in a low-metallicity environment (Kocevski et al. 2023; Übler et al. 2023).

Although less luminous, but much more common than quasars, these early AGNs are likely playing an important role in the evolution of their host galaxies by exerting feedback processes (Koudmani et al. 2022). An intriguing example is the detection of a prominent AGN-driven outflow in the most distant AGNs, GN-z11 (Maiolino et al. 2024a), which is observed ejecting gas and metals in its circumgalactic medium (Maiolino et al. 2024b), while also heating and ionising it (Scholtz et al. 2024). These phenomena may result in rapid suppression of star formation and lead to the early emergence of quiescent galaxies, or contribute to short-term quenching and to the burstiness of star formation (Carnall et al. 2023a,b; Looser et al. 2024, 2023; Dome et al. 2024; Strait et al. 2023).

There have been differing claims about the number of AGNs in early galaxies, with fractions ranging from 1% to 5% (Harikane et al. 2023; Matthee et al. 2024). Even more unclear is the potential contribution of AGNs to the re-ionisation of the Universe, with some estimates claiming that they could contribute up to 50% (Giallongo et al. 2019; Harikane et al. 2023) and others indicating that they are unlikely to contribute significantly (Matthee et al. 2024).

In this paper, we present the discovery of a sample of 12 new broad-line AGNs at z > 4 in the first DEEP tier and two of the MEDIUM tiers of the JADES survey (Eisenstein et al. 2023), by using the NIRSpec MOS spectroscopic observations. These observations are deeper than previous observations and were performed with multiple dispersers providing different spectral resolutions. Therefore, these datasets enable us to unveil AGNs with a diversity of broad line widths and often in a lower luminosity regime than previous surveys, either associated with lower-mass BHs, lower accretion rates, or more obscured AGNs. We show that this data uncovers the properties of the early phases of BH formation and their connection with their host galaxies, some of which nicely confirm expectations from models and simulations, and others of which are unexpected and prompt further theoretical modelling.

Throughout this work, we use the AB magnitude system and assume a flat ΛCDM cosmology with Ω_m = 0.315 and H₀ = 67.4 km/s/Mpc (Planck Collaboration VI 2020). With this cosmology, 1″ corresponds to a transverse distance of 5.84 proper kpc at z = 6.

2. Sample, observations, and data processing

2.1. Observing strategy and target selection

The data used in this paper have been obtained as part of the JADES survey (Eisenstein et al. 2023). This survey combines nearly 800 hours of NIRCam, NIRSpec-MOS, and MIRI observations in parallel mode, in the GOODS-S and GOODS-N fields. This is a multi-tiered survey reaching different depths (down to AB ∼ 30.5 in imaging and AB ∼ 29 in spectroscopy) in multiple bands, multiple dispersers, and over different areas (for a total of 175 arcmin²). An extensive description of the survey is given in Eisenstein et al. (2023). Here, we only discuss briefly the three specific spectroscopic tiers that are used in this paper: Deep/HST in GOODS-S, Medium/JWST in GOODS-N, and Medium/HST in GOODS-N.

A detailed description of the target selection and of the spectroscopic observations is given in Eisenstein et al. (2023) and Bunker et al. (2024). Here, we only summarise that in the Deep/HST tier spectroscopic targets were selected giving higher priority to the highest redshift candidates, according to their photometric redshifts, and primarily relying on the Lyα dropout signature, and then gradually lower priorities to galaxies at lower redshifts. The targets selected in Deep/HST were primarily obtained from previous HST imaging. However, a number of high-priority targets were also selected from NIRCam imaging (Rieke et al. 2023) obtained shortly before the NIRSpec observations. We note that the parent HST targets list was leveraging previously published catalogues, which may have discarded high-z AGNs based on their ‘stellarity’; in other words, their point-like appearance. Moreover, the initial selection from NIRCam images may have discarded some high-z AGNs because of their colours being similar to brown dwarfs. Therefore, the Deep/HST tier is likely biased against (type 1) AGNs.

The observations in Deep/HST were obtained with three dithered configurations of the Micro Shutter Array (MSA) (Jakobsen et al. 2022; Ferruit et al. 2022; Böker et al. 2023), each with 3-shutters nodding. The low-resolution prism, the three medium-resolution gratings (G140M/F070LP, G235M/F170LP, G395/F290PL), and the high-resolution G395H/F290LP grating were used, for a total observing time of up to ∼27 hours with the prism and up to ∼7 hours with each of the gratings (the specific exposure for each target depends on whether they could be accommodated in the three dithered MSA configurations or not). Overall, spectra for about 250 galaxies were obtained in Deep/HST (Bunker et al. 2024).

In the case of the Medium/JWST GOODS-N tier, the targets were primarily selected from NIRCam images, with a similar set of selection criteria and priorities as for Deep/HST, although targeting brighter magnitudes on average. In contrast to previous observations in JADES, this was the first tier that was not biased against the selection of AGNs (this is the reason why we chose it for this paper). Actually, three targets were specifically selected because of their imaging (very compact) and colour properties, similar to those of other AGNs identified in other JWST observations (Übler et al. 2023; Harikane et al. 2023; Matthee et al. 2024). Specifically, to identify these specific candidate AGNs, we initially considered well-detected sources (F444W > 250 nJy) with F090W < F200W < 3×F444W. From this initial selection, we inspected the SEDs (to reject unreliable measurements and obvious brown dwarfs). We further inspected the cut-out images to reject contaminated sources and sources with proper motions. A considerable fraction of these candidates (> 30%) has one or more blue sources or features within 0.1–0.2 arcsec, similar, for example, to GS-3073 (Übler et al. 2023). The final sample prioritises objects that have nearby, blue companions first, because these objects were considered to have the lowest risk of being brown-dwarf contaminants.

The observing strategy in the Medium/JWST tier was very similar to that used for Deep/HST, but with shorter overall exposures, resulting in 2.6 hours in each of the five dispersers (prism, the three medium-resolution gratings, and G395H). Four different pointings were planned in GOODS-N for this tier, but one of them failed twice because of a telescope guide problem and MSA shorts; hence, the fourth pointing is planned for a repeat observation in 2024. Overall, in the three successful pointings, 712 targets were observed.

Finally, we also explored spectra from the Medium/HST GOODS-N tier. This consisted of four pointings (two dithers) with targets selected from HST imaging, with selection criteria similar to Deep/HST, but brighter magnitudes. In this tier, we used the prism (1.7 hours on source) and the three medium-resolution gratings (0.8 hours on source for each of them). About 660 sources were observed in this tier.

GN-z11 is part of the JADES sample and was observed both in the Medium/HST and Medium/JWST tiers in GOODS-N. It was specifically targeted because previously identified as a galaxy at z ∼ 11 based on HST and ground-based observations (Oesch et al. 2016). The first spectrum was obtained by Bunker et al. (2023) (within Medium/HST), and additional MSA (Medium/JWST) and IFS observations were obtained by Maiolino et al. (2024a,b) and Scholtz et al. (2024), which identified it as a type 1 AGN (specifically an analogue of NLSy1). Although specifically observed because its previously known properties, it is part of the JADES sample and it is included in the analysis of this paper.

2.2. Data processing

The processing of the MSA data is also described in Bunker et al. (2024) and other MSA-JADES papers (e.g. Carniani et al. 2024). A detailed description of the data processing will be presented in Carniani et al. (in prep.). Here, we only recall that we have used the pipeline developed by the ESA NIRSpec Science Operations Team and the NIRSpec Guaranteed Time Observations (GTO) Team. As we are primarily interested in the detection of the broad component of the Balmer lines emitted from the unresolved broad-line region of AGNs, we used the spectra extracted from the central 3 pixels of each 2D spectrum (i.e. the central 0.3″), which maximises the signal-to-noise (S/N) for compact sources. We also mention that the pipeline automatically corrects for the wavelength dependent slit losses (by also taking into account the source position in the slit), by assuming a point-like source. The latter assumption is certainly appropriate for the radiation coming from the AGN, but more broadly also for their host galaxies, as these systems are extremely compact, as will be discussed in the next sections.

2.3. Spectral resolution

We finally mention that the nominal resolution of the NIRSpec dispersers (R ∼ 30−300 for the prism, R ∼ 700−1300 for the medium-resolution gratings, and R ∼ 1900−3500 for the high-resolution gratings) applies only in the case of uniformly illuminated shutters. This is rarely the case for galaxies at z > 4, which are generally more compact than the shutters’ width (0.2″), even when convolved with the telescope’s point spread function (PSF) for most wavelengths, and it is certainly true for the broad-line region, which is totally unresolved. As a consequence, the effective resolution of MSA observations is generally significantly higher, and primarily driven by the telescope’s PSF rather than the shutters’ width.

Although the resolution slightly depends on the position on each quadrant (due to the slightly varying PSF) and on the position of the target in the shutter, for a point source the effective resolution spans the following ranges (de Graaff et al. 2024): R ∼ 60−500 for the prism, R ∼ 1200−2500 for G140M, R ∼ 1100−2300 for G235M, R ∼ 1300−2000 for G395M, and R ∼ 3600−5400 for G395H.

3. Identification of broad-line active galactic nuclei

In this section we describe the identification of broad-line AGNs and their basic properties. As was mentioned, we focus on galaxies at z > 4, as this is the new regime probed specifically by JWST.

3.1. Criteria for the detection of a broad-line region

The presence of a broad-line region (BLR) is assessed via the detection of a broad component of either Hα or (at z > 7) Hβ line emission, without a broad counterpart in the forbidden transitions (in particular the bright [OIII]5007).

We fit the Balmer lines with two Gaussian component, narrow and broad. The [NII] doublet around Hα is forced to have the same width and velocity as the narrow component of Hα and, similarly, the [OIII] doublet near Hβ is forced to have the same width and velocity as of the narrow Hβ. For each of the [NII] and [OIII] doublets the intensity ratios are fixed by to the corresponding Einstein coefficients.

While for each line we simultaneously fit the spectra in the available dispersers, we do not simultaneously fit the Hα+[NII] and Hβ+[OIII] groups because of two main reasons: small wavelength calibration uncertainties (primarily associated with uncertainty in the location of the target in the slit, both because of residual astrometric uncertainties and because of MSA target acquisition uncertainties of ∼0.05″, Jakobsen et al., in prep.; Carniani et al., in prep.), may result in slight artificial wavelength shifts and dispersion solutions (and also slightly different resolutions) for the two groups; secondly, due to the different PSF at the two wavelengths, the spectrum may sample slightly different regions of the host galaxy. While we do not fit the Hα+[NII] and Hβ+[OIII] groups simultaneously, we checked that the two separate fits are fully consistent within uncertainties.

In order to claim the detection of a BLR, we require the second, broader component of the Balmer lines to be at least a factor of two broader than the narrow component and have a significance of at least 5σ. Furthermore, we require that the Bayesian information criterion (BIC) parameter (Liddle 2007), defined (in the case of Gaussian noise) as:

$$\begin{array}{c}\hfill \mathit{BIC}={\chi }^2+klnn\end{array}$$

(where k is the number of free parameters and n is the number of data points), for the fit with the broad component is at least a factor of six smaller than the value for fit with only the narrow component; that is,

$$\begin{array}{c}\hfill BI{C}_{only-Narrow}-BI{C}_{Broad+Narrow}=\mathrm{\Delta }BI{C}_{\mathit{NB}}>6.\end{array}$$

We conservatively mark two cases with 6 < ΔBIC_NB < 10 as ‘tentative’ (IDs 3608 and 62309) although their broad components are detected at > 5σ; their removal from our analysis would not change the conclusions.

Finally, we note that while the broad lines of low- and intermediate-luminosity type 1 AGNs are often well fitted with a Gaussian profile (e.g. Marziani et al. 2019), for very luminous quasars it has been suggested that a double power law profile may reproduce better the profile of the BLR permitted lines (Nagao et al. 2006). Our targets are certainly not in the category of luminous quasars and indeed the broad Hα is not fit better with a double power law profile. This aspect is not of particular interest in the context of the detection of the broad components, but it is relevant for the interpretation of the complex Hα profiles, which will be discussed in the next section.

3.2. Ruling out outflows

For many objects, the case for a BLR is pretty much obvious, with a prominent and nearly symmetric broad component of Hα as in classical type 1 AGNs. However, in cases in which the broad wings of the Balmer lines are fainter or asymmetric, these could in principle be associated also with galactic outflows, and indeed the JADES survey has revealed a number of galactic outflows by inspecting the emission line profiles (Carniani et al. 2024). However, high-velocity ionised gas in the host galaxy, especially if associated with outflows, should be seen even more prominently in the profile of metal lines, especially the strong [OIII]5007 transition. Indeed, the higher excitation of the gas in outflows, along with the fact that galactic outflows are naturally more metal enriched than the host galaxy, typically make the [OIII] line stronger than both Hα and Hβ, even in outflows at high redshift (Holden & Tadhunter 2023; Carniani et al. 2015; Marshall et al. 2023). Therefore, a requirement for the identification of BLR is that the broad component should not be detected in the [OIII] line.

We finally note that it is very unlikely that, despite being brighter, an [OIII] outflow component is not detected, while seen in Hα, because of dust extinction. Indeed, being out of the galactic plane, and also quite extended, is generally less obscured by dust in the galactic disc. Dust extinction is even less prominent at z > 4, where most galaxies show little or no dust reddening (Fiore et al. 2023; Sandles et al. 2024). Specifically, the galaxies with identified broad component of Hα in our sample have Hα/Hβ Balmer decrements generally consistent with the case B value of 2.8, or only slightly higher by a factor of less than 1.5; this indicates, that even in the extreme case of a putative outflow obscured at the same level as the galaxy ISM, the associated broad wings of [OIII] (typically a factor of at least a few stronger than Hα) would still be detectable.

3.3. Broad-line active galactic nuclei identified in JADES

In the three JADES sub-tiers analysed by us we identify 12 new broad-line AGNs at z > 4 (in addition to GN-z11 reported by Maiolino et al. 2024a, which required a dedicated analysis). These are all identified via the detection of a broad component of Hα. Although we identified some possible cases of Hβ with broad component at z > 7, these did not pass our criterion for the identification of a BLR.

Figure 1 shows the spectral region around Hα for the 12 newly identified AGNs. For sake of simplicity in most cases we show only the G395M spectra, which are those showing the broad Hα visually seen more clearly, while generally the higher-resolution grating spectra (if available for Hα) have a consistent profile, although noisier. The only exception is ID073488 for which we directly show the G395H spectrum as it has very high S/N (but we shall also show the corresponding G395M spectrum in the next section). In Fig. 1 all spectra are fitted with a narrow and a single broad Hα component, as detailed below; however, for three objects we shall discuss a more complex Hα profile fitting in the next section.

Fig. 1. Medium-resolution spectra of the 12 new JADES galaxies with evidence for a broad component of Hα ascribed to the BLR of an AGN. We show a zoom around Hα. The solid blue line shows the spectrum along with the errors (light blue shaded area). The dashed red line shows the total multi-component fit; the dashed blue and purple lines show the narrow and broad components of Hα, respectively, while the dashed green lines show the components fitting the [NII] doublet (often undetected).

In Fig. 1 the observed spectra are shown with a light blue line (and errors with light blue shading), while the dashed lines show the various components. Specifically, the narrow component of Hα is shown with a dashed dark blue line, the [NII] doublets (which we recall are forced to have the same width as the narrow component of Hα and to have doublet ratios fixed by the Einstein coefficients) are shown with dashed green lines (often undetected), and the broad component of the Hα lines is shown with a dashed purple line. The total modelled profile, including all the above components and a power-law continuum, is shown with a dashed red line.

The flux, widths and shifts of the main lines of interest resulting from the fit are listed in Table 1, and in Table 2 for the three cases with more complex Hα profiles. The ΔBIC_NB values, as was described above, are also reported in Table 1 (in Table 2 the meaning of ΔBIC is different, as is discussed in the next section).

For the mentioned two tentative detections, IDs 3608 and 62309, the significance of the broad component of Hα is 8.1σ and 5.4σ, respectively, and the ΔBIC_NB are 6.2 and 8.0, respectively. However, for all other AGNs in the sample the detection of the broad component of Hα is between 11σ and 67σ, while the ΔBIC_NB are between 50 and 330.

In Appendix A we also show the spectra around the Hβ+[OIII] region, primarily to illustrate the absence of a broad counterpart of [OIII], hence excluding an outflow origin. A broad component of Hβ is only detected in the case of ID 954 at z = 6.759, which is the most luminous AGN in our sample.

Grating spectra of [OIII] are not available for IDs 8083 and 53757. However, in these two cases the identification of a BLR is unambiguous given that the Hα broad component profile is nearly symmetric (in contrast to outflow profiles, which are generally blueshifted) and the flux of the broad component is similar or higher than the flux of the narrow component, which would imply that the ionised gas in the outflow is more massive than the whole ionised ISM in the host galaxy.

3.4. NIRCam images of broad-line active galactic nuclei

Fig. 2 shows the NIRCam images (Rieke et al. 2023) of the selected broad-line AGNs (each thumbnail is 3″ × 3″ in size). In most cases these targets are very compact, often dominated by a central point source, as is expected for type 1 AGNs. Some of them have red colours and a point-like appearance similar to those identified by Matthee et al. (2024) and found to be hosting reddened AGNs. We shall see in Sect. 7 that indeed the AGNs in these targets tend to be reddened, although not by a large amount. Along with the central point source tracing the AGN, more extended structures are often also seen.

Fig. 2. NIRCam images of the selected broad-line AGNs. The following false-colour coding was adopted: blue – F115W; green – F200W; and red – F444W. Each thumbnail is 3″ × 3″ in size.

We have quantified the presence of an underlying host galaxy by using the ForcePho software to perform such a point source (i.e. the central AGN) and host-galaxy decomposition in our sample. ForcePho (Johnson et al., in prep.) fits multiple PSF-convolved profiles simultaneously to all spectral bands by sampling the joint posterior distribution via Markov Chain Monte Carlo (MCMC). The fitting of the components was done in the individual NIRCam tiles and observations, even before their combination and mosaicing, which allows a much more accurate control of the PSF and avoids issues with correlated noise in the mosaics. This software and methodology has already been successfully employed in similar cases using the same set of NIRCam images (e.g. Tacchella et al. 2023; Baker et al. 2023; Robertson et al. 2023).

In those cases for which a host galaxy could be detected by ForcePho, the resulting radii and Sérsic indices are reported in Table 3. Clearly in most cases the host galaxies of these AGNs are extremely compact, typically with effective radii of a few 100 pc. The Sérsic indices are typically disc-like (n ∼ 1), with the exception of ID 8083 which has an extremely high Sérsic index, and which may indicate the presence of an early, compact spheroid.

These parameters will be useful both for assessing the dynamical mass and for interpreting the velocity dispersion in Sect. 8.

4. Candidate merging black holes

4.1. Identification of candidate dual active galactic nuclei

Among the galaxies for which we have identified a broad component of the Balmer lines, three cases stand out for the peculiar profile of Hα, which is not fitted properly with a single broad component. The first (and most distant) one is ID 10013704 at z = 5.919, in Deep/HST GOODS-S, whose medium- and high-resolution spectra are shown in Fig. 3. In addition to the narrow component, having the same width as [OIII] (140 km/s), the medium-resolution grating shows a clear, broad component of Hα with a width of 2400 km/s (Fig. 3, top left panel). Not surprisingly, such a broad component is nearly lost in the noise in the high-resolution spectrum (Fig. 3, bottom left panel). However, the medium-resolution spectrum also shows a prominent redshifted wing of the narrow component. Such a wing is seen also in the high-resolution spectrum, where, rather than a wing, it is more clearly resolved as a bump, slightly redshifted by ∼250 km/s and a width of 415 km/s. The central and right panels of the same figure show the medium- and high-resolution spectra after subtracting the narrow and one of the two broad components, to better illustrate the significance of the remaining broad component.

Fig. 3. Spectra around Hα of ID 10013704, the highest-redshift AGN showing an indication for a dual BLR. The top and bottom panels show the medium- and high-resolution spectra, respectively. The line coding is the same as in Fig. 1, but in this case the dashed violet line shows the intermediate broad component that is needed to properly reproduce the observed profile. The central and rightmost panels show the spectrum from which the narrow components, as well as one of the two broad components have been removed, to better highlight the significance of the other broad component.

Table 2 gives the best fit parameters for the two components of Hα that are broader than the narrow component, showing that the very broad component is detected at > 10σ and the second broad component is detected at > 6σ. In this case the ΔBIC given in the table is not the difference of the BIC with only a narrow component and after adding a broad component (which is very large), it is instead the difference of the BIC with only one broad component and with two broad components:

$$\begin{array}{c}\hfill \mathrm{\Delta }BI{C}_{1B2B}=BI{C}_{1BLR}-BI{C}_{2BLRs}\end{array}$$

In the case of ID 10013704, ΔBIC_1B2B = 77, clearly indicating the strong need of a second broad component.

Neither the very broad component nor the second broad and slightly redshifted bump are seen in the bright [OIII]5007 profile (a detailed comparison of the [OIII] and Hα profiles is discussed in Appendix B and shown in Fig. B.1), implying that it is very unlikely that either of them is associated with an outflow, as galactic outflows, as is discussed in Sect. 3.2, are generally more metal enriched than the host galaxy and, especially if AGN-driven, more prominent in [OIII]. Moreover, asymmetries in the wings associated with outflows are typically blueshifted, as a consequence of dust extinction in the host galaxy disc, which absorbs preferentially the receding (redshifted) component of the outflow. It is also very unlikely that the bump is tracing star formation in a merging companion. The lack of broad [OIII] emission would imply extremely low metallicity (< 0.01 Z_⊙, Curti et al. 2024; Laseter et al. 2024; Vanzella et al. 2023), but the width of the hump (415 km/s) would indicate that the putative merging galaxy, despite being extremely metal-poor, is much more massive than the primary galaxy, which would completely contrast with any formulation of any mass-metallicity relation (Maiolino & Mannucci 2019; Curti et al. 2020).

We suggest that the second broad, redshifted bump in ID 10013704 is tracing the BLR of a second, fainter AGN, probably associated with a secondary BH in the process of merging with the BH in the primary galaxy.

We note that morphologically ID 10013704 shows the presence of a possible weak tidal tail (Fig. 2), indicating that it very likely went through a recent merger, and that the BH of the secondary galaxy is now approaching the nucleus of the primary galaxy, while accreting gas and being detectable as a secondary AGN¹.

A similar case is seen in the spectrum of ID 73488 at z = 4.133, in the Medium/JWST tier in GOODS-N, as illustrated in Fig. 4. Also in this case there are a clear narrow (68 km/s) and a very broad (∼2400 km/s) component of Hα. However, the line profiles, especially in the high-resolution spectrum, cannot be explained with only two components and reveal the presence of an intermediate component with a width 460 km/s (the results of the double component fitting are shown in Fig. 4 and reported in Table 2). In this case, the velocity shift of the intermediate component is much smaller, but the S/N is much higher than in the case of ID 10013704. Unfortunately, in this case, we do not have the high-resolution grating covering [OIII], yet this transition is observed in a spectral region where the G235M grating has a fairly high resolution (resolution FWHM ∼ 150 km/s, or σ ∼ 63, for compact sources), and the presence of the intermediate component in [OIII] should be clearly seen if present (at S/N > 20 even in the unlikely case of F(Hα) = F([OIII])), but is undetected. Therefore, following the same arguments as for ID 10013704, we suggest that this is an additional case of two BLRs associated with a dual AGN.

Fig. 4. As Fig. 3 but for the additional candidate dual BLR AGN ID 73488.

Also in the case of ID 73488 both broad components are detected at very high significance (> 35σ), as reported in Table 2. Introducing the second broad component of Hα the BIC improves by more than 500 (ΔBIC_1B2B = 547; Table 2).

We note that, being the additional broad component less shifted in velocity relative to ID 10013704, this could be a case in which the profile of Hα is better reproduced with a double power law. However, in Appendix C we illustrate that this is not the case. While in that appendix we rule out the double power law scenario more quantitatively, here we simply notice that the need for two BLR separate components is visually clear from the Hα profile in the high-resolution spectrum (Fig. 4 bottom left), which shows a clear inflection of the profile at about ±400 km/s from the line centre.

Finally, in Fig. 5 we show an additional case of candidate dual BLR in the galaxy ID 53757. Unfortunately, in this case we do not have any grating spectrum of [OIII]. However, the intensity of the broad component of Hα (more than two times the flux of the narrow component) makes it very unlikely that it is associated with outflow, as it would imply more gas mass in the outflow than in the host galaxy. Moreover, the nearly symmetric profile of Hα (if anything slightly redshifted, as is discussed below), argues against the outflow interpretation (which generally requires a blueshift of the high velocity gas). Also in this case, as for ID 73488, the Hα broad profile, especially the inflection seen in the high-resolution spectrum strongly suggest the presence of a second BLR. Additionally, the second broad Hα component is slightly redshifted (by ∼100 km/s) relative to the narrow component, as in the case of ID 10013704. The result of the fit with two BLRs is shown in Fig. 5 and the resulting parameters reported in Table 2: the second broad component, with width ∼700 km/s, is detected at > 5σ, while the broader component, with FWHM ∼ 2800 km/s, is detected at 10σ. Introducing the second broad component results into a ΔBIC_1B2B = 39.

Fig. 5. As Fig. 3 but for the additional candidate dual BLR AGN ID 53757.

In terms of morphologies, similarly to ID 10013704, both ID 73488 and ID 53757 show the presence of weak asymmetric features, which are likely remnants of recent mergers, hence fitting in the scenario in which these systems may be hosting dual BHs, in the process of merging.

We note that these three cases are not ‘double-peaked’ broad lines that have sometimes been invoked in the past as possible signatures of dual AGNs, and whose interpretation is ambiguous since BLRs in disc-like configurations can potentially display a double-peaked profile (Eracleous et al. 1997; Eracleous & Halpern 2003; Krolik et al. 2019). Instead, what we detect in these three objects are broad components of the Hα line with very different profiles, on top of each other.

Other recent works on the profile of the Hβ of quasars have been fitted with multiple components (Yang et al. 2023b), however these are systems that (beside being in a totally different luminosity regime than our targets) have a Hβ profile that is heavily blended with very broadened [OIII] lines and with strong FeII multiplet emission. Therefore the profile of the Hβ from the BLR cannot be really disentangled properly from the other emission features, including Hβ outflowing components. Bosman et al. (2024) fit the Hα of a z ∼ 7 quasar with a double Gaussian, however: (1) it may well be that this is a double BLR too (JWST is indeed revealing that high-z quasars are nearly ubiquitously in merging systems, so the possibility of dual BH is very high); (2) they have not attempted a double power law profile fit, which is generally more appropriate for quasars; (3) they do not have access to [OIII], hence they cannot assess whether one of the two broad components is due to outflowing gas. In our case, the absence of any broad component in the [OIII] profile rules out that any of the two BLR components is actually associated with an outflow.

In summary, while it might be possible to envisage peculiar geometries of a single BLR which can mimic these profiles, we are not aware of models predicting a narrower profile on top of a broader profile around a single BH. Even if possible from future models, we think that the interpretation of a dual AGN is a plausible one.

It is interesting that in galaxies ID 1093 (z = 5.59) and ID 3680 (z = 5.27) the broad component is significantly blueshifted. These could potentially be additional cases in which an accreting BH is merging, while the BH of the primary galaxy is not actively accreting. Alternatively, in these cases the two putative BHs already merged and the resulting BH received a recoil velocity kick, which is expected to happen after BH coalescence (Blecha et al. 2011, 2016; Civano et al. 2010; Chiaberge et al. 2018; Morishita et al. 2022).

In all three cases presented in this paper the putative dual AGN must be located within the MSA shutter (0.2″) hence they must be separated by less than ∼1 kpc. There is no clear evidence for a double nucleus in the NIRCam images of these galaxies, indicating that the separation is likely less than about 500 pc (which is the NIRCam projected resolution at 2 μm at z = 5). However, the nucleus with the smaller BH has lower luminosity (as is discussed in the next section) and it is possible that it is outshone by the larger BH in the NIRCam images, even if at separation larger than 500 pc. However, we are clearly probing a regime of dual AGNs much closer and much less luminous than dual quasars found at lower redshifts (z ∼ 0.5 − 3) via imaging and spatially resolved spectroscopy (e.g. Mannucci et al. 2022, 2023; Ciurlo et al. 2023; Scialpi et al. 2024). Yet, NIRSpec-IFU spectroscopy is starting to reveal dual AGNs spatially resolved on scales of a few 100 pc at higher redshift (Übler et al. 2024), providing further support that dual accreting BHs on relatively small scales may not be uncommon.

4.2. Fraction of merging black holes

We have found three candidate merging BH out of the 11 AGNs in the redshift range 4 < z < 6, which is where we have the best statistics. This is excluding the two AGNs with shifted BLR, which may also be merging BHs (in which the primary BH is inactive), but which may also be recoiled BHs after merging.

As is discussed in the next section, the smaller BHs have masses of about 10⁶ M_⊙, while the companion larger BHs have masses of a few to several times 10⁷ M_⊙.

Multiple simulations and semi-analytical models have predicted merging BHs and dual AGNs at various cosmic epochs (Di Matteo et al. 2023a,b; Chen et al. 2023; Volonteri et al. 2020, 2022; Barausse et al. 2020; Barai et al. 2018; Valentini et al. 2021; Vito et al. 2022; Di Mascia et al. 2021; Mannerkoski et al. 2022). It is however difficult to compare our findings of dual, close-pair AGNs with predictions from simulations, as most of them provide predictions of dual AGNs and/or merging BHs at lower redshifts (z < 4) and/or at larger separations (> 1 kpc) and/or more massive BHs (M_BH > 10⁷ M_⊙) (Chen et al. 2023; Volonteri et al. 2020, 2022; Barausse et al. 2020; Di Matteo et al. 2023a). Alternatively, other works explore a wider range of BH masses and redshifts, but only provide the expected merger rates (Sesana et al. 2004, 2007), for the goal of establishing the detectability of their gravitational wave signals, which are difficult to compare with our results.

However, it is useful to compare our findings with simulations in other regimes published so far. At redshifts z ∼ 2–3, the simulations mentioned above predict a fraction of dual AGNs of only about 1–2%, while this fraction increases to 10–20% if considering BH pairs in which only one of the two BH is accreting (hence identified as AGN). This may appear in tension with our fraction of ∼25% dual AGNs (i.e. in which both BH are accreting). However, assessing in simulations whether a BH is accreting or not, at the time of observation, is regulated by assumption at the sub-grid physics level, and therefore there is probably scope to increase the fraction of expected active dual BH in simulations to levels similar to those inferred by us.

Moreover, while the smallest BHs in our putative merging systems have masses of only 10⁶ M_⊙ or a few times 10⁵ M_⊙, the simulations do not provide predictions on the merging fraction for BHs smaller than 10⁷ M_⊙; so the fraction of dual AGNs predicted in simulations would probably increase significantly if including lower-mass BHs.

Finally, the fraction of dual AGNs and BH mergers is seen to increase at high redshift in some of these simulations.

It is important that future simulations provide predictions specifically in the luminosity, M_BH, and redshift ranges explored by us, in order to compare exactly the same regimes.

It is also important that future spectroscopic surveys, possibly including IFS observations, assess the fraction of merging BHs, with higher statistics and expanding the luminosity and BH mass ranges, also to guide the comparison with ongoing and future gravitational wave experiments (e.g. Amaro-Seoane et al. 2012, 2023; Agazie et al. 2023)

4.3. Comparison with previous surveys

In their sample of ten broad-line AGNs, from NIRSpec medium-resolution spectroscopy, Harikane et al. (2023) do not find evidence for peculiar Hα profiles that could be ascribed to dual BLRs. However, they lack information from high-resolution spectroscopy, which we have seen to be crucial to identify and explore the presence of dual BLRs. Specifically, in all cases (possibly with the exception of ID 10013704) we would have not identified the presence of a peculiar BLR profile, if it was not for the profile seen in the high-resolution grating.

Matthee et al. (2024) used NIRCam slitless spectroscopy, at resolution R ∼ 1600, to identify 20 broad-line AGNs. Their Hα broad lines do show some peculiar profiles that could be in principle be ascribed to dual BLRs. We note that the spectral profile measured through slitless spectroscopy is convolved with the spatial distribution of the line emission. If the spatial profile of a line emitting source is not point-like, then it can be difficult to discriminate an intrinsically complex spectroscopic profile from a profile induced by spatially resolved distribution of the nebular emitting gas. However, there are two cases in their sample with double peaked broad Hα profiles, which are intriguing. Matthee et al. (2024) interpret these as due to Hα absorption. However, Hα is not a resonant line and the n = 2 is not a metastable level. As a consequence, it is difficult to see Hα in absorption, as it requires extremely high densities, temperatures of ∼10⁴ K and high column densities. Indeed, although a dozen cases of Hα absorption have been seen in some type 1 AGNs and quasars (e.g. Williams et al. 2017; Shi et al. 2016; Zhang et al. 2015; Schulze et al. 2018), these are extremely rare at low/intermediate redshifts, less than about 0.1% of the AGN population. It might be unlikely that 10% of the type 1 AGNs found by Matthee et al. (2024) are characterised by such a rare phenomenon. It is possible that the peculiar broad Hα profiles might be instead associated with dual AGNs. Yet, before making any claim on the Hα profile of these targets, slit or IFS spectroscopy would be needed, so to avoid any potential profile artifacts resulting from the slitless spectroscopy, as was discussed above. Within this context we note that the detailed spectral analysis of some AGNs discovered by JWST is actually revealing Balmer absorption lines that are clearly not associated with a double BLR (e.g. Juodžbalis et al. 2024b), but their occurrence has yet to be assessed. This clearly prompts for more extensive and more detailed spectroscopic observations (especially at higher resolution) of the type 1 AGNs that are being discovered by JWST.

5. Stellar masses and dust extinction

We have used the low-resolution prism spectra from NIRSpec, along with the multi-band photometry from NIRCam, to model the ultraviolet (UV) to optical rest-frame spectral energy distribution with BEAGLE (Chevallard & Charlot 2016). We have added a power-law continuum component to BEAGLE to account for the contribution of an AGN. We fit the spectra with a delayed-exponential star formation history, with a burst of constant star formation lasting 10Myrs prior to observation. We model the dust attenuation of the star-forming component with the two-component dust law of Charlot & Fall (2000). The simple model of the emission from the accretion disc consists of a single power-law component and is parameterised by the slope and the fractional contribution to the luminosity at a rest-frame wavelength at 1500 Å. Since BEAGLE does not include broad-line region models, we mask all emission lines with significant emission in the prism spectra, limiting the fits to the shape of the continuum. Since the accretion disc emission is likely to be reddened, we included extinction of the power-law component with an SMC dust-law Pei (1992), which is often found to be appropriate for high-z AGNs (Richards et al. 2003; Reichard et al. 2003). To avoid full degeneracies between reddenning and the power-law slope, we tried fitting with two fixed power-law slopes, allowing the attenuation to vary freely. Specifically, we use both the power-law of a standard Shakura-Sunyaev accretion disc – β = −7/3 ≈ −2.33 (with Fλ ∝ λ^β) (Shakura & Sunyaev 1973) – as well as the power-law inferred from the Sloan composite quasars template obtained by Vanden Berk et al. (2001); that is, β = −1.56. The latter is redder than the Shakura-Sunyaev slope as probably it already incorporates some dust reddening (in addition to the convolution with the accretion disc turnover at high energies in the individual quasar spectra). Of course, quasars and AGNs with redder slopes are observed, but generally ascribed to various degrees of reddening (Richards et al. 2003). Therefore our choice of two power laws and dust reddening, should cover most of the observed cases, while miniminzing degeneracies. Both power laws generally give acceptable fits. In the following we use the average stellar masses and attenuation inferred from the two cases, and the errorbars will reflect the difference between the two cases (in addition to the errors associated with the individual fitting).

Note that in the case of candidate merging BHs the AGN component was modelled with a single (reddened) power law, as attempting to include two separate power laws would result in strong degeneracies.

Table 3 lists the stellar masses and dust attenuation inferred from BEAGLE, as was discussed above. In Appendix D we show the BEAGLE fits to the prism spectra. The inferred stellar masses span from ∼10⁸ M_⊙ up to a few times 10¹⁰ M_⊙. We shall discuss how, when compared with the dynamical masses, in the case of ID 954 the stellar mass is probably overestimated (likely because of this is the most luminous AGN in our sample and therefore it is more difficult to disentangle host galaxy and AGN components), while for ID 11836 the stellar mass might be underestimated (but see discussion in Sect. 8.2).

The inferred dust extinctions are similar to those inferred by Harikane et al. (2023) for the AGNs in CEERS and GLASS, but significantly lower than those inferred by Matthee et al. (2024). The latter result may indicate that, as a consequence of the JADES selection function, we may be missing type 1 AGNs that are very reddened. This has implications for the AGNs census, as will be discussed later on.

In Table 3 we also include GN-z11, whose properties were inferred in Maiolino et al. (2024a) and Tacchella et al. (2023). As was mentioned, GN-z11 will be included in the rest of the analysis, as part of the JADES sample in GOODS-N (Medium-HST and Medium-JWST tiers).

6. Black hole masses and accretion rates

We can estimate the BH masses in these systems by using the local virial relations, which link BH masses with the width of the broad lines and the continuum or line luminosity. We specifically use the relation provided by Reines et al. (2013) and Reines & Volonteri (2015), which provides the BH mass in terms of width and luminosity of the broad component of Hα:

$$\begin{array}{cc}\hfill log\left(\frac{M_{\mathrm{BH}}}{M_{\odot }}\right)& =6.60\hfill \\ \hfill & +0.47log\left(\frac{L_{\mathrm{H}\alpha }}{{10}^{42}\mathrm{erg}/\mathrm{s}}\right)+2.06log\left(\frac{{\mathrm{FWHM}}_{\mathrm{H}\alpha }}{{10}^3\mathrm{km}/\mathrm{s}}\right)\hfill \end{array}$$(1)

The advantage of using this relation is that we can consistently compare with the local scaling relations provided by Reines & Volonteri (2015). When comparing with the results from other surveys, for consistency we shall re-calculate the BH masses by using this relation. The Hα luminosity used in Eq. (1) was corrected for dust extinction, as inferred in the previous section.

Obviously, the local virial relations, such as Eq. (1), are derived locally and there is no guarantee that they apply also to high redshift AGNs and quasars. These relations are associated with the small scale (< pc) physics and dynamics around the BH, and there is no reason to think that this would change with redshift. The only potential concern is the drop in metallicity and the associated reduction of the dust content. Indeed, it is thought that the radius of the BLR is primarily set by the dust sublimation radius, which gives the squared root dependence on the luminosity in the virial relations. However, in the dusty ‘torus,’ the medium remains optically thick to the UV radiation even if the dust-to-gas ratio is very low, because of the extremely large column densities (typically > 10²³ cm⁻², Risaliti et al. 1999), while the sublimation radius remains set by the dust micro-properties, which do not seem to change drastically in these low-mass systems at z ∼ 4–7 (Witstok et al. 2023a,b). Even if there was some change in the dust properties (which may be possible in much more massive and more luminous systems Maiolino et al. 2004; Gallerani et al. 2010), the weak (squared root) dependence of the BH mass on the AGN luminosity in the virial relations would not affect strongly the estimates of the BH masses at high-z.

The inferred BH masses are reported in Table 3 and span from 8 × 10⁷ M_⊙ down to 4 × 10⁵ M_⊙. The uncertainties on the BH masses include the propagation of the errors of the quantities involved in calculating them; however, we also take into account the scatter of the virial relations, which contribute to about 0.3 dex (in quadrature) to the BH mass uncertainty.

The BH associated with BLR1 of ID 10013704, at z = 5.9, has the lowest BH mass, 4 × 10⁵ M_⊙. Interestingly, this is in the range expected for DCBH (10⁴ − 10⁶ M_⊙). Obviously, this does not mean that such a small BH is a DCBH, but it is nonetheless encouraging that with JWST we are starting to probe this regime. More statistics in this range from future data, and in comparison with expectations from models, may provide constraints on the BH seeding scenarios.

We can potentially also infer the AGN bolometric luminosity from the monochromatic luminosity of the AGN continuum at a given UV or optical rest frame wavelength and then using bolometric corrections (e.g. Netzer 2019; Duras et al. 2020; Saccheo et al. 2023). However, this requires a proper deblending of the (reddened) AGN component from the light of the host galaxy, which subject to significant uncertainty. Moreover, in the case of dual AGNs, it is not possible to disentangle the contribution to the continuum from the two AGNs. An alternative is to use the scaling relations between the luminosity of the broad component of Hα and the AGN bolometric luminosity, which also allow the disentangling of the two companions in the case of dual AGNs. Specifically, we have used the scaling relation between broad Hα (extinction corrected) and AGN bolometric luminosity provided by Stern & Laor (2012), as given by their Eq. (6). The resulting AGN luminosities are reported in Table 3.

We can in principle also compare the inferred bolometric luminosity with the Eddington luminosity of the BHs, although one should be aware that, since the Hα luminosity has been used to calculate both quantities, there are unavoidably spurious correlations, although the BH mass (used to infer the Eddington luminosity) primarily depends (quadratically) on the width of the line and its dependence on the Hα luminosity is only with the power of 0.5. Aware of these caveats, Table 3 reports the inferred L/L_Edd ratios, and these are also visually reported in Fig. 6, together with the associated BH masses. We warn that the apparent anti-correlation between L/L_Edd ratio and BH mass in Fig. 6 is primarily resulting from the fact that the BH mass is included at the denominator of the y-axis quantity. Hence, Fig. 6 should only be used for a quick visualisation of the distribution of Eddington ratios and BH masses. Results from previous JWST spectroscopic surveys (Übler et al. 2023; Harikane et al. 2023; Matthee et al. 2024) are also shown with blue symbols, as is indicated in the legend (we only show results with a significance larger than 3σ and for which the relevant information is available).

Fig. 6. Distribution of BH masses and Eddington ratios (L/LEdd) for the broad-line AGNs in JADES (large golden circles). The dashed orange lines connect candidate dual AGNs. We also show the results from other JWST surveys with blue symbols (see legend, only detections at > 3σ are shown). Note that the apparent anticorrelation is probably spurious, as the BH mass is at the denominator of the Y-axis quantity. The plot has simply the purpose of visually illustrating the distribution of the two quantities.

Despite the uncertainties, it is clear that most BHs in our sample accrete at sub-Eddington rates, mostly L/L_Edd < 0.5. However, the smallest BHs, with masses around or below 10⁶ M_⊙, tend to accrete at Eddington or super-Eddington rates. One of these is GN-z11, already discussed in Maiolino et al. (2024a), and the other are two of the candidate small merging BHs.

As has already been discussed in Maiolino et al. (2024a), these small BH accreting at high-z are very similar to local/low-z Narrow Line Seyfert 1’s, which are indeed found to have low-mass BH, hosted in low-mass galaxies, accreting at super-Eddington rates (Greene & Ho 2004; Mathur et al. 2012)

Our finding of such small and vigorous BHs at high redshift, may suggest that the early phases of BH accretion happen through Eddington or super-Eddington phases, as suggested by various models (Trinca et al. 2022; Bennett et al. 2024). It is also interesting to note that simulations of dual/merging BHs predict that the smaller BH might often accrete at a higher rate than the more massive one (Chen et al. 2023).

However, it is important also to be aware of possible observational biases. For instance, the detection of small BHs with accretion rates significantly below Eddington would likely make them undetectable in our spectroscopic observations. Therefore, although the finding of low-mass BH at high-z with Eddington/super-Eddington accretion rates is exciting, the lack of their sub-Eddington counterparts might be an observational effect.

Finally, we note that the lack of BHs with masses below a few times 10⁵ M_⊙ is also an observational effect, as such small BHs would have a broad component of Hα (only a few 100 km/s) which would be difficult to disentangle from the narrow component of the host galaxy and from even low-velocity outflows.

7. Overmassive black holes relative to the host galaxy stellar mass

The top left panel of Fig. 7 shows the BH mass as a function of stellar mass of the host galaxy. The small red circles show the distribution of local active galaxies obtained (consistently with the same calibration) by Reines & Volonteri (2015); their best-fitting relation is shown by the solid black line, while the shaded grey region shows the dispersion and uncertainty of the fit. This is probably the best local relation for comparing the M_BH − M_star with our results, as Reines & Volonteri (2015) use the same method as ours for estimating the BH masses in AGNs, and also because their host galaxies are mostly star-forming, late-type systems as in our sample (given that the bulk of the host galaxies of the AGNs in our sample have Sérsic index of ∼1, as is shown in Table 3). The red diamonds show additional measurements provided by Greene et al. (2020); the best-fit relation provided by them is shown with the blue line and blue shaded region (giving the dispersion and uncertainty of the relation); we caution that Greene et al. (2020) uses different BH estimations relative to our prescription, so the comparison with our results has this additional uncertainty factor. For completeness, we also show the local relation by Kormendy & Ho (2013) with dash-dotted line, which is however for a sample made primarily of early type galaxies, hence probably not directly comparable with our sample in terms of stellar mass (Sturm & Reines 2024) given that our AGNs are primarily hosted in late type galaxies, but may be more adequate for comparing with dynamical masses, as is discussed in the next section. The blue squares indicate high redshift broad-line AGNs recently discovered with JWST by various AGN surveys, where we report only those detections that are at least more significant than 3σ and which have a stellar mass reported (Kocevski et al. 2023; Harikane et al. 2023; Übler et al. 2023; Ding et al. 2023; Bogdán et al. 2024; Goulding et al. 2023). The new JADES broad-line AGNs, including GN-z11, are shown with large golden circles. The vertical, dashed orange lines connect the candidate merging BHs.

Fig. 7. BH mass as a function of host galaxy properties, specifically: Stellar mass (top left), velocity dispersion (bottom left), and dynamical mass (bottom right). The JADES results presented in this work are shown with large golden circles. Dashed orange vertical lines connect candidate merging BHs. Blue symbols indicate measurements from other JWST surveys at high-z, as is indicated in the legend (only detections at > 3σ are shown). Grey triangles are measurements of high-z QSOs using ALMA data. Small red symbols show the distribution of local galaxies, as is indicated in the legend; straight lines show the local relation fits (with shaded regions providing the scatter and slope uncertainty), using the samples that are best matching our high-z systems, in each panel (see text for details). In the top left panel, small orange circles show the maximum stellar masses estimated for a few JADES AGNs (see text) and the dotted violet lines show constant values of the MBH/Mstar ratio. In the bottom left panel, small squares indicate the effect of correcting for rotation velocity broadening within the slit.

Most BHs at these early epochs are significantly over-massive relative to the stellar mass in their host galaxies, when compared with the local scaling relation. This was already found by previous surveys (blue squares), but for the newly identified AGNs this phenomenon becomes even more prominent. In some extreme cases we even find BH masses approaching the stellar masses of their host galaxies (see also Bogdán et al. 2024; Goulding et al. 2023).

The strongest deviations occur at low stellar masses (M_star < 10⁹ M_⊙), although we caution that at such low masses the local scaling relation is actually poorly explored and we mostly rely on the extrapolation from higher masses.

Such a strong offset may be partly due to the stellar masses being significantly underestimated. However, even if we use masses obtained without accounting for the AGN in the continuum modelling (i.e. assuming that all continuum emission is due to stars), the JADES AGNs are still located well above the local relation; this is illustrated by the orange circles in the top left of Fig. 7, which show the maximum stellar mass inferred in this way (for sake of clarity we only show those few cases for which the maximum stellar mass exceeds the mass inferred with the standard stellar+power law method by more than 0.1 dex).

BHs being overmassive relative to the local M_BH − M_star relation may be partly due to selection effects. Specifically, given a scatter around the M_BH − M_galaxy relation, the sensitivity limit of quasar/AGN surveys may favour the detection of more massive BHs (which, for a given average L/L_Edd, are more luminous) (Willott et al. 2005; Lauer et al. 2007). Although this effect was thought to be less important at low AGN luminosities (Izumi et al. 2019), especially below the quasar regime, it may still play a role, as has also recently been highlighted by Volonteri et al. (2023). The BH selection bias on this relation was explored more recently and more in detail by Li et al. (2024) and Juodžbalis et al. (2024a) specifically for the JWST-selected broad-line AGNs. Li et al. (2024) claim that the offset can entirely explained with a M_BH − M_star relation consistent with the local one, but a much larger dispersion (about 1 dex) and selection effects on the BH mass. However, Juodžbalis et al. (2024a) illustrate that the finding of overmassive BHs that are dormant (i.e. with a very low accretion rate) does not fit in the scenario outlined by Li et al. (2024). Additionally, in the next section we shall show that the JWST-discovered broad-line AGNs are fairly consistent with the local M_BH − σ relation (where σ is the velocity dispersion) – if the offset on the M_BH − M_star were dominated by selection effects then the same offset should be seen on the M_BH − σ relation. There are also other issues related to the space density of AGN that makes unlikely the selection bias being the dominant effect in the M_BH − M_star offset; however, this is beyond the scope of this paper and will be discussed in a separate, dedicated work (Maiolino et al., in prep.).

Finally, it is interesting to note that BHs above the local scaling relation at early epochs, and even approaching M_BH ∼ M_star, are expected by various theoretical models, especially scenarios envisaging BHs accreting at super-Eddington rates or with very high efficiency, and are also expected to be associated with heavy seeds (i.e. DCBHs) (Trinca et al. 2022; Schneider et al. 2023; Koudmani et al. 2022).

8. Whether the M_BH−σ and M_BH − M_dyn relations are fundamental and universal

8.1. The black hole–velocity dispersion relation

In the local Universe, studies of the BH scaling relations with the properties of the host galaxies, have unambiguously found that the tightest relation is with the central (stellar) velocity dispersion (Kormendy & Ho 2013; Terrazas et al. 2017; Piotrowska et al. 2022), suggesting that this is the most fundamental (causal) relation, while other relations may simply be an indirect by-product.

As has been discussed, most of the AGNs in our sample also have high-resolution spectra covering both [OIII]5007 and Hα, or at least one of these two bright emission lines. For compact sources the spectral resolution in terms of velocity dispersion is about 25 km/s (FWHM ∼ 60 km/s) in the spectral range of interest (de Graaff et al., in prep.). We can therefore accurately measure the velocity dispersion of the gas in the host galaxies of these AGNs (taking into account the instrumental resolution, although this is generally a small correction). This is not necessarily the same as the stellar velocity dispersion used in the local scaling relations; however, various studies have shown that the ionised gas velocity dispersion can be used as a good proxy, although it may require some small correction (Bezanson et al. 2018). Following a similar approach as in Übler et al. (2023) for another JWST identified AGN at z = 5.55, we implement an upward correction of 0.12–18 dex (Bezanson et al. 2018) to the gas velocity dispersion in order to get a close estimate of the stellar velocity dispersion.

The bottom left panel of Figure 7 shows the resulting relation between BH mass and velocity dispersion. As in the top left panel, the red small circles indicate local galaxies from Kormendy & Ho (2013), Bennert et al. (2021) and Woo et al. (2015), while the solid black line shows the best fit relation from Bennert et al. (2021), with the shaded region showing the uncertainty and dispersion. Red diamonds show the new local measurements reported by Greene et al. (2020), while the blue line and shaded region show the local relation and dispersion obtained by them. The previous JWST measurement on a single AGN from Übler et al. (2023) is shown with a blue diamond, while our new JADES AGNs are shown with large golden circles. Also in this case, the vertical, dashed orange lines show the candidate merging BHs.

The AGNs at z ∼ 4–6 are mostly located on the same relation as local galaxies, most of them fully consistent with the local scatter, at least when compared with the local relation obtained by Bennert et al. (2021). There is a small offset relative to the relation obtained by Greene et al. (2020), although still within the scatter, and anyhow the offset is far smaller than for the M_BH − M_star relation. The only exception is the more massive component of the putative merging system ID 73488 (the leftmost JADES point in the figure), which is significantly above the relation, but whose smaller companion is, interestingly, on the relation. However, even Bennert et al. (2021) report some local cases with a similar large deviation (see their lowest velocity dispersion point in the figure).

The fact that the M_BH − σ relation at z > 4 is essentially consistent with the local relation, confirms that this relation is not only fundamental, but also universal, as it holds across most cosmic epochs, at least out to z ∼ 6. The implications of this finding for the co-evolution of galaxies and BHs are profound. It is beyond the scope of this paper to explore the physics behind this relation and why it is so stable across the Universe. Here, we only comment that, combined with the previous finding of the dissolving relation between M_BH and M_star at high-z, these results indicate that the BHs at the centre of galaxies have essentially no knowledge of, nor connection with, the star formation history of their host galaxies, while they are connected with the mass assembly history.

The tight relation between BH mass and velocity dispersion is often interpreted as a causal link between the formation of BHs and the merging history of the central part of the galaxy, responsible for the formation of the spheroidal component (Kormendy & Ho 2013). The finding that such a relation was already in place around the epoch of re-ionisation, suggests that early merging and the early spheroidal formation (Baker et al. 2023) was linked to the early evolution of BHs. Our finding of candidate merging BHs, presented in Sect. 4 consistently supports this scenario, although quantifying the scenario requires a close comparison with simulations.

We conclude by noting that our measurements of the line-of-sight gas velocity dispersion include contributions from both the random motions from the gas as well as, potentially, projected rotational velocities of the gas within the slit. This velocity measurement is therefore a good approximation of the second moment of the velocity distribution of the system, which enters the virial theorem (e.g. Cappellari et al. 2006). As a consequence, the measured velocity dispersion is a good proxy of the gravitational potential and dynamical mass (as we shall discuss in the next section). However, measurements of the velocity dispersions in local galaxies (and used in the local M_BH − σ relation) span a wide range of projected apertures, from 100 pc scale to kpc scale (hence overlapping with the scales probed by us), and may or may not include the contribution from unresolved rotation. Therefore, it is not clear whether our measurements should be corrected for the effect of unresolved rotation in order to be compared with the local relation. However, we also explore the extreme case in which we should only consider the intrinsic velocity dispersion. Carniani et al. (2024) estimated the effect of line broadening due to potentially unresolved rotation curves, when information on the mass of the galaxy, radius and inclination is available. Based on the parameters of the host galaxy that we could derive from imaging and from prism spectroscopy (Sects. 3.4 and 5), we estimate the line broadening from putative rotation and subtract it in quadrature to the observed gas velocity dispersion. The resulting values are shown with large offset orange squares in the bottom left panel of Fig. 7. In most cases the effect is minimal, although there are a few cases for which we could only infer an upper limit on the intrinsic velocity dispersion (marked with orange upper limits). We also note that the previous result from Übler et al. (2023) on an AGN at z = 5.55 (blue square) was obtained via spatially resolved NIRSpec IFS and in that case it is directly verified that the small amount rotation has a negligible contribution to the velocity dispersion (and it is accounted for). Summarising, the result of the M_BH − σ relation at z ∼ 4–7 being mostly consistent with the local one remains valid even if considering only the intrinsic velocity dispersion, although there might be a few outliers.

8.2. The black hole–dynamical mass relation

Given that the velocity dispersion is also linked to the dynamical mass through the virial theorem, it is interesting to explore also the relation between BH mass and dynamical mass.

To infer the latter we need the information on the radius of the galaxy and on its mass profile (probed by the Sérsic index), which were inferred through the morphological analysis in Sect. 3.4. We then estimate the dynamical mass by following the same approach as in Übler et al. (2023) through the equation

$$\begin{array}{c}\hfill M_{\mathrm{dyn}}=K(n)K(q)\frac{{\sigma }^2{R}_e}{G},\end{array}$$(2)

where K(n) = 8.87 − 0.831n + 0.0241n² with Sérsic index n, following Cappellari et al. (2006), K(q) = [0.87 + 0.38e^{−3.71(1 − q)}]², with axis ratio q following van der Wel et al. (2022), and R_e is the effective radius.

We note that, although the broadening associated with unresolved rotation is likely not contributing significantly to the observed velocity dispersion, the advantage of the dynamical mass is that Eq. (2) does require σ to include both intrinsic and unresolved rotation components (Cappellari et al. 2013). It therefore does not suffer from the potential issues discussed in the previous section. However, it is also true that, given the uncertainties on the several quantities involved in the determination of the dynamical mass (and in particular the radius), this quantity is more uncertain than the simple velocity dispersion measurement.

The resulting dynamical masses are reported in Table 3, limited to those cases for which we have enough information to derive this quantity. In most cases the dynamical mass is larger than the stellar mass. In some cases the two are consistent with each other within the uncertainties. However, in some cases the dynamical mass is significantly larger than the stellar mass, even by more than one order of magnitude. This is not totally surprising, as it is known that high redshift galaxies (already at z ∼ 1–4) can have gas fraction even higher than 90%, especially in the low-mass regime explored here (Tacconi et al. 2020; Santini et al. 2014; Scoville et al. 2017; Liu et al. 2019; Dessauges-Zavadsky et al. 2020; Zhang et al. 2021). The specific case of ID 11836 is remarkable, as the inferred dynamical mass is about two orders of magnitude higher than the stellar mass. It might be that in this case the stellar mass is significantly underestimated owing the difficulty to account for the AGN contribution, and indeed not including the power-law component gives a stellar mass that is ∼15 times higher; however, it is also true that previous studies at z < 4 have found some galaxies with inferred gas masses that are nearly two order of magnitudes higher than the stellar masses (Scoville et al. 2017; Zhang et al. 2021; Liu et al. 2019).

ID 954 is also another extreme case in which the inferred dynamical mass is about one order of magnitude lower than the stellar mass. This is the most luminous AGN in our sample, and it is likely that the spectral decomposition struggled to disentangle the stellar component.

The bottom right panel of Fig. 7 shows the BH mass as a function of dynamical mass, where the values obtained for the JADES broad-line AGNs are shown with golden large circles. In the local Universe there is shortage of proper equivalent measurements of the scaling relation in terms of dynamical masses. This was attempted by Beifiori et al. (2012); however, most of the BHs in their sample have upper limits on their masses and the sample has a limited mass range. We follow Übler et al. (2023) by taking, as a good approximation, the M_BH − M_bulge relation obtained by Kormendy & Ho (2013), where their measurement of M_bulge in ellipticals is close to a dynamical mass measurement. The values for local galaxies obtained by them are shown with small red squares and their best-fit relation is shown with a dash-dotted black line and the uncertainties and scatter with the shaded region. For completeness we also show the relations by Reines & Volonteri (2015) (solid black line) and Greene et al. (2020) (solid blue line), although, as is discussed, these are likely less appropriate for comparing with the dynamical mass. We also plot the estimates of M_BH and M_dyn for a sample of quasars at z ∼ 6.5 (grey triangles), for which the dynamical masses were inferred from the velocity dispersion of far-IR/submm transitions (typically [CII]158 μm) with ALMA (Izumi et al. 2019). We also report the previous measurement obtained with the NIRSpec IFU observation of an AGN at z = 5.5 by Übler et al. (2023).

The M_BH − M_dyn relation at high-z is, unfortunately, less populated than the M_BH − σ relation (due to the lack of information to derive the dynamical mass in a number of galaxies). The larger scatter of the M_BH − M_dyn relative to the M_BH − σ relation is probably partly a consequence of the additional uncertainties in deriving M_dyn. However, the interesting result is that the broad-line AGNs found in JADES are not strongly offset from the local relation, in contrast with the case of the M_BH − M_star relation. Broad-line AGNs at z > 4 are generally scattered around the local relation. There is a tendency for more AGNs being offset above the local relation, but not by more than two times the local dispersion, and consistent with the deviations seen in other local galaxies. Part of the better agreement may be a consequence of using here the Kormendy & Ho (2013) relation (which, as was discussed above, is more adequate in this case), instead of the Reines & Volonteri (2015) relation. However, even considering only the Kormendy & Ho (2013) local relation, it remains true that the offset is much larger on the M_BH − M_star than on the M_BH − M_dyn relation.

More accurate measurements are required, possibly with IFS spectroscopy. However, based on the current results, at face value we can say that the M_BH − M_dyn is generally consistent with the same relation in local galaxies.

Overall, these findings indicate that, at these early epochs, the BH mass follows well the mass assembly of the host galaxy, in a similar way as local galaxies. However, at these early epochs most of the mass is still in gas, which has not yet been converted to stars, and therefore explaining the large offset in the M_BH − M_star relation.

9. High-z active galactic nuclei are elusive in standard diagnostic diagrams

For AGNs at z < 7 the JADES spectra cover both [OIII] and Hβ as well as [NII] and Hα. Therefore, from the narrow components of these lines it is possible to locate these AGNs on the so-called BPT diagnostic diagrams (Baldwin et al. 1981; Veilleux & Osterbrock 1987). Generally in these high-z AGNs the [NII] doublet is very faint and in many cases undetected, so we can often only set an upper limit on the [NII]/Hα ratio. The resulting distribution of the JADES AGNs on the [NII]-BPT diagram is shown with golden circles in Fig. 8, in which the local distribution of galaxies from the SDSS survey is shown with shaded contours (where the lowest contour includes 99% of the local galaxies). The solid and dashed lines indicate the demarcation between AGN and star-forming galaxies provided by Kewley et al. (2001) and Kauffmann et al. (2003), respectively, for local galaxies.

Fig. 8. [NII]-BPT diagram. Contours show the distribution of galaxies in the local Universe from SDSS (DR7), showing the AGN sequence (right) and the star-forming sequence (left). The solid and dashed lines are the corresponding dividing lines from Kewley et al. (2001) and Kauffmann et al. (2003), respectively. The JADES broad-line AGNs presented in this paper are shown with golden circles, while blue diamonds show the values obtained by Übler et al. (2023) for an AGN at z = 5.5. Stars show predictions by AGN photoionisation models from Nakajima & Maiolino (2022) colour-coded by metallicity.

As has already been found by other studies based on JWST spectroscopic data (Kocevski et al. 2023; Übler et al. 2023; Harikane et al. 2023), AGNs at z > 4 are completely offset from the AGN locus in the local Universe, and are mostly overlapping with the region that is locally occupied by star-forming galaxies. Clearly, these diagnostic diagrams cannot discriminate AGNs from star-forming galaxies in these early systems.

10. Very low-metallicity host galaxies

One possible explanation for the offset in the BPT diagram is that the narrow emission lines are not primarily associated with the AGN NLR, but actually dominated by star formation in the host galaxy. Recently, Maiolino et al. (2024c) have suggested that the covering factor of the BLR clouds in these JWST-discovered AGNs could be very high, and hence leave few ionising photons escaping to produce an NLR on larger scales. However, we know that in at least a few cases the narrow lines are certainly dominated by the AGN. For instance, in the broad-line AGN GS-3073, explored in detail and with a very high S/N by Übler et al. (2023), the narrow lines also have the clear detection of strong HeII4686 (with high Equivalent Width, EW, typical of AGNs), and other high ionisation lines, such as [ArIV] and coronal lines, which are typically tracing the NLR of AGNs. Despite this, GS-3073 is clearly completely offset from the AGN local locus on the BPT diagram and is instead located on the local SF sequence (blue diamonds in Fig. 8).

Another likely possibility, also pointed out by previous works (Kocevski et al. 2023; Übler et al. 2023; Harikane et al. 2023), is that the NLR of high-z AGNs is characterised by low metallicities. This is indicated by the starred symbols in Fig. 8, which are the result of the photoionisation models for the NLR of AGNs obtained by Nakajima & Maiolino (2022) and colour-coded by metallicity. Clearly, as the metallicity of the NLR decreases from super-solar to sub-solar, the expected location on the BPT diagram shifts from the local AGN locus to overlap the local star-forming locus, and beyond. Most of the AGNs on the local star-forming sequence can be explained in terms of a NLR characterised by a metallicity of about 0.2 Z_⊙, although there is large dispersion in the models.

It is interesting that two of the AGNs in JADES are even below the local star formation sequence. In these cases, the metallicity of the NLR is probably particularly low, less than 0.1 Z_⊙. These targets are interesting candidates for follow-up with IFS observation to investigate their very low-metallicity environment.

11. Fraction of active galactic nuclei at high redshift

Previous JWST studies assessing the fraction of type 1 AGNs have reached different conclusions. Matthee et al. (2024) used NIRCam slitless spectroscopy to assess the fraction of broad-line AGNs in the EIGER and FRESCO surveys, and infer that AGNs in their sample are less than 1% of star-forming galaxies at z ∼ 5. It is important to consider that, due to the limited sensitivity of the slitless surveys used by them, they probe Hα broad-line luminosities higher than 8 × 10⁴² erg/s and resulting in AGN bolometric luminosities L_bol > 5 × 10⁴⁴ erg/s. Harikane et al. (2023) uses slit spectroscopy from the CEERS, ERO and GLASS surveys to search for AGNs, and estimate that 5% of galaxies at 4 < z < 7 host a broad-line AGN. Their slit spectra allow them to probe Hα broad lines with luminosities down to 10⁴¹ erg/s and infer AGN bolometric luminosities down to L_bol > 5 × 10⁴³ erg/s. Therefore, the higher fraction of broad-line AGNs inferred by Harikane et al. (2023) is likely a consequence of the lower luminosity range probed by them.

In the case of JADES, assessing the fraction of broad-line AGNs based on the current JADES NIRSpec spectra is not simple, as the selection function of the NIRSpec targets is quite complex, and varying in different tiers. In Deep/HST only a fraction of the targets was selected based on NIRCam data. Moreover, as is discussed in Sect. 2, Deep/HST, Medium/HST and some of the other early JADES observations have likely been biased against AGNs. These issues may be responsible for the fact that only two broad-line AGNs have been found in Deep/HST and two in Medium/HST GOODS-N.

The Medium/JWST GOODS-N tier is the first JADES tier mostly based on NIRCam-selected targets and in which high-z AGNs may have not been discarded because of their peculiar colours, and hence not biased against AGNs. The fraction of broad-line AGNs found in this tier is probably more representative of the population of AGNs at high redshift. However, in the case of Medium/JWST a few objects (specifically 73 488, 77 652 and 61 888) were targeted just because of their compact morphology and peculiar colours (resembling the type 1 AGNs, such as the one identified by Übler et al. 2023) suggesting their AGN nature (Sect. 2), and were indeed confirmed as such in our analysis; therefore, in the case of this sub-sample the selection has obviously been biased in favour of AGNs.

Aware of all these caveats, we attempt an assessment of the fraction of AGNs in the Medium/JWST GOODS-N tier in the redshift range 4 < z < 6, in which the statistics of targeted galaxies with confirmed spectroscopic redshift is high enough. The fraction of broad-line AGNs (with L_AGN > 10⁴⁴ erg/s, i.e. the luminosity range to which we are sensitive) in this subsample is 11%. If we exclude the galaxies that were specifically selected because of their AGN-like properties in imaging, then the fraction drops to 7.5%. However, removing the AGNs targeted because of their imaging and photometric properties is actually biasing the sample against AGNs (as in the early JADES tiers), as their properties are actually in the range of the JADES selection function. Therefore, actual fraction of galaxies at 4 < z < 6 hosting broad-line AGNs with L_AGN > 5 × 10⁴³ erg/s must be between the two estimates given above; in other words, between 7.5% and 11%. This is higher than what inferred by Harikane et al. (2023), who probed a similar sensitivity range. However, the statistics is still modest in both studies (nine targets in this redshift range in both studies). Moreover, the selection functions in the two samples are different and the selection criteria in GLASS and CEERS surveys may have penalised the selection of type 1 AGNs.

As is discussed in Sect. 5, the AGNs in our sample have significantly lower dust extinction than the broad-line AGNs found by Matthee et al. (2024) via slitless spectroscopy, and this may be a consequence of the JADES targets selection function (which has various priority classes that are selected also based on the UV luminosity) that has likely penalised reddened AGNs. So the actual fraction of broad-line AGNs, including the significantly reddened population in our luminosity range, is certainly higher. It is difficult to assess accurately this fraction, as it would require a dedicated, deep MOS spectroscopic survey purely selected based on the red photometry (e.g. F444W).

Finally, our estimated fraction of AGNs refers only to the class of type 1, while the fraction of obscured, type 2 AGNs is certainly higher. In the local Universe, type 2 AGNs are about a factor of four more numerous than type 1 AGNs (e.g. Maiolino & Rieke 1995). At high redshift, the fraction of type 2 AGNs has been more difficult to assess, due to sensitivity issues (e.g. Merloni et al. 2014; Maiolino et al. 2007; Netzer et al. 2016). These aspects are explored more extensively, by using the JADES data, in a separate paper (Scholtz et al. 2023).

It is important to compare our results with the expectations of some models and simulations. In particular, using the semi-analytical model CAT, Trinca et al. (2023) predict that, at the magnitude limit probed by the Medium/JADES tier, we should detect a few dozen type 1 AGNs in the redshift range of 5 < z < 6. This is not too far from the seven that we have detected in the same redshift range, taking into account that we have only used about 1/3 of the Medium/JWST tier (the GOODS-N component) and that the Medium/HST tier may be biased against AGNs, as was discussed above. Moreover, we may have still missed a significant fraction of AGNs, because too faint relative to their host galaxy, an issue that has been extensively discussed in Volonteri et al. (2023) and Schneider et al. (2023). It is interesting to note that Trinca et al. (2023) expect that most of these AGNs should host BHs in the mass range 10⁶ − 10⁸ M_⊙, as found by us. According to their model most of these BH are formed out of ‘heavy’ seeds; that is, from DCBHs.

12. Contribution of active galactic nuclei to the ultraviolet luminosity function

We can explore the contribution of the galaxies that host broad-line AGNs to the UV luminosity function. Assuming that within the 4 < z < 6 redshift bin, selection effects within each UV luminosity bin are secondary (i.e. star-forming galaxies and AGNs in a given UV luminosity bin do not have significantly different probability of being selected for spectroscopy), then the contribution of galaxies hosting broad-line AGNs in each bin of the UV luminosity function can be inferred from the fraction of AGN in that bin in our spectroscopically targeted sample. As a reference, we take the functional form provided by Bouwens et al. (2021) to describe the UV luminosity function of galaxies at z ∼ 5. As was discussed above, we have the ambiguity of excluding or including the three AGNs specifically targeted because of their AGN appearance in the NIRCam images. We take the mid-point of the two extreme cases and take the deviation from each of the two cases as a contribution of the errorbars, in quadrature to the Poissonian noise and cosmic variance (estimated following Somerville et al. 2004, where the bias factor is function of redshift and mass).

The UV luminosity function of broad-line AGNs with L_AGN > 5 × 10⁴³ erg/s and 4 < z < 6 inferred by us from the JADES survey is shown with golden circles in Fig. 9 and reported in Table 4 (in the last column of the same table we also provide the minimum density of broad-line AGNs, assuming the most conservative and extreme case, that the AGNs identified in the current spectroscopic sample are the only AGNs in the entire volume surveyed; in other words, that there are no other AGNs missed among the vast majority of galaxies that have not been targeted spectroscopically). To obtain a reasonable statistics in UV luminosity bins, we only adopt two bins, one centred at M_UV = −18 and a second one centred at M_UV = −19.5. In the lower luminosity bin broad-line AGNs and their hosts contribute 12% to the luminosity function of galaxies, while in the higher luminosity bin the contribution is 8%. Certainly these contributions are expected to be higher if one could probe lower luminosity AGNs.

Fig. 9. UV luminosity function of galaxies at z = 5 from Bouwens et al. (2021) (solid blue line) and inferred contribution of galaxies hosting broad-line (type 1) AGNs with Lbol > 1044 erg/s, inferred by us from the JADES survey (golden circles). Results from other surveys are also reported, as is indicated in the legend. The dashed blue line shows the galaxy luminosity function scaled downward by a factor of ten and fitting the JADES points. The orange-shaded region shows the range of possible extrapolated luminosity functions for QSOs from Niida et al. (2020).

We do not attempt to fit a functional form of the luminosity function to these two data points, as the statistics is still too low and errorbars still too large, which would result in different functional forms being largely unconstrained. However, we can reproduce the values for broad-line AGNs by simply scaling down the Bouwens et al. (2021) UV luminosity function of galaxies by a factor of ten, as is illustrated by the dashed blue line.

In Fig. 9 we also compare our finding with results from other JWST surveys: Harikane et al. (2023) blue squares; Matthee et al. (2024) green triangles; Kocevski et al. (2023) violet pentagon. The red triangles show the luminosity function inferred by Giallongo et al. (2019) based on X-ray surveys. We also show the range of extrapolations of the quasar luminosity function at z ∼ 5 as inferred by Niida et al. (2020).

Our results are, within uncertainties, consistent with the finding of Harikane et al. (2023), not surprisingly given that they probe a similar range of AGN luminosities. The luminosity function inferred by Matthee et al. (2024) is much lower, which is likely a consequence of the higher luminosities probed by them. The low value inferred by Kocevski et al. (2023), which uses the same CEERS spectra as Harikane et al. (2023), is probably a consequence of the low statistics in that early study. It is interesting that our estimated density of AGN is higher than inferred by Giallongo et al. (2019) by using some of the deepest X-ray data, indicating that JWST is currently being much more effective in finding AGNs than current X-ray surveys and indicating that future X-ray missions are necessary to find and characterise this population of high-z AGNs and, most importantly, find the obscured counterparts.

We finally note that the range of possible extrapolations of the QSO luminosity function at z ∼ 5 (orange shaded region) is clearly below the AGN luminosity function estimated by us, indicating that we are probing a quite different population relative to luminous quasars, and not simply their low-luminosity tail.

Given that the yield of ionising photons (ξ_ion) by AGN is larger than for star-forming galaxies, these findings imply that the contribution to re-ionisation of galaxies hosting broad-line AGNs (with L_AGN > 10⁴⁴ erg/s) can potentially be significant.

We note that this result does not necessarily imply that AGNs contribute substantially to the re-ionisation of the Universe, as the UV luminosity that we measure is the sum of the contribution from the AGN and host galaxy. In order to obtain the specific contribution from BH accretion we should disentangle in each of the selected galaxies the contribution of the (dust-reddened) accretion disc from the light emitted by the young stellar population. The decomposition attempted in Sect. 5 is appropriate to infer the properties of the stellar population, however it is inappropriate to use it to extrapolate the AGN contribution to the extreme UV. A more detailed modelling (also involving the nebular lines) is needed for this goal. However, if most of the UV emission turns out to be dominated by the AGN, then (given their large escape fraction) these could potentially contribute to a large fraction of the photon budget required for the re-ionisation of the Universe (Madau et al. 2024).

13. Summary and conclusions

We have used three tiers of the JADES NIRSpec survey, specifically Deep/HST (in GOODS-S), Medium/HST, and Medium/JWST (in GOODS-N), to search for broad-line AGNs at z > 4. The combination of depth and the use of dispersers providing three different resolutions has enabled us to find this class of AGNs more efficiently and has allowed us to explore different regimes relative to previous studies.

In addition to the previously discovered GN-z11 at z = 10.6, we have identified twelve new broad-line AGNs at z > 4. In these cases, the Hα line emission shows a broad component (in addition to a narrow component tracing the ISM in the host galaxy) that does not have a counterpart in [OIII]5007, and hence cannot be ascribed to outflowing gas and is most likely tracing the broad-line region of an AGN. Our analysis of the 13 broad-line AGNs reveals the following findings:

• In three cases, the Hα profile requires an additional intermediate-width component (FWHM ∼ 400–700 km/s). We interpret this additional component as tracing a secondary accreting BH with smaller mass, in the same galaxy, which will likely merge with the larger BH. The finding in the other two galaxies of broad lines that are significantly shifted relative to the narrow component may indicate that these are also BHs in the process of merging, but in which the more massive BH is not accreting at the time of observation; alternatively, these could be BHs recoiled from a recent merger. However, we cannot exclude that in these cases the complex broad-line profile is due to a complex geometry of the BLR. Follow-up observations with IFU spectroscopy may help to further assess the merging nature of these objects (as it has already successfully happened for other samples, Übler et al. 2024).

• By using local virial relations, we have inferred BH masses that are in the range between 4 × 10⁵ M_⊙ and 8 × 10⁷ M_⊙. Interestingly, the lowest-mass BH is in the regime of DCBHs, which is one of the favoured scenarios for the heavy seeds of supermassive BHs. This does not imply that this is a DCBH, as it may have formed from another kind of seed and gained its mass through various evolutionary paths; however, this result shows that we are now capable of probing this regime that is potentially populated by DCBHs. More detections, and hence more statistics, of BHs in this mass range will help to test different seeding scenarios, especially in merging systems.

• Although estimating the intrinsic bolometric luminosity is difficult, most BHs in our sample seem to be accreting at sub-Eddington rates, mostly with L/L_Edd < 0.5. However, small BHs, with masses of M_BH ∼ 10⁶ M_⊙, tend to accrete at Eddington or super-Eddington rates. This might be a consequence of selection effects; small BHs become detectable only when they are accreting at a very high rate. However, this finding also provides support to scenarios that envisage phases of super-Eddington accretion in the early phases of BH formation.

• We have found that BHs at 4 < z < 11 are over-massive relative to their host galaxies, when compared to the local M_BH − M_star relation. We even find cases that are approaching M_BH ∼ M_star. While selection effects might be partially responsible for this finding, the result may indicate that at early epochs BHs may form and grow faster than the stellar population in their host galaxies. We also note that high M_BH/M_star ratios at high redshift are an expectation of models that envisage super-Eddington accretion at early epochs and/or heavy seeds (i.e. DCBHs).

• The high-resolution spectra allowed us to estimate the velocity dispersion in the host galaxy. We find that the M_BH − σ relation of AGNs at z > 4 is generally consistent with the local relation, with only a few exceptions. In the local Universe, the M_BH − σ relation (in which M_BH scales as ∼σ⁴) is the tightest of all BH scaling relations with the host galaxy properties, and therefore considered the most fundamental (while other relations are possibly an indirect by-product). Our finding that high-z BHs follow the same relation confirms that the M_BH − σ is more fundamental than other scaling relations, and that it is also universal; that is, it holds at least out to z ∼ 6.

• We have also attempted to estimate the dynamical masses, based on measurements of the radius of the host galaxies. The JADES broad-line AGNs are broadly consistent with the local M_BH − M_dyn relation (within its 2σ scatter), although in this case the uncertainties are larger than for the M_BH − σ relation.

• The large scatter and strong deviation of the M_BH − M_star relation relative to local galaxies, together with the fact that the M_BH − σ and M_BH − M_dyn relations at z > 4 are instead consistent with the local relations, suggest that BH formation is little connected to the formation of stars in the host galaxy, while it is tightly connected to the mass assembly history of the central spheroidal component.

• We have found that the location of the narrow components of these high-z AGNs on the [NII]-BPT diagram is completely offset from the local locus of AGNs, while overlapping with the local star-forming sequence. This confirms that some of the standard nebular optical diagnostics for identifying AGNs are ineffective at high redshift.

• We have shown that the offset of high-z AGNs on the BPT diagram is consistent with the fact that these systems, and their narrow-line region, are metal-poor (typically Z ∼ 0.2 Z_⊙). We even find two targets that are located below the local star-forming sequence; these are likely AGNs whose NLR have very low metallicities, below 0.1 Z_⊙.

• We estimate that the fraction of broad-line AGNs with L_bol > 10⁴⁴ in galaxies at 4 < z < 6 is about 10%.

• The luminosity function of galaxies hosting broad-line AGNs with L_bol > 10⁴⁴ at 4 < z < 6 is consistent with the luminosity function of galaxies, in the same redshift range, scaled down by a factor of ten.

• The contribution of galaxies hosting broad-line AGNs to the re-ionisation of the Universe is larger than 10%. Establishing the specific contribution of AGNs (without their host galaxies) requires disentangling the AGN and stellar contributions in these galaxies more precisely, as well as larger statistics.

Acknowledgments

We thank Marta Volonteri, Raffaella Schneider, Alessandro Trinca, Debora Sijiacki, Martin Haehnelt, Jake Bennett, Sophie Koudmani and the anonymous referee for useful comments. FDE, JS, RM, TL, WB acknowledge support by the Science and Technology Facilities Council (STFC), by the ERC through Advanced Grant 695671 “QUENCH”, and by the UKRI Frontier Research grant RISEandFALL. RM also acknowledges funding from a research professorship from the Royal Society. ECL acknowledges support of an STFC Webb Fellowship (ST/W001438/1) S.C & G.V. acknowledge support by European Union’s HE ERC Starting Grant No. 101040227 - WINGS. HÜ gratefully acknowledges support by the Isaac Newton Trust and by the Kavli Foundation through a Newton-Kavli Junior Fellowship. S.A. and M.P. acknowledge support from Grant PID2021-127718NB-I00 funded by the Spanish Ministry of Science and Innovation/State Agency of Research (MICIN/AEI/ 10.13039/501100011033). AJB & GCJ acknowledge funding from the “FirstGalaxies” Advanced Grant from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 789056) E.E. and DJE are supported as a Simons Investigator and by JWST/NIRCam contract to the University of Arizona, NAS5-02015 BER acknowledges support from the NIRCam Science Team contract to the University of Arizona, NAS5-02015. M.P. acknowledges support from the research project PID2021-127718NB-I00 of the Spanish Ministry of Science and Innovation/State Agency of Research (MICIN/AEI/ 10.13039/501100011033), and the Programa Atraccion de Talento de la Comunidad de Madrid via grant 2018-T2/TIC-11715 The research of CCW is supported by NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. The authors acknowledge use of the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315.

References

Appendix A: [OIII]+Hβ spectra

Figure A.1 shows a zoomed-in view of the spectra around Hβ and [OIII], primarily to illustrate that the broad component is not seen in the [OIII] line, despite being brighter than Hα, and the S/N being very high. In two cases, namely ID 8083 and ID 53757, [OIII] is not covered by grating spectra (because in the detectors gap); however, the Hα line profile is so broad and symmetric that it is very unlikely that is due to an outflow. Moreover, in these cases the flux of the broad Hα is higher than the narrow component, if the former was due to an outflow then it would imply that the mass of ionised gas in outflow is larger than in the ISM of the whole galaxy, which is very unlikely.

We generally show only the medium-resolution grating, with the exception of the candidate dual AGN ID 10013704, for which we show both the low- and high-resolution spectra, to illustrate that neither of the two BLR components is present.

The fainter Hβ generally does not show a broad component. The broad component seen in Hα is not expected to be seen in the Hβ, especially with the modest dust extinction inferred towards the central region, as is discussed in the text. The only exception is ID 954 (the most luminous AGN in the sample) for which the BLR is seen also for Hβ.

Fig. A.1. Spectra around [OIII] and Hβ for ten out of the twelve new JADES galaxies with evidence for a broad component of Hα. The solid blue line shows the spectrum (not continuum subtracted) along with the errors (light blue shaded area). The dashed red line shows the total multi-component fit; the dashed blue and purple lines show the narrow and broad components of Hβ, respectively (the broad component of Hβ is detected only in ID 954), while the dashed green lines show the components the [OIII] doublets. The main purpose of these spectra is to illustrate that there is no [OIII] counterpart of the Hα broad component. For simplicity, for all galaxies in the sample, we only show the medium-resolution grating spectrum, with the exception of the candidate dual AGN ID 10013704, for which we also show the high-resolution spectrum.

Appendix B: [OIII] profile of the dual broad-line region in ID 10013704

In this appendix we perform a closer analysis of the [OIII] profile of ID 10013704. The left panels of Fig. B.1 show a version of the medium (top) and high (bottom) resolution spectra zoomed around the [OIII]5007 line. We overplot the three Gaussian components used to reproduce the Hα profile, rescaled so that the Hα narrow component flux matches the [OIII] flux. Clearly, none of the two broad components has a counterpart in the [OIII] profile.

The medium-resolution [OIII] profile has the hint of a redshifted component, but mostly consistent with the noise, not seen in the high-resolution spectrum, and anyway not matching the profile of either of the two Hα broad components. This is more clearly seen in the right panels, where we show the same spectra after subtracting the narrow [OIII] component. The marginal redshifted residual is consistent with other features in the noise and not matching either of the two broad Hα components. If confirmed with higher-resolution data, this faint component may be associated with the remnant of the merging galaxy that hosted the smaller BH.

Fig. B.1. Analysis of the [OIII]5007 line in the spectra of ID 100013704 excluding the presence of a broad component. Left panels: medium (top) and high (bottom) resolution spectra of ID 10013704 (candidate to host a dual BLR), zoomed around the [OIII]5007 line. The dashed lines indicate the Gaussian components used to fit the Hα profile, rescaled so that the narrow component matches the flux of [OIII]. The right panels show the same spectra, where the [OIII] (narrow) profile has been subtracted. Clearly, there is no [OIII] counterpart of the Hα broad components.

Appendix C: Attempt to fit a double power law to the Hα profile of ID 73488

In this appendix, we show that a double power-law profile (sometimes used to fit the broad components of high-z quasars, e.g. Nagao et al. 2006), cannot reproduce the broad Hα profile of ID 73488. The fit (in addition to the standard narrow component) is shown with dashed lines in Fig. C.1, both for the low (left) and high (right) resolution spectra. The double power law fails to properly fit the broad component of Hα, as is highlighted by the strong systemic residuals, which are shown in the bottom panels.

The ΔBIC of this fit, relative to the single Gaussian fit, is 273, much higher than the ΔBIC_1B2B = 494 obtained from the BIC difference between the fit with one Gaussian and two Gaussians. This indicates that the two-Gaussians models describes the data much better than the double power law model with a very high level of confidence.

Fig. C.1. Medium (left) and high (right) resolution spectra of ID 73488, zoomed around Hα with the fit overlaid using a narrow (Gaussian) component (dashed blue line) and a double power-law profile (dashed purple line) to describe the broad component. The total fit is shown with the dashed red line. The double power-law profile fails to properly fit the broad component, as is highlighted by the strong residuals, which are shown in the bottom panels.

Appendix D: Beagle spectral fits

Fig. D shows the low-resolution prism spectra of the 12 AGNs in our sample (GN-z11 is discussed separately in Maiolino et al. 2024a), with the Beagle best fit, adopting the composite model including both a stellar population and a reddened AGN power-law. The nebular emission lines are masked, hence the fit is applied only to the continuum sections of the observed spectrum, shown in red. The Beagle best fit is shown with the blue line.

Fig. D.1. Low-resolution prism spectra (dark and light red line) of the targets in our sample, with overlaid the Beagle best fit to the continuum (dark blue line).

All Tables

All Figures

Fig. 1. Medium-resolution spectra of the 12 new JADES galaxies with evidence for a broad component of Hα ascribed to the BLR of an AGN. We show a zoom around Hα. The solid blue line shows the spectrum along with the errors (light blue shaded area). The dashed red line shows the total multi-component fit; the dashed blue and purple lines show the narrow and broad components of Hα, respectively, while the dashed green lines show the components fitting the [NII] doublet (often undetected).

In the text

	Fig. 2. NIRCam images of the selected broad-line AGNs. The following false-colour coding was adopted: blue – F115W; green – F200W; and red – F444W. Each thumbnail is 3″ × 3″ in size.
In the text

Fig. 3. Spectra around Hα of ID 10013704, the highest-redshift AGN showing an indication for a dual BLR. The top and bottom panels show the medium- and high-resolution spectra, respectively. The line coding is the same as in Fig. 1, but in this case the dashed violet line shows the intermediate broad component that is needed to properly reproduce the observed profile. The central and rightmost panels show the spectrum from which the narrow components, as well as one of the two broad components have been removed, to better highlight the significance of the other broad component.

In the text

	Fig. 4. As Fig. 3 but for the additional candidate dual BLR AGN ID 73488.
In the text

	Fig. 5. As Fig. 3 but for the additional candidate dual BLR AGN ID 53757.
In the text

Fig. 6. Distribution of BH masses and Eddington ratios (L/LEdd) for the broad-line AGNs in JADES (large golden circles). The dashed orange lines connect candidate dual AGNs. We also show the results from other JWST surveys with blue symbols (see legend, only detections at > 3σ are shown). Note that the apparent anticorrelation is probably spurious, as the BH mass is at the denominator of the Y-axis quantity. The plot has simply the purpose of visually illustrating the distribution of the two quantities.

In the text

Fig. 7. BH mass as a function of host galaxy properties, specifically: Stellar mass (top left), velocity dispersion (bottom left), and dynamical mass (bottom right). The JADES results presented in this work are shown with large golden circles. Dashed orange vertical lines connect candidate merging BHs. Blue symbols indicate measurements from other JWST surveys at high-z, as is indicated in the legend (only detections at > 3σ are shown). Grey triangles are measurements of high-z QSOs using ALMA data. Small red symbols show the distribution of local galaxies, as is indicated in the legend; straight lines show the local relation fits (with shaded regions providing the scatter and slope uncertainty), using the samples that are best matching our high-z systems, in each panel (see text for details). In the top left panel, small orange circles show the maximum stellar masses estimated for a few JADES AGNs (see text) and the dotted violet lines show constant values of the MBH/Mstar ratio. In the bottom left panel, small squares indicate the effect of correcting for rotation velocity broadening within the slit.

In the text

Fig. 8. [NII]-BPT diagram. Contours show the distribution of galaxies in the local Universe from SDSS (DR7), showing the AGN sequence (right) and the star-forming sequence (left). The solid and dashed lines are the corresponding dividing lines from Kewley et al. (2001) and Kauffmann et al. (2003), respectively. The JADES broad-line AGNs presented in this paper are shown with golden circles, while blue diamonds show the values obtained by Übler et al. (2023) for an AGN at z = 5.5. Stars show predictions by AGN photoionisation models from Nakajima & Maiolino (2022) colour-coded by metallicity.

In the text

Fig. 9. UV luminosity function of galaxies at z = 5 from Bouwens et al. (2021) (solid blue line) and inferred contribution of galaxies hosting broad-line (type 1) AGNs with Lbol > 1044 erg/s, inferred by us from the JADES survey (golden circles). Results from other surveys are also reported, as is indicated in the legend. The dashed blue line shows the galaxy luminosity function scaled downward by a factor of ten and fitting the JADES points. The orange-shaded region shows the range of possible extrapolated luminosity functions for QSOs from Niida et al. (2020).

In the text

Fig. A.1. Spectra around [OIII] and Hβ for ten out of the twelve new JADES galaxies with evidence for a broad component of Hα. The solid blue line shows the spectrum (not continuum subtracted) along with the errors (light blue shaded area). The dashed red line shows the total multi-component fit; the dashed blue and purple lines show the narrow and broad components of Hβ, respectively (the broad component of Hβ is detected only in ID 954), while the dashed green lines show the components the [OIII] doublets. The main purpose of these spectra is to illustrate that there is no [OIII] counterpart of the Hα broad component. For simplicity, for all galaxies in the sample, we only show the medium-resolution grating spectrum, with the exception of the candidate dual AGN ID 10013704, for which we also show the high-resolution spectrum.

In the text

Fig. B.1. Analysis of the [OIII]5007 line in the spectra of ID 100013704 excluding the presence of a broad component. Left panels: medium (top) and high (bottom) resolution spectra of ID 10013704 (candidate to host a dual BLR), zoomed around the [OIII]5007 line. The dashed lines indicate the Gaussian components used to fit the Hα profile, rescaled so that the narrow component matches the flux of [OIII]. The right panels show the same spectra, where the [OIII] (narrow) profile has been subtracted. Clearly, there is no [OIII] counterpart of the Hα broad components.

In the text

Fig. C.1. Medium (left) and high (right) resolution spectra of ID 73488, zoomed around Hα with the fit overlaid using a narrow (Gaussian) component (dashed blue line) and a double power-law profile (dashed purple line) to describe the broad component. The total fit is shown with the dashed red line. The double power-law profile fails to properly fit the broad component, as is highlighted by the strong residuals, which are shown in the bottom panels.

In the text

	Fig. D.1. Low-resolution prism spectra (dark and light red line) of the targets in our sample, with overlaid the Beagle best fit to the continuum (dark blue line).
In the text