Undermassive Host Galaxies of Five z ∼ 6 Luminous Quasars Detected with JWST

All our quasars were observed in three NIRCam filters chosen from F210M, F360M, F410M, F430M, and F480M (see Table 1), with exposure times ranging from 265.3 to 2623.8 s depending on the quasar and filter. We used Module B SUB400P (field of view = 125 × 125, shortwave (SW) band; 25'' × 25'', longwave (LW) band) to image the quasar field using either the RAPID or BRIGHT2 readout modes. We adopted 4 × 4 dithering patterns to improve the PSF sampling and mitigate cosmic rays and detector artifacts, using primary dither type INTRAMODULEBOX and subpixel dither type STANDARD. Reference PSFs were obtained for all quasars by observing nearby bright stars. Stars were observed using the same module, dither types, and filter configurations as their corresponding quasars, but all stellar observations used the RAPID readout mode to avoid saturation. All stars were observed for 530.7 s in the SW bands and 265.3 s in the LW bands.

Table 1. Measured Properties of the High-redshift Quasar Sample

Quasar	Redshift	Band	Quasar Flux	Measured Galaxy Flux	Upper Limit Galaxy Flux	Inferred Galaxy Mass
			(μJy)	(μJy)	(μJy)	log(M⊙)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
J0731+4459	5.01	F410M	194.0 ± 0.2	0.24 ± 0.21	4.32
		F430M	112.0 ± 0.2	0.85 ± 0.22	5.12	<10.7
		F210M	95.1 ± 0.3	0.78 ± 0.28	1.56
J1340+2813	5.36	F410M	173.2 ± 0.2	0.88 ± 0.18	5.28	<10.8
		F430M	192.7 ± 0.2	1.06 ± 0.23	6.36
		F210M	6.81 ± 0.04	0.12 ± 0.07	0.60
J2239+0207	6.25	F360M	7.43 ± 0.05	0.34 ± 0.05	1.36	${10.0}_{-0.3}^{+0.4}$
		F480M	11.32 ± 0.09	0.45 ± 0.08	1.32
		F210M	120.5 ± 0.2	0.39 ± 0.23	0.78
J1148+5251	6.40	F360M	117.1 ± 0.2	0.45 ± 0.18	2.70	<10.7
		F410M	103.2 ± 0.1	0.48 ± 0.11	1.92
		F210M	42.0 ± 0.1	0.12 ± 0.12	0.48
J1120+0641	7.09	F360M	53.2 ± 0.1	0.20 ± 0.07	0.80	<10.3
		F480M	48.3 ± 0.1	0.31 ± 0.13	0.62

Note. The measured quasar and host galaxy fluxes of our five high-z quasars. The measured galaxy fluxes in column (5) were obtained by normalizing our reference PSF to match the flux of the quasar at the center of the image before subtracting, therefore subtracting the flux to zero near the center of the host galaxy (see Section 3.1). The uncertainties in column (5) reflect purely statistical errors, without taking into account systematics. We then adjusted the normalization of the reference PSF to determine how much flux could be attributed to the host galaxy near the image center, and adopted the maximum flux with the subtracted image and radial profile consistent with a Sérsic profile rather than a PSF as the upper limit on the galaxy flux (column (6)). We then calculate masses (column (7)) by applying Equation (1) to the limits on the host galaxy flux.

Download table as:  ASCIITypeset image

We processed our NIRCam data using the JWST pipeline v1.9.6 (Bushouse et al. 2023), roughly following the procedures recommended in the STScI JWebbinars. ¹ The pipeline parameter reference file is registered in the JWST Calibration Reference Data System (CRDS) as jwst_1084.pmap. We added a custom step to Stage 2 of the pipeline to characterize and subtract 1/f noise, a striping pattern in the images caused by the detector readout, from each frame. We produced final mosaic images of each quasar and star in Stage 3 of the pipeline by aligning and stacking all frames obtained in Stage 2, and resampled the images to smaller pixel scales (00147 pixel⁻¹ in the SW filters; 00300 pixel⁻¹ in the LW filters) using the drizzling algorithm in the Resample step of the pipeline to improve the accuracy of the PSF subtraction. Finally, we used the Photutils Background2D function on the mosaiced and resampled images to perform a global background subtraction.

We also obtained NIRSpec IFU data of three quasars in our sample, with a 3'' × 3'' field of view composed of IFU elements of 01 × 01 on the sky. J0731+4459 and J1340+2813 were observed with the G235M/F170LP disperser/filter combination, offering wavelength coverage from 1.66 to 3.07 μm with a nominal spectral resolution ∼ 1000; J2239+0207 was observed in the PRISM mode with a wavelength coverage from 0.60 to 5.30 μm with a nominal resolution ∼ 100. All these observations were carried out with the SPARE-CYCLING dither type in LARGE size at points 1, 2, 3, and 4 to mitigate cosmic rays and detector artifacts during data reduction. The NRSIRS2 readout was selected, with a total integration time of 5894 s for J0731+4459, 7353 s for J1340+2813, and 8870 s for J2239+0207.

We processed the NIRSpec IFU data with JWST pipeline v1.12.0 (Bushouse et al. 2023; with CRDS pipeline parameter reference file jwst_1088.pmap) following the standard steps. This includes the first stage, Detector1Pipeline, to apply detector-level corrections (e.g., dark current and bias subtraction, persistence correction, and cosmic-ray removal) to the raw data of individual exposures and produce the uncalibrated 2D spectra; and Spec2Pipeline to assign a world coordinate system information to the data, apply a flat-field correction, make a flux calibration, and construct 3D data cubes from the 2D spectra obtained at each dither location. The third stage, Spec3Pipeline, combines the individual data cubes at different dither positions to produce the final merged data cube with outlier rejections, drizzling, etc., to remove any additional artifacts and improve spatial sampling. We found that the default pipeline did not always reject all obvious outliers; we therefore added an additional step to mask out the bad pixels in the 2D spectra manually after a careful visual inspection of individual cubes, and redid the final cube construction. Finally, we extracted the quasar spectrum using a circular aperture with a radius of 035. For the background subtraction, we placed the same circular aperture at random locations of the cube that do not contain real structures, computed the median background spectrum, and subtracted it from the quasar spectrum.

Our quasar observations have small fields of view with no bright stars, so building an empirical PSF from stars in the field (as in, e.g., Ding et al. 2023; Yue et al. 2023) was not a viable option. Instead, for each quasar we observed a nearby PSF star in the same bands immediately after observing the quasar (to reduce any temporal variation in the instrument's PSF between the images) and subtracted it from the quasar image to reveal any underlying host galaxy emission.

Our approach has the advantage of being immune to different wave front errors and PSF variation over the field of view. That is, the quasar and star are placed at the same location on the detector to minimize the effects from any spatial variation of the PSF across the detector. We additionally chose medium bands for all NIRCam observations to eliminate differences in the profiles of the quasar and stellar PSFs due to their different spectral shapes across the filter (the flat continuum of the quasar versus the declining Rayleigh–Jeans tail of the star). In a wide bandpass, the quasar's redder color in the filter could tend to broaden its PSF; if we then subtracted the narrower PSF of the star, we might observe a spurious galaxy detection.

However, we still considered the possibility of a discrepancy in the PSF shapes. We investigated the magnitude of any potential PSF mismatch by generating simulated stellar and quasar PSFs in each of our observed bands using WebbPSF. We used a Rayleigh–Jeans spectrum to generate the stellar PSFs, and a general quasar template (Lyu et al. 2016), shifted to the appropriate redshift, for the quasars. We constructed differential radial profiles of the stellar and quasar PSFs, and compute the normalized difference between the two as a function of radius. As in S23, we find no significant difference between the simulated quasar and stellar PSFs in any band. The normalized difference between the simulated stellar and quasar PSF at any given radius, $\tfrac{F(\mathrm{star})-F(\mathrm{QSO})}{F(\mathrm{star})}$, is significantly less than 1% in all bands and at almost all radii, except that normalized differences up to 5% occur very near the PSF core where the annulus contains fewer pixels. We therefore do not expect any significant residuals to be introduced in the PSF subtraction as a result of the different spectral shapes of the quasars and their PSF stars, particularly outside the PSF core in the region where we search for evidence of the host galaxy.

The process of PSF subtraction begins by aligning the background-subtracted star and quasar images in a given band. We perform this alignment by hand, iterating until minimal PSF artifacts are present in the subtracted images. Next, we normalize the stellar PSF image to match the flux of the quasar near the center of the image (within ∼015), where the quasar light is most dominant. As we iterate the alignment and scaling, we examine the azimuthally averaged radial profiles of the star and quasar, shown for all quasars in Figure 1. If the star and quasar are properly aligned and scaled, their profiles will exhibit identical shapes at small radii. Indeed, as shown in Figure 1, the profiles of the quasars and their PSF stars are virtually identical in all bands at small radii. Far from the center (r ≳ 04), deviations between the stellar and quasar PSFs are due to noise. The excellent agreement in the shape of the star and quasar radial profiles underscores the results from our WebbPSF simulations—despite their different spectral shapes, the quasars and their corresponding PSF stars exhibit nearly identical PSFs near the image centers. Any deviations between the PSF shapes of a quasar and star at slightly greater radii (∼02–04 from the PSF center) are likely caused by extended emission from the host galaxy, rather than a PSF mismatch.

Zoom In Zoom Out Reset image size

Once the stellar reference PSF is aligned and scaled to the quasar PSF, we subtract the stellar PSF and examine the residuals for evidence of extended emission. The PSF-subtracted images of all quasars are shown in Figure 2 in order of ascending redshift. The cleanliness of the images, with excess noise under the core of the PSF but few other residuals, demonstrates the quality of the subtractions. Only the F410M image of J0731+4459 displays a higher degree of PSF mismatch than the other images. Most of the images do not display evidence of extended emission, with the possible exceptions of J0731+4459 in F430M, J1340+2813 in F410M and F430M, and—as discussed in S23—J2239+0207 in F360M and F480M.

Zoom In Zoom Out Reset image size

We performed tests on the images of J2239+0207 to compare the results of our manual alignment technique to computerized alignment. We used the photutils.centroid.quadratic function to determine the centroids of the quasar and stellar images, and calculated the shift necessary to align the reference PSF to the quasar image. The differences between the shifts obtained from manual and computerized alignment were of order 0.1–0.2 pixels. We then subtracted the PSF from both the manually and automatically aligned images, and measured the resulting flux in a 04 radius aperture. The PSF-subtracted fluxes from manual versus computerized alignment were nearly identical, differing at the 0.001%, or 10 pJy, level. Our manual alignment technique is therefore reliable, and produces results identical to those from computerized alignment.

As a first estimate, we first measure the PSF-subtracted flux by placing a 04 radius aperture on each PSF-subtracted image (we do not apply aperture corrections to these measurements, but based on the measured effective radii of typical z ∼ 6 galaxies from JWST; Morishita et al. 2024; this aperture should capture ≳80% of each galaxy's flux). The host galaxy of J0731+4459 (z = 5.01, top row of Figure 2) is undetected in the noisy F410M image (measured flux at the <3σ level) and marginally detected (∼4σ) in F430M. J1340+2813 (z = 5.36, second row) is undetected in F210M and marginally detected (4–5σ) in F410M and F430M. J2239+0207 (z = 6.25, third row) is undetected in F210M but returns fluxes at the 6.8 and 5.5σ levels in F360M and F480M, respectively (as discussed in S23), our only apparently unambiguous detection with this PSF normalization method. J1148+5251 and J1120+0641 (bottom two rows, the two highest-redshift quasars in our sample) are undetected in all bands. That is, the only images in Figure 2 where the host galaxy emission can be unambiguously seen are in the F360M and F480M bands for J2239+0207. These results are consistent with the radial profiles (Figure 1), where only J2239+0207 shows obvious deviations between the radial profiles of the quasar and its PSF star.

We report the flux in all images within a 04 radius aperture in column (5) of Table 1. This method of normalizing the reference PSF necessarily removes any flux from the host galaxy near the center of the image, and may lead to oversubtraction, potentially making the values underestimates, not true upper limits. Also, there are additional systematic errors associated with the PSF subtraction. We address these issues in the next section (Section 3.2), where we determine upper limits to the possible host galaxy fluxes to accompany the potentially oversubtracted fluxes derived here.

Given how faint the host galaxies are, the dominant errors in measuring them (or obtaining limits) are the systematic ones associated with the PSF subtraction, not statistical errors. In addition, because the host galaxies of our quasars are very faint, the usual practice of subtracting the PSF down to the projected galaxy emission is not feasible. Instead, as described in the preceding section, we scaled the stellar PSF to match the quasar PSF at the center, which necessarily reduces the residual flux to zero near the center of the image. Placing an aperture on the PSF-subtracted image therefore returns essentially a lower limit on the total galaxy flux, because any flux at the very center of the galaxy (behind the PSF artifacts) has been subtracted. To estimate the amount of "missing" flux and provide an alternate, more conservative measurement of the galaxy flux and therefore its mass, we iterate the normalization of our reference PSF and fit Sérsic profiles to the residual flux in the PSF-subtracted images.

We begin by modifying the normalization of the reference PSF. When we reduce the normalization of the reference PSF to the quasar PSF, the resulting profile is broadly consistent with a Sérsic profile for all of our quasars, especially at r > 02 (see Figure 3). Early JWST results (e.g., Yang et al. 2022; Ormerod et al. 2024) have found that Sérsic indices of n ∼ 1 are typical for galaxies at z ∼ 6; indeed, we find that this is the Sérsic index most consistent with the profiles of our PSF-subtracted images. While a higher Sérsic index may provide a slightly better fit to the very inner, extremely noisy part of the image (at r < 01), any index above n ∼ 1.5 does not fit the region from 015 to 04 well. The quasar PSF, on the other hand, is highly inconsistent with a Sérsic profile. Therefore, by decreasing the flux of the reference PSF before subtracting and examining the resulting modified PSF-subtracted image and radial profiles until the PSF pattern becomes visible in the "subtracted" image and the radial profile is no longer well fit by an n = 1 Sérsic profile, we can place an upper limit on the flux attributable to the galaxy rather than to the PSF. We adjust the normalization to a range of values, such that the measured flux in the 04 aperture returns from 2 to 20 times the flux obtained when normalizing the stellar PSF to the quasar PSF at the center of the image. An example of this process is shown in Figure 3, for the F360M image of J1340+2813.

Zoom In Zoom Out Reset image size

We examine the radial profiles of these adjusted PSF-subtracted images, and fit the region of the profile relatively unaffected by central PSF artifacts (at radii ≳ 015), with an n = 1 Sérsic profile. We set the effective radius of the model Sérsic using Equations (5) and (6) of Morishita et al. (2024), determining an "expected" mass from the local M_BH–M_* relation and the black hole masses of our quasar sample; these effective radii are between 1 and 1.5 kpc for our quasars.

We find only small adjustments to the normalization are needed to account for oversubtraction. For most bands, when tuned to return more than 4–6 times the original, normalized-to-zero-at-center flux, the PSF pattern begins to become apparent in the subtracted images, and the core of its radial profile becomes inconsistent with a Sérsic profile (see Figure 3). The exception is the F410M image of J0731+4459, which displays worse stellar-quasar PSF agreement than the rest of the sample and therefore a noisier radial profile; we can increase its flux by approximately a factor of 20 before the radial profile becomes inconsistent with the Sérsic. We report the upper limit on the J0731+4459 F410M flux in Table 1, but do not use this band to constrain the galaxy mass. We take the maximum normalization that returns (a) a subtracted image without obvious residual PSF patterns, and (b) a radial profile that can be fit by an n = 1 Sérsic profile, as the upper limit on the flux of the galaxy as listed in column (6) of Table 1. The resulting PSF-subtracted images are shown in Figure 4.

Our five quasars possess high-mass SMBHs ($\mathrm{log}({M}_{\mathrm{BH}}\,[{M}_{\odot }])=8.8$ to 9.8); the local M_BH–M_* relation therefore predicts large host galaxy masses, between $\mathrm{log}({M}_{* }\,[{M}_{\odot }])=11$ and 12. Such massive hosts would have been easily detected. However, obtaining realistic host galaxy masses is nontrivial, as each quasar has been observed in no more than three closely spaced bands. It is therefore not feasible to determine the host masses by performing, e.g., spectral energy distribution (SED) fitting because of the degeneracies in models due to our small number of data points, combined with the low signal-to-noise ratio. We must therefore turn to another method.

The Spitzer Space Telescope observed hundreds of high-redshift AGN and inactive galaxies during its mission, including many galaxies in legacy fields with additional data spanning the ultraviolet to far-infrared or radio, and therefore with robust masses from SED fitting. We can therefore use these galaxies to examine the relationship between their IRAC Band 1 (3.6 μm) fluxes, whose bandpass closely resembles NIRCam F360M, and their stellar masses from SED fitting. At the redshifts of our quasars, F360M probes the rest-frame optical and should therefore scale with the galaxy mass—we can use this to build an empirical relationship between 3.6 μm flux and mass to constrain the masses of our quasar host galaxies. This method is based on the assumption underlying the Magorrian relation that quasar host galaxies, as a population, should resemble field galaxies at the same redshift.

McLure et al. (2011), Curtis-Lake et al. (2013), and Jiang et al. (2016) model massive z ∼ 5–8 galaxies with high-quality observations spanning a wide range of wavelengths, including IRAC Band 1 measurements. ² These samples include high-mass, highly star-forming galaxies similar to those we expect to host our luminous quasars. Because these sample galaxies—as well as our quasar hosts—lie at a range of redshifts, the IRAC Band 1 filter probes a range of rest-frame wavelengths. This may introduce biases if the galaxies' spectra are not flat, requiring complicated K-corrections. Reassuringly, however, when we create model SEDs with reasonable star formation histories and the small levels of reddening typical of galaxies at these redshifts (E(B − V) ∼ 0.05), we find that they are virtually flat in frequency units at the wavelengths of interest (rest-frame ∼ 450–600 nm). The necessary K-corrections are therefore negligible.

Zoom In Zoom Out Reset image size

We correct the IRAC Band 1 fluxes of the McLure et al. (2011), Curtis-Lake et al. (2013), and Jiang et al. (2016) data to a single redshift (z = 6), and plot their masses as a function of IRAC Band 1 flux in Figure 5. The errors on the McLure et al. (2011) and Curtis-Lake et al. (2013) galaxies are directly from the references, but the quoted errors in Jiang et al. (2016) are much smaller and appear not to include systematic uncertainties. We therefore increase the error bars on the Jiang et al. (2016) points by 0.15 dex, making them similar to the errors quoted in the other references to ensure those points do not dominate our fits.

Zoom In Zoom Out Reset image size

The points in Figure 5 clearly fall into two distinct populations. Very young galaxies, with a dominant stellar population age (from SED fitting) less than 10 Myr have a larger 3.6 μm flux for a given mass; these are plotted as open circles in Figure 5. Older galaxies, on the other hand, which we expect to resemble the galaxies hosting our quasar sample more closely, produce less 3.6 μm flux at a given mass (filled circles in Figure 5). We want to place maximally conservative upper limits on the masses of our host galaxies; including the very young galaxies biases our fit, tending to return a lower galaxy mass for a given flux. Using such a fit may lead us to conclude mistakenly that our quasar host galaxies lie above the local M_BH–M_* relation even if they do not. A more conservative result is obtained by rejecting these very young galaxies when we fit the relationship between stellar mass and 3.6 μm flux.

We therefore fit the relationship between the 3.6 μm flux density (adjusted to z = 6) and galaxy mass measurements by least squares minimization, holding the slope to 1 (i.e., the mass should be proportional to the 3.6 μm flux density) and excluding galaxies with a dominant stellar population age < 10 Myr. The reduced χ² for this fit is 1.98; if the single most discordant point (the value at $\mathrm{log}({F}_{\nu }\ [\mu \mathrm{Jy}])=-0.31$, $\mathrm{log}({M}_{* }\,[{M}_{\odot }])=9.17$) is rejected, the reduced χ² is 1.50. The equation of this best-fit line is

$$\begin{eqnarray}&&\mathrm{log}({M}_{* }\,[{M}_{\odot }])=\mathrm{log}({F}_{\nu }\ [\mu \mathrm{Jy}])+10.20,\end{eqnarray} \tag{ 1 }$$

where F_ν is the 3.6 μm flux adjusted to z = 6 in units of μJy and the galaxy mass is in solar units.

With this empirical relation, we can place lower and upper limits on the masses of our quasar host galaxies from the lower and upper limits on their fluxes. As with the galaxies in the models, we begin by correcting our lower limits (from aperture photometry, column (5) of Table 1) and upper limits (from varying the normalization in Section 3.2, column (6) of Table 1) on the galaxy fluxes to z = 6; we again make the assumption that the galaxy SED is flat at the wavelengths of interest, around rest-frame 500 nm. We can then use these K-corrected fluxes and Equation (1) to obtain a mass range for each quasar host galaxy. For the galaxies with an F360M flux (J2239+0207, J1148+5251, and J1120+0641), we utilize that flux to calculate the mass. For J1340+2813, we use the band at the most similar wavelength, F410M. J0731+4459 also has an F410M flux, but because its F410M PSF-subtracted image displays more severe PSF-subtraction artifacts than its F430M image, suggesting some degree of mismatch between its quasar and stellar PSF in that band, we conservatively use F430M to calculate its mass instead. Again, the rest-frame optical SED is very flat for these galaxies and using F410M or F430M instead of F360M should not probe a significantly different part of the SED. The masses are reported in column (7) of Table 1.

We obtained NIRSpec spectra of three of our quasars: J0731+4459 (z = 5.01), J1340+2813 (z = 5.36), and J2239+0207 (z ∼ 6.25) to constrain the strength of the Balmer break and measure black hole masses from Hα and Hβ. All three quasar spectra are shown in Figure 6 in order of increasing redshift. The prism spectrum of J2239+0207 (bottom row of Figure 6) is extensively discussed in S23: the spectral resolution and signal-to-noise ratio are low, but we are nonetheless able to use Hα to place a lower limit on the black hole mass of $\mathrm{log}({M}_{\mathrm{BH}}\,[{M}_{\odot }])\gt 8.5$, consistent with previous measurements from Mg ii and C iv. The redshifts of both J0731+4459 and J1340+2813 place Hβ within the NIRSpec wavelength range, which we use to calculate an additional black hole mass estimate.

Zoom In Zoom Out Reset image size

We do not observe an obvious narrow component of Hβ in the spectra of J0731+4459 and J1340+2813. However, both spectra display a [Fe ii] emission feature blended with Hβ, which must be fit simultaneously with the Hβ line. We therefore fit two Gaussian profiles—one for the broad component of Hβ and one for the [Fe ii] feature—superimposed on a local linear continuum. We fix the central wavelengths to the line centers and allow the amplitudes and widths of both lines to vary. For J0731+4459, we obtain a FWHM ∼ 3400 km s⁻¹. This line width can be converted to a black hole mass using Equation (5) in Vestergaard & Peterson (2006)

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{log}{M}_{\mathrm{BH}}=\mathrm{log}\left\{{\left[\displaystyle \frac{\mathrm{FWHM}({\rm{H}}\beta )}{1000\ \mathrm{km}\,{{\rm{s}}}^{-1}}\right]}^{2}{\left[\displaystyle \frac{\lambda {L}_{5100{\rm{\mathring{\rm A} }}}}{{10}^{44}\ \mathrm{erg}\,{{\rm{s}}}^{-1}}\right]}^{0.5}\right\}\\ \,\,\,+(6.91\pm 0.02)\end{array}\end{eqnarray} \tag{ 2 }$$

where we use the flux density in F410M (Table 1) as a proxy to determine the continuum luminosity at 5100 Å. We find a black hole mass $\mathrm{log}({M}_{\mathrm{BH}}\,[{M}_{\odot }])={9.3}_{-0.3}^{+0.2}$.

For J1340+2813, we measure a FWHM of 5100 km s⁻¹, and using Equation (2) obtain mass $\mathrm{log}({M}_{\mathrm{BH}}\,[{M}_{\odot }])={9.4}_{-0.4}^{+0.2}$. This value is consistent with that reported by Jun et al. (2015).

J0731+4459 and J1340+2813 also have rest-frame 3650 Å positioned within the NIRSpec medium wavelength range. While their host galaxies are only marginally detected in the NIRCam images, their spectra might reveal a Balmer break if there is a massive host galaxy sufficiently compact to hide behind the PSF-subtraction artifacts.

To test for this possibility, we fit a Balmer break spectrum to the quasar spectra, varying the strength of the break. We created a model of the Balmer break using PopStar models (Mollá et al. 2009) assuming a Chabrier initial mass function. The strength of any Balmer break depends on the age of the stellar population dominating the SED; a very young population will have a very weak feature. For a baseline model, we have assumed constant star formation for 500 Myr. To fit the resulting complex spectrum to the quasar spectrum, we first mask all absorption and emission features in both spectra and fit the model PopStar spectrum with two lines, shortward and longward of the Balmer break, and a cubic function connecting the two to represent the break itself. We fit this parameterization of the Balmer break to the spectra of our quasars, with the normalization of the model, the magnitude of the break, and an overall slope (to represent any effects of dust extinction) as free parameters. We allowed negative break values as a test for fitting residuals.

The best-fit strengths of the Balmer break are approximately 0.8% and 1.4% the flux of the quasar continuum in J0731+4459 and J1340+2813, respectively, as shown in Figure 6, and both are negative (i.e., the fit returns a reverse break, with the spectrum brighter at blue wavelengths than red). We do not believe this is significant, as in both cases, no break at all is consistent with the spectra and fitting a single flat continuum returns a nearly identical reduced χ².

Placing an upper limit on the Balmer break is therefore difficult, because the best-fit value of the break is negative and consistent with zero; clearly the contribution of the stellar continuum to the spectrum is very small. As a rough estimate of the maximum possible contribution, we take upper limits to the nominal break values as the negatives of the best fits (so the breaks are positive). These stellar continuum fluxes permit galaxy flux densities for J0731+4459 and J1340+2813 of 2.5 and 5.4 μJy, respectively, in the 3–4 μm bands, values consistent with the upper limits in Table 1 for the same quasars. This result, and the lack of obvious stellar contribution when visually inspecting the quasar spectrum, implies that there is no significant flux (relative to the quasar flux) from the galaxy in the central regions of the PSF, and rules out the possibility of compact host galaxies sufficiently massive to preserve the local M_BH–M_* relation hiding behind PSF-subtraction residuals in the core of the PSF.

A theoretical example of the growth of a galaxy and SMBH together is provided by Li et al. (2008; see their Figure 1). In these particular models, the central SMBH and its host galaxy have reached their final masses by z = 5, with nearly all growth (particularly of the host galaxy) complete even earlier, by z = 6. As such, between z = 6 and z = 5, the host galaxy is slightly overmassive compared to the black hole.

Habouzit et al. (2022) discuss the results of six large cosmological simulations with different approaches to SMBH–galaxy coevolution, providing a more complete view. Such studies are limited by the currently achievable simulation volumes (see Habouzit et al. 2022 for a discussion), but the two with adequate statistics to extend to high-luminosity quasars indicate that at z = 5 they reside in host galaxies in the 0.3 to 1 × 10¹¹ M_⊙ range, with relatively little scatter. ³ However, the quasars we have studied have SMBH masses well above the simulation results. ⁴ J2239+0207 lies at the very high end of the simulated distribution (S23), but our other four quasars have much larger SMBH masses: for J0731+4459, 2 × 10⁹ M_⊙ (this paper); J1120+0641, 1.35 × 10⁹ M_⊙ (Yang et al. 2021); J1148+5251, 7 × 10⁹ M_⊙, (Barth et al. 2003; Willott et al. 2003; Shen et al. 2019); ⁵ and J1340+2813, 2.5 × 10⁹ M_⊙, (Jun et al. 2015; this paper). The additional quasars analyzed for host galaxies by Yue et al. (2023) also all have very massive SMBHs, as listed in that paper, ranging from 1.9 × 10⁹ to 1.2 × 10¹⁰ M_⊙. The AGN luminosities are proportional to their SMBH masses (Shen et al. 2019), placing all of these quasars high into the "bright" category (Habouzit et al. 2022). According to the predictions from existing simulations, we therefore expect these SMBHs to lie in hosts of 1 × 10¹¹ M_⊙ (and larger), significantly more massive than the upper limits we have calculated for our quasar sample. We now consider whether the apparent contradiction to these predictions we have found might arise from other effects than the intrinsic quasar/host behavior.

4.1.1. Effect of Intense Nuclear Star Formation

4.1.2. Host Galaxy Extinction

Even with our conservative upper mass limits, all five of our quasars lie significantly above the local M_BH–M_* relation: the central SMBHs are overmassive compared to their host galaxies. Similar behavior is found for the EIGER quasar sample, whose members also lie above the local M_BH–M_* relation (Yue et al. 2023; pink circles in Figure 7, though they have significantly larger radii for their host galaxies than we find for ours). In contrast, the two quasars whose host galaxy emission was detected via subtraction of an empirical PSF by Ding et al. (2023; purple diamonds in Figure 7) have host galaxies that are, within their uncertainties, roughly as massive as predicted from their black hole masses by the local M_BH–M_* relation of Kormendy & Ho (2013).

Zoom In Zoom Out Reset image size

Overall, these results strongly suggest that the SMBHs in the high-redshift quasars are overmassive compared with their host galaxies. However, it is critically important to consider the selection effects at play when drawing conclusions from high-redshift quasar samples. In particular, our sample of luminous z ∼ 6 quasars and hosts is subject to two significant biases.

The first is simple: more massive host galaxies are significantly easier to detect and constrain after performing PSF subtraction. If z ∼ 6 quasars of a given black hole mass display a range of host galaxy masses, the more massive hosts (i.e., those galaxies lying to the right of the M_BH–M_* relation) are more likely to be detected. This bias, however, is contrary to our results. No z ∼ 6 quasar host galaxies detected with JWST to date lie to the right of the local relation (as can be seen in Figure 7).

The situation is complicated, however, by the selection bias first described by Lauer et al. (2007). Due to the steep slope at the bright end of the galaxy luminosity function, the most massive black holes are more likely to lie in lower-mass galaxies (where the black holes are overmassive compared to the local M_BH–M_* relation), than in high-mass galaxies where the massive black hole is proportional. Because massive black holes can accrete at higher rates and produce brighter AGN, this bias may then also manifest in flux-limited quasar samples. Therefore, to compare our sample to local AGN and draw conclusions about the M_BH–M_* relation at high redshift accurately, we must use a similarly biased sample of the most luminous quasars in the local Universe.

As discussed in more detail in S23, we create this comparison z ∼ 0 sample from the Palomar-Green (PG) sample of Type 1 quasars at z ≤ 0.5. This provides a valid comparison because the high-redshift quasars are remarkably similar in all the relevant properties to very luminous local Type 1 quasars (Fan et al. 2023). By examining PG quasars with L_AGN ≥ 2 × 10¹² L_⊙ (shown in Figure 7 as open triangles), we recreate the same luminosity bias present in our high-z sample and others in the literature (e.g., Ding et al. 2023; Yue et al. 2023). The stellar and black hole masses we retrieve for the PG quasars are derived using methods as similar as possible to those used for our z ∼ 6 quasar sample: stellar masses from rest-frame optical photometry and black hole masses from single-epoch spectra (see S23 for more details). If these z = 0 luminous quasars behaved like the z ∼ 6 JWST samples, there would be a significant population lying to the left of the local relation: this is not the case. The local quasars lie on or to the right of the local relation, with only a very small minority of the sample (four quasars out of 37 total) slightly to the left. Indeed, a two-sample Kolmogorov–Smirnov test comparing the distance above/below the local M_BH–M_* relation for the high-L_AGN PG quasars and the z ∼ 6 samples of this work, Ding et al. (2023), and Yue et al. (2023) returns a p-value of 4 × 10⁻⁶, which indicates these samples are drawn from distinct populations (despite both being subject to the same luminosity bias). The Lauer et al. (2007) bias therefore cannot be entirely responsible for the increasing number of high-redshift quasars far from the local M_BH–M_* relation.

It is therefore unlikely that the position of our quasars above the local M_BH–M_* relation is only a consequence of extinction effects on our estimates of the host masses, or of selection biases. Instead, we find that even the most massive SMBHs at z ∼ 6 share the behavior found for lower-mass ones (e.g., Goulding et al. 2023; Harikane et al. 2023; Kokorev et al. 2023; Maiolino et al. 2023), namely to show a large scatter in M_BH/M_* with a strong tendency toward overmassive SMBHs compared with the local Magorrian relation. The small number of SMBHs discovered so far at even higher redshifts (Übler et al. 2023; Bogdan et al. 2024; Goulding et al. 2023; Kokorev et al. 2023; Pacucci et al. 2023) share this behavior.

The spread of M_BH–M_* ratios observed—from galaxies consistent with the local relation to cases an order of magnitude less massive than predicted—aligns with existing measurements of [C ii] dynamical masses using ALMA and other submillimeter telescopes (see the small gray diamonds in Figure 7). Indeed, these ALMA results eliminate the possibility that the local M_BH–M_* relation holds for gas masses rather than stellar masses at z ∼ 6; ALMA dynamical masses at this redshift also tend to indicate galaxies less massive than predicted by the local M_BH–M_* relation. However, because the gas dynamics can be influenced by nongravitational forces and the mass estimates are dependent on assumptions about the gas distribution around the SMBH, these results should be interpreted with caution until stellar masses can be determined (Inayoshi et al. 2020). An illustrative case is the host galaxy of J2239+0207, with a stellar mass more than an order of magnitude less than its dynamical mass (S23). Such galaxies are not unheard of in the high-redshift Universe (see, e.g., Riechers et al. 2013, for an example of a z = 6.34 nonactive galaxy with a stellar mass nearly an order of magnitude less than its dynamical mass).

Clearly, not all galaxies—or even all AGN—at z ∼ 6 will lie significantly above the local relation; even some luminous quasars with detected host galaxies are consistent with it. However, the number of luminous quasars observed to have overmassive black holes relative to the stellar masses of their host galaxies (a rarity at z = 0), taken together with similar results from submillimeter telescopes, implies that there is additional scatter in the M_BH–M_* relation at high redshift; the relationship between galaxy mass and black hole mass is not nearly as tight as in the local Universe. 1 Gyr after the Big Bang, SMBHs and their host galaxies do not yet march in lockstep.

These results appear to contradict predictions that the growth of galaxy mass will precede that of SMBH mass such as in Li et al. (2008). They are also problematic for the some of the simulations summarized by Habouzit et al. (2022), which predict reasonably coordinated coevolution. These results suggest that the growth of the SMBHs is not purely due to galaxy mergers, where the overall host galaxy mass grows and much of the interstellar gas falls to the center where it can be accreted by the SMBH. Nonetheless, the presence of metals and dust around these quasars (e.g., Maiolino et al. 2003; Venemans et al. 2018) shows that they must lie in evolved populations of stars.

The mystery of the rapid SMBH growth therefore persists (Fan et al. 2023). Solutions may well require ways that the growth of SMBHs could be turbocharged to proceed much more rapidly than in the conventional major merger scenario, as discussed for example by Inayoshi et al. (2020).