# Detection of the Low-Stellar Mass Host Galaxy of a $z\sim 6.25$ Quasar with
JWST

Meredith A. Stone Steward Observatory, University of Arizona, 933 North Cherry
Avenue, Tucson, AZ 85721, USA  Jianwei Lyu (吕建伟) Steward Observatory,
University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA  George
H. Rieke Steward Observatory, University of Arizona, 933 North Cherry Avenue,
Tucson, AZ 85721, USA  Stacey Alberts Steward Observatory, University of
Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA

(Received April 10, 2023; Revised July 5, 2023; Accepted July 28, 2023)

###### Abstract

We characterize the stellar mass of J2239+0207, a $z\sim 6.25$ sub-Eddington
quasar (M${}_{1450}=-24.6$), using dedicated JWST/NIRCam medium-band
observations of a nearby PSF star to remove the central point source and
reveal the underlying galaxy emission. We detect the host galaxy in two bands
longward of the Balmer break, obtaining a stellar mass of $\sim 10^{10}$ M⊙,
more than an order of magnitude less than this quasar’s existing measured [C
II] dynamical mass. We additionally calculate the mass of J2239+0207’s central
supermassive black hole using JWST/NIRSpec IFU observations, and determine
that the black hole is $\sim 15$ times more massive than predicted by the
local $M_{\mathrm{BH}}$-$M_{*}$ relation, similar to many high-redshift
quasars with dynamical masses determined via millimeter-wave line widths. We
carefully consider potential selection effects at play, and find that even
when $z\sim 6$ quasars are compared to a local sample with similarly
determined dynamical masses, many of the high-redshift quasars appear to
possess overmassive black holes. We conclude $z\sim 6$ quasars are likely to
have a larger spread about the $M_{\mathrm{BH}}$-$M_{*}$ relation than
observed in the local Universe.

††facilities: JWST(NIRCam, NIRSpec)††software: Astropy (Astropy Collaboration
et al., 2013, 2018), Matplotlib (Hunter, 2007), NumPy (van der Walt et al.,
2011), photutils (Bradley et al., 2022), WebbPSF (Perrin et al., 2014)

##  1 Introduction

Active galactic nuclei (AGN) are energetic enough to greatly impact the star
formation, morphologies, and other properties of their host galaxies, but the
nature of these effects—especially at high redshift—is still poorly
understood. Understanding the interplay between galaxies and their central
supermassive black holes (SMBHs), including their coeval growth and triggering
of the AGN phase, is critical to building a complete picture of galaxy
formation and evolution.

Studying powerful AGN—particularly quasars—and their host galaxies is
difficult. Black hole masses are easiest to measure up to high redshift in
Type 1, face-on AGN via single-epoch spectra, where the Doppler-broadened
emission lines from gas orbiting the black hole are visible in the spectrum.
However, in these systems the bright quasar tends to outshine the emission
from its host, making the characterization of the host galaxy challenging.

At high redshift, the properties of quasar host galaxies can be explored via
the rest-frame far-infrared and submillimeter regimes, where the contribution
from the AGN is low and the emission is dominated by cold gas and dust tracing
regions of star formation (e.g., Lyu et al., 2016). The width of the [C II]
$158\mu$m fine-structure cooling line in particular, accessible at high
redshift with ALMA and other ground-based submillimeter telescopes, can be
used to probe the dynamical mass of quasar hosts (e.g. Wang et al., 2013;
Venemans et al., 2016; Decarli et al., 2018; Izumi et al., 2019) in systems
where the stellar emission is difficult to distinguish from the quasar (Shao
et al., 2017; Pensabene et al., 2020; Neeleman et al., 2021; Shao et al.,
2022). However, this method is subject to a number of assumptions and
uncertainties regarding the nature of the galaxy and its orientation on the
sky. To improve on this approach, submillimeter observations can be used to
create spatially-resolved line maps and model the dynamics of the host
galaxies. Dynamical masses determined in these ways are often compared with
local relations based on more directly determined stellar masses when
comparing high-redshift quasar samples to the local Universe via e.g. the
$M_{\mathrm{BH}}/M_{*}$ relation (Izumi et al., 2019; Pensabene et al., 2020;
Yue et al., 2021).

Another, more direct method to detect the underlying stellar emission from
quasar host galaxies and estimate their stellar masses is to subtract a
theoretical or empirical point-spread function (PSF) from quasar images.
However, these efforts have mainly been limited to $z\lesssim 2-3$ (e.g.
Gehren et al., 1984; Hutchings et al., 2002; Dunlop et al., 2003; Jahnke et
al., 2004; Kim et al., 2017; Yue et al., 2018; Li et al., 2021), using ground-
based optical telescopes or the Hubble Space Telescope to probe the rest-frame
optical to near-infrared emission of the quasar where the host galaxy
contribution is the strongest. These studies have revealed that SMBH masses
($M_{\mathrm{BH}}$) and the masses of their host galaxies ($M_{*}$, or the
central bulge mass $M_{\mathrm{bulge}}$) are correlated in the low- and
intermediate-redshift Universe (Kormendy & Ho, 2013; Reines & Volonteri,
2015). The question follows: does this correlation persist to high redshift? A
steeper relation at high-z—i.e. overmassive black holes relative to their host
galaxies—implies that black holes build up their masses before their host
galaxies do, whereas greater scatter at high-z implies that the coeval growth
of supermassive black holes and galaxies has not yet settled to its local
behavior.

Most of the small number of distant ($z\gtrsim 3$) quasars whose host galaxy
emission has been detected or constrained (e.g. Peng et al., 2006; McLeod &
Bechtold, 2009; Mechtley et al., 2012; Targett et al., 2012; Schramm et al.,
2019; Marshall et al., 2020, 2023) lie above the established
$M_{\mathrm{BH}}$-$M_{*}$ relation, with overmassive black holes given their
stellar mass; on the other hand, studies of high-$z$ quasars with dynamical
masses measured from [C II] have implied $M_{\mathrm{BH}}$-$M_{*}$ relations
both above and consistent with the local Universe (Walter et al., 2004; Wang
et al., 2013; Venemans et al., 2016; Decarli et al., 2018; Izumi et al., 2019;
Farina et al., 2022). However, this observed trend may be due in part to
selection effects: at such high redshifts, studies are biased towards common,
lower-mass galaxies with overmassive black holes (Lauer et al., 2007; Zhang et
al., 2023). Additionally, at these redshifts, optical telescopes probe rest-
frame ultraviolet emission of galaxies, rather than the rest-frame optical
that directly traces the stellar mass. It is therefore important to study
underluminous quasars accreting below the Eddington limit—and in particular
their rest-frame optical emission— to mitigate this bias and explore the
potential scatter of the high-redshift $M_{\mathrm{BH}}$-$M_{*}$ relation.

JWST’s NIRCam instrument (Rieke et al., 2005, 2023) will greatly expand the
number of high-redshift quasars with host galaxy detections by probing stellar
emission out to $z\sim 9$ at high sensitivity and resolution. High-$z$ quasar
samples being built with NIRCam will help to fill out the high-$z$
$M_{\mathrm{BH}}$-$M_{*}$ relation to better understand the details of galaxy
formation and evolution in the first $\sim 1$ Gyr after the Big Bang (see e.g.
Ding et al., 2022, for other early results). To demonstrate this possibility,
we designed a JWST GTO program (ID 1205; Rieke et al. 2017) to obtain NIRCam
imaging and/or NIRSpec IFU observations of six $z\sim 5-6.5$ quasars spanning
a range of host galaxy star formation rate, AGN luminosity, and obscuration.

In this paper, we report our first results for the $z\sim 6.25$ quasar HSC
J2239+0207 with multi-band NIRCam imaging and NIRSpec prism IFU observations.
Discovered by the Subaru High-$z$ Exploration of Low-Luminosity Quasars
(SHELLQ) project (Izumi et al., 2019), J2239+0207 is a moderately-luminous
quasar ($M_{\rm 1450,optical}=-24.7$ or $L_{\rm bol}\sim 2.5\times
10^{12}~{}L_{\odot}$) that is accreting at $\sim$30% of the Eddington limit,
though its precise Eddington ratio is uncertain due to the spread in measured
BH masses: $M_{\mathrm{BH}}=6\times 10^{8}$ and $1\times 10^{9}\;M_{\odot}$
from C IV $\lambda 1549$ and Mg II $\lambda 2798$ respectively (Onoue et al.,
2019). The quasar’s redshift situates the stellar $4000\,\mbox{\AA}$ break at
approximately 3 $\mu$m, amidst the long-wavelength NIRCam filters. Due to its
relatively low Eddington ratio and high dynamical mass ($\log\,$Mdyn
(M${}_{\odot})$$=11.5$, Izumi et al., 2019), the host galaxy of J2239+0207
should be easily detectable in the rest-frame optical. J2239+0207 is therefore
an ideal case for host galaxy detection and characterization of a quasar at
the end of the Epoch of Reionization with JWST.

This paper is organized as follows. Section 2 summarizes the JWST observations
of J2239+0207 and the relevant data reduction processes. We outline the
procedures to extract the host galaxy signal from the NIRCam images and to
determine the black hole mass from the NIRSpec spectrum, and present our
measurements in Section 3. After characterizing its host galaxy mass, we
determine the location of J2239+0207 relative to the local
$M_{\mathrm{BH}}$-$M_{*}$ relation and, considering various selection biases,
discuss its implications on the cosmic coevolution of supermassive black holes
and their host galaxies in Section 4. A final summary is given in Section 5.
Throughout this work, we assume $H_{0}$ = 69.6, $\Omega_{M}$ = 0.286, and
$\Omega_{\Lambda}$ = 0.714.

##  2 Observations and Data Reduction

###  2.1 NIRCam

HSC J2239+0207 was observed in three NIRCam medium-band filters, with total
exposure times of 2227.2 seconds in F210M and 1113.6 seconds in F360M and
F480M. We used Module B SUB400P, (FOV 12.5$\arcsec\times 12.5\arcsec$, SW;
$25\arcsec\times 25\arcsec$, LW), to image the quasar field using the BRIGHT2
readout mode. To improve the PSF sampling and mitigate cosmic rays and
detector artifacts, we adopted 4$\times$4 primary and sub-pixel dithering
patterns with dither types INTRAMODULEBOX and STANDARD, respectively. To
obtain a reference PSF, the nearby bright star 2MASS 2M22390990+0207329 was
also imaged in the same filters and observing setup with the RAPID readout
mode and shorter exposure times to avoid saturation (424 seconds in F210M, 265
seconds in F360M, and 159 seconds in F480M). We selected these medium-band
filters to minimize changes in the point spread function (PSF) due to the
different shapes of the spectral energy distributions of the quasar and star
across the photometric bandpass: differences in the color of the quasar and
star have negligible effect (see Section 3.1).

We processed this data using the JWST pipeline version 1.8.4, following
roughly the standard procedures as recommended by the STScI
JWebbinars111https://github.com/spacetelescope/jwebbinar_prep/tree/main/imaging_mode.
The pipeline parameter reference file is jwst_1015.pmap, as registered in the
JWST Calibration Reference Data System (CRDS). During Stage 2, we added a
custom step to measure and subtract the image striping caused by the 1/f noise
from the detector readout in each frame. The final mosaic images were produced
by aligning and stacking all the frames obtained in Stage 2 using the Stage 3
pipeline. To improve the accuracy of PSF subtraction, these images were re-
sampled with smaller pixel scales at 0.0147$\arcsec$/pixel in F210M, and
0.0300$\arcsec$/pixel in F360M and F480M with the drizzling algorithm in the
Resample step. Finally, we conducted a global background subtraction using the
Photutils Background2D function. The background-subtracted quasar images in
the three bands are shown in Figure 1: all three bands are dominated by the
central PSF.

Figure 1: NIRCam images of J2239+0207 in our three science bands: F210M
(left), F360M (center), and F480M (right). All three images have a field of
view of $\sim 1.8$″. The F360M image, with a different stretch and greater
field of view, is shown inset into the center panel to show the emission-line
companion (circled) visible in both our NIRSpec spectrum and NIRCam image at a
distance of $\sim 1\arcsec$ from the quasar.

###  2.2 NIRSpec

To aid the analysis and interpretation of the host galaxy signals, we also
obtained NIRSpec IFU data of J2239+0207 with the PRISM/CLEAR disperser-filter
combination. This IFU data has a $3\arcsec\times 3\arcsec$ field of view with
the size of each IFU element at $0.1\arcsec\times 0.1\arcsec$, and wavelength
coverage from 0.6 $\mu$m to 5.3 $\mu$m with nominal spectral resolution
$\sim$100 but increasing significantly toward longer wavelengths. To mitigate
cosmic ray and detector artifacts, we adopted the SPARE-CYCLING dither type in
LARGE size at Points 1, 2, 3, 4. With the NRSIRS2 readout, the total
integration time was 8870 seconds.

We reduced the NIRSpec IFU data using the same JWST pipeline version and CRDS
file as the NIRCam data. Initially we followed the standard
pipeline,222https://github.com/spacetelescope/jwebbinar_prep/tree/main/ifu_session
but found it has several issues: for example, a few IFU elements at the quasar
center had been mistakenly identified as outliers by the default pipeline and
automatically masked out. We therefore performed a visual inspection of the
data cube before and after the outlier detection to differentiate the true
outliers from the mistakenly flagged bright quasar pixels. We produced a final
data cube for science analysis by skipping the automatic outlier detection
step, and instead manually masking out the IFU elements we believe to be true
outliers.

In this work, we use the IFU data to measure the quasar black hole mass and
determine if any host galaxy signal detected via PSF subtraction of the NIRCam
data is strongly contaminated by gas emission lines. We therefore extract a
quasar spectrum using a circular aperture at a radius of 0.35$\arcsec$, as
shown in Figure 2. Further analysis of the IFU data is reserved for a future
paper.

Figure 2: Observed NIRSpec spectrum of the HSC J2239+0207 nucleus with some
prominent emission lines identified. As a comparison, we also plot an average
spectral template derived from the low-$z$ quasars in the DR7 edition of the
SDSS Quasar Catalog (Schneider et al., 2010).

##  3 Analysis and Results

###  3.1 Retrieving Host Galaxy Signals via PSF Subtraction

Figure 3: Comparison of the observed radial profiles of J2239+0207 and the PSF
reference star in our three bands. Top: The radial profiles of the PSF star
(red dashed line) and quasar (solid blue line): the star image has been scaled
to match the quasar flux in the central pixels. The pixel sizes are
0$\farcs$015 for F210M and 0$\farcs$030 for the other two bands. The quasar
PSF does not appear broader than the stellar PSF in any band; the risk of a
spurious galaxy detection due to PSF mismatch is therefore low. Bottom: The
normalized difference in flux between the star and quasar profiles. Slight
mismatches between the shape of the quasar and stellar PSFs near the core
produce the dip and spike observed in F360M and F480M (orange highlighted
region). The host galaxy is visible via the excess flux at a radius of $\sim
1-1.5$ kpc/$\sim 0.2\arcsec$ in both F360M and F480M, but not in F210M (yellow
highlighted region). Finally, at radii greater than $\sim 0.4\arcsec$, noise
dominates and any excess observed is likely not meaningful (green shaded
region).

As there are not a large number of bright stars in the quasar field suitable
for building an empirical PSF, we obtained a dedicated observation of a PSF
star. The NIRCam images of the star were taken immediately after observing the
quasar to minimize temporal variation of the PSF between the observations. We
take the stellar image as the PSF and scale and subtract it from the quasar
directly, rather than building an empirical PSF from stars in the image, or
fitting a model that relies on assuming certain properties of the underlying
host galaxy. Wolff et al. (2023) evaluated a number of state-of-the-art PSF
subtraction techniques in a search for the debris disk around $\epsilon$
Eridani, and find this technique to be optimal. The PSF star and quasar are
also placed at the same location on the detector, which minimizes PSF
differences between the star and quasar due to distortions across the
detector. Any differences in sub-pixel alignment between the star and quasar
PSFs are mitigated by the dithering process.

Because the quasar and the PSF star have different spectral shapes within the
bands (the relatively flat optical spectrum of the quasar compared to the
declining Rayleigh-Jeans tail of the star) it is possible that the shapes of
their PSFs will differ. This effect will tend to introduce a spurious galaxy
detection, as the quasar’s redder spectrum will produce a broader PSF than the
star. The choice of NIRCam medium bands for this analysis helps to mitigate
this effect, but we also investigate any possible differences by comparing
NIRCam PSFs generated with WebbPSF (Perrin et al., 2014) using as input a
Rayleigh-Jeans spectrum and our NIRSpec spectrum of J2239+0207. The resulting
WebbPSFs display functionally identical radial profiles in F360M and F480M at
all radii, and are very similar in F210M: the median normalized difference
($\frac{\mathrm{Quasar}-\mathrm{Star}}{\mathrm{Star}}$) between the quasar and
star profiles at all radii is $\sim 3\%$ in F210M and $\sim 0.5\%$ in F360M
and F480M. The difference in flux obtained within a 0.4″-radius aperture
between the simulated star and quasar is $\sim 0.3\%$ in F360M and F480M and
$\sim 1.5\%$ in F210M. We can therefore confirm that the differences in the
PSF of the quasar and star have a negligible effect on this analysis,
particularly away from the PSF core.

To subtract off the quasar point-spread function from the images, we first
align the background-subtracted quasar and PSF star images in each band. We
normalize the stellar PSF image in each band, scaling it to match the flux of
the star to the flux of the quasar in the central pixels, where the quasar
flux is dominated most heavily by the point source rather than any underlying
galaxy emission. We subtract the scaled star image from the quasar image and
examine the residuals, verifying that this simple method of normalizing the
PSF does not leave obvious over- or under-subtraction artifacts in the
resulting PSF-subtracted images (see Figure 4).

The radial profiles of the quasar and the star (the latter scaled to match the
quasar in the central pixels) in all three bands are shown in Figure 3. The
quasar and stellar profiles are very similar within $\sim 2$ kpc of the center
of the image where we expect to observe galaxy emission, with only slight
evidence of a broader quasar PSF relative to that of the star very near the
PSF core, consistent with the WebbPSF predictions. However, the quasar F360M
and F480M profiles display low-level excess emission ($\sim 10$% of the quasar
signal at the relevant radii) at a distance of $\sim 1$ kpc or 0.2″(bottom
panels); no such strong excess is seen in F210M. For comparison, the maximum
deviation at any radius between the simulated WebbPSF quasar and stellar
radial profiles was $\sim 7\%$ in F360M and F480M.

The PSF-subtracted quasar images are shown in Figure 4. F210M, which at the
quasar redshift probes shortward of the 4000Å break, shows bright PSF
subtraction artifacts at the center but exhibits no evidence of underlying
galaxy emission, consistent with the match of the quasar and star radial
profiles in the leftmost panel of Figure 3. As hinted by the radial profiles,
some faint signal is visible around the quasar in F360M and F480M, which both
lie longward of the Balmer break. We measure galaxy fluxes within a circular
0.4$\arcsec$-radius ($\sim 2$ kpc) aperture of 0.34 $\mu$Jy in F360M and 0.44
$\mu$Jy in F480M (6.8 and 5.5$\sigma$ respectively), which correspond to
respective host-to-quasar flux ratios of 4.6% and 3.9%. We can also put a
secure $3\sigma$ upper limit on the F210M flux: all fluxes are listed in Table
1.

These fluxes likely slightly underestimate the flux of the host galaxy,
because the galaxy flux near the PSF core is subtracted along with the PSF. In
Section 4.1, we use the outskirts of the galaxy where the PSF artifacts are
not severe to model the galaxy as a Seŕsic profile and determine the degree to
which we are underestimating the host galaxy flux.

Figure 4: The quasar image in F210M (left), F360M (center) and F480M (right) after subtraction of the stellar PSF, normalized by the total (Poissonian + pixel) error budget. The $0.4\arcsec$-radius aperture used to extract the host galaxy flux is shown in white. All three images are shown with the same scaling and stretch. F210M does not display significant emission beyond the central PSF artifacts, but extra extended emission is visible in F360M and F480M. Orange, yellow, and green shaded regions, as in Figure 3, mark the radii at which the emission is dominated by PSF artifacts, host galaxy emission, and noise, respectively. All images have a field of view of $\sim 1.8$″.  Table 1: Measured quasar and galaxy properties for J2239+0207  Band | Quasar | Galaxy | Approximate | Inferred
---|---|---|---|---
| Flux | Flux | Galaxy SizeaaEffective radius of best-fit Sérsic profile (Section 4.1) | Mass
| ($\mu$Jy) | ($\mu$Jy) | ″ |  $\log$(M⊙)
F210M | $6.81\pm 0.04$ | $<0.15$ | N/A | N/A
F360M | $7.43\pm 0.05$ | $0.34\pm 0.05$ |  $\sim 1$ kpc | $10.0^{+0.3}_{-0.5}$
F480M | $11.32\pm 0.09$ | $0.44\pm 0.09$ |  $\sim 1$ kpc | N/A

A companion to the quasar at a distance of $\sim 1$″is visible in F360M (see
Figure 1, central panel inset). A preliminary spectrum extracted from the
NIRSpec data shows this companion only in emission lines. Its redshift
relative to the quasar suggests that it is either an outflow, or a galaxy
passing close by. More details of this companion will be presented in a future
paper.

###  3.2 Black hole mass from the NIRSpec spectrum

To measure the emission lines in the quasar spectrum, we adopt a continuum
based on the classic Elvis general result (Elvis et al., 1994), slightly
adjusted in slope to match the slope of the quasar’s continuum. It provides a
good fit, with no indication of abnormalities, and is a useful baseline
against which to identify emission lines. Prominent lines are indicated in
Figure 2.

Figure 5: Results from modeling the H$\alpha$ \- [NII] line blend. The figure
shows the normalized $\chi^{2}$ as a function of the ratio of peak line
strengths for H$\alpha$ to [NII]6548Å. H$\alpha$ must be about seven times
stronger than the shorter wavelength [NII] line and the preferred H$\alpha$
FWHM is about 3050 km s-1.

The H$\alpha$ line in our spectrum (Figure 2) can provide an independent
estimate of the black hole mass (the spectral resolution at H$\beta$ is too
low to resolve the line well, and our efforts to model it were unsuccessful).
At the observed wavelength of H$\alpha$, the prism spectrum has a
resolution333https://jwst-docs.stsci.edu/jwst-near-infrared-
spectrograph/nirspec-instrumentation/nirspec-dispersers-and-filters of $\sim$
1200 km s-1. The FWHM of the blended H$\alpha$ line is $\sim$ 3000 km s-1, so
the spectral resolution is adequate to measure its width. To evaluate the
contribution of the [N II] lines to the width of H$\alpha$, we build a model
of the three lines (H$\alpha$ and both [N II] lines), with the width of the [N
II] lines set to the spectral resolution and the ratio of their fluxes to the
theoretical value of 3. By varying the width of the H$\alpha$ line and its
strength relative to the [N II] lines, we can compare the resulting line
profile to our observed spectrum and calculate a $\chi^{2}$ value. The results
are shown in Figure 5. The preferred width is $\sim 3050\pm 250$ km s-1. To
check this result, we also measure the half width at half maximum blueward of
the line center assuming $z=6.2497$ (Izumi et al., 2019). Because the
$\lambda$ 6543 line is at least six times weaker than H$\alpha$ (see Figure
5), we expect the half width at half maximum not to be significantly affected
by blending. The result is 1450 km s-1, or a full width of 2900 km s-1, in
agreement within the expected errors.

We adopt 3050 km s-1. Removing the broadening due to the spectral resolution
yields a final width of 2800 km s-1. There may be a narrow component to the
H$\alpha$ line—at our resolution there is no way to tell. The presence of a
narrow component would make our width estimate too small: this width is
therefore a lower limit. The single-epoch black hole relation is calibrated
with H$\beta$, and there are subtle differences in the typical measured FWHMs
of the two lines. We use the relation in Greene & Ho (2005) (their equation 3)
to convert to an equivalent H$\beta$ FWHM of $\gtrsim$ 3090 km s-1. We can
then use the relations in Vestergaard & Peterson (2006) (their equation 5) to
convert the line widths to an estimate of the black hole mass, as:

$\log{M_{\mathrm{BH}}}=\log\Biggl{\\{}\\!\biggl{[}\frac{\mathrm{FWHM(H}\beta)}{1000\;\mathrm{km\,s}^{-1}}\biggr{]}^{2}\biggl{[}\frac{\lambda
L_{5100\mbox{\scriptsize\AA}}}{10^{44}\;\mathrm{erg\,s}^{-1}}\biggr{]}^{0.5}\Biggr{\\}}\\\
+(6.91\pm 0.02)$ (1)

We take the continuum luminosity at 5100 Å from the flux density in F360M,
7.34 $\mu$Jy, and correct it for the emission line contribution using the
spectrum. The resulting mass is $\gtrsim(3.5\pm 0.7)\times 10^{8}$ M⊙.

Onoue et al. (2019) estimated the mass of the black hole in J2239+0207 using
both the Mg II and C IV lines. They measured the FWHM of the lines to be,
respectively, $4670^{+910}_{-700}$ and $4630^{+1040}_{-1260}$ km s-1,
corresponding to mass estimates of $11^{+3}_{-2}$ and $6.3^{+2.0}_{-2.5}$
$\times 10^{8}$ M⊙. The latter value includes a correction for outflows
indicated by the blueshift of the line (Coatman et al., 2017); without this
correction the value would be $8.9^{+2.8}_{-3.4}\times 10^{8}$ M⊙. They also
find a redshift of z = 6.2499. Taking the three mass estimates together (3.5,
6.3, and 11 $\times 10^{8}\;M_{\odot}$), for a net value of $\sim 6\times
10^{8}$ M⊙, the Eddington ratio is $\sim 0.3^{+0.2}_{-0.15}$ (for
L${}_{bol}=2.5\times 10^{12}$ L⊙)444Estimated from M1450 (Izumi et al., 2019).
As shown by the differences among the three values, the uncertainty is
dominated by systematic errors.

##  4 Discussion

###  4.1 Host Galaxy

A SMBH of the mass discussed above would be expected to lie within a massive
host galaxy. Assuming the local relation between black hole and galaxy bulge
mass, the J2239+0207 black hole of $6\times 10^{8}$ M⊙ predicts a galaxy of
$1.2\times 10^{11}$ M⊙ (Kormendy & Ho, 2013).

To compare with the values from our observations, we need to relate host
galaxy mass to flux. Our low SNR detections of the host galaxy and small
number of filters are inadequate to support detailed modeling of its
characteristics via e.g., spectral energy distribution fitting. However, we do
test to see if our measurements are consistent with a plausible SED. The host
galaxy has a strong far infrared excess (Izumi et al., 2019), which indicates
a SFR of $\sim$ 280 M⊙/yr assuming a Kroupa IMF (Kennicutt & Evans, 2012): the
galaxy is in a starburst stage. We run a model555Our simulations utilize the
PopStar models (Mollá et al., 2009) for a Chabrier IMF, which we combine to
provide SEDs for different SFHs. We have modified the code to include emission
lines, tied to the level of star formation. These models do not include
massive binary stars, which would extend the duration of strong ionizing
fluxes by a factor of $\sim 1.5-2$ (e.g. Stanway et al., 2016), but would not
change our results in a significant way. assuming the galaxy is halfway
through a 20 Myr long burst that will account for a 10% increase in its mass,
which was built up by similar bursts distributed uniformally over the past 500
Myr. The measurements are consistent with the predicted SED provided the
extinction is $A_{V}\geq 0.5$ mag. This model is far from unique; for example,
if the high rate of star formation has a duration of 2 Myr rather than 10 Myr,
then models with $A_{V}>$ 0.2 are consistent with the measurements.

As illustrated by these examples, model-fitting to our three measurements is
highly degenerate. Two detections of moderate signal to noise and one upper
limit are just too limited to constrain models further. We therefore use
observations of massive galaxies at similar redshift but without active nuclei
to provide insights to the probable host galaxy properties. These galaxies
should include the appropriate range of star formation histories and levels of
extinction. In addition, many of them are measured in Spitzer/IRAC Band 1,
which coincides closely with our F360M observations. The model discussed in
the preceding paragraph demonstrates that this band will give a reliable mass
estimate: for example, if the recent star formation in the model is eliminated
so there has been none for the last 100 Myr, the mass required to produce the
same flux at 3.6 $\mu$m is only a factor of 1.6 higher.

To evaluate the range of host galaxy properties and probable mass, we draw on
masses from modeling of galaxies with redshifts between 5 and 8 (Labbé et al.,
2010; McLure et al., 2011; Curtis-Lake et al., 2013). We adjust the IRAC Band
1 measurements from these references to z = 6.25 and K-correct them. The
latter corrections are based on a model without the current starburst but
allowing various levels of extinction, which we set according to the E(B-V)
determined from the modeling in the cited papers. We impose a slope of 1, i.e.
the mass should be proportional to the 3.6 $\mu$m luminosity. The data suggest
a non-unity slope, but $\chi^{2}$ is only increased by a factor of 1.6 by
forcing the slope to unity. The result is:

$\log(M_{*})=\log(F(\nu))+10.05$ (2)

where the flux is in $\mu$Jy and $M_{*}$ in M⊙. The scatter around this
relation is 0.23 dex. The sample used for this calibration also gives insight
to the expected level of extinction; 1 galaxy (3%) has an indicated $A_{V}$ =
1.1, three more (10%) have values consistent with $A_{V}=$ 1, and the rest
(87%) have indicated values of $A_{V}\leq$ 0.5.

Figure 6: left: The PSF-subtracted F360M (top) and F480M (bottom) galaxy
images, with the galaxy emission visible around the central PSF subtraction
artifacts. When a Sérsic profile with a measured flux of 3 $\mu$Jy—expected
for a stellar mass of $10^{11}$ M⊙—and effective radius of $\sim 1$ kpc is
subtracted from each PSF-subtracted image (center), the galaxy emission is
strongly oversubtracted. Conversely, subtracting a profile with a flux of 0.4
$\mu$Jy (F360M) or 0.5$\mu$Jy (F480M) effectively removes the extended galaxy
emission without evidence of oversubtraction (right). This indicates that the
stellar mass of the galaxy is likely closer to $10^{10}$ M⊙.

Equation 2 returns a galaxy mass of $\log{M_{*}\,(M_{\odot})}=9.6\pm 0.3$ from
the F360M flux measured for the J2239+0207 host galaxy: this uncertainty
includes both the uncertainty on the F360M flux from our sampling with a
random set of apertures and the scatter in our empirical IRAC Band-1-stellar-
mass relation. There are, however, a number of adjustments needed to this
value. As already discussed, for consistency with our F210M upper limit,
$A_{V}\geq 0.5$ is needed; we adopt $A_{V}$ = 1 at the rest wavelength of the
F360M band (0.5 $\mu$m). This would put the host galaxy in the top 13% of
extinction among our calibrator sample. It would make the mass estimate
$\log{M_{*}\,(M_{\odot})}=10$. However, our model indicates that about 15% of
the flux in F360M comes from the H$\beta$ and [OIII] emission lines, which
will be absent in the calibration galaxies since they have much lower levels
of star formation. In addition, as discussed above, the contribution of the
starburst in the host galaxy raises its flux in the F360M band by as much as a
factor of 1.6 over that from less active galaxies such as those used for the
calibration. Both of these effects would result in an overestimate of the host
galaxy mass; we reflect this possibility with asymmetric error bars and quote
the mass as $\log{M_{*}\,(M_{\odot})}=10.0^{+0.3}_{-0.5}$.

The stellar mass is therefore significantly lower than the [C II] dynamical
mass estimated for the same galaxy by Izumi et al. (2019),
$\log{M_{\mathrm{dyn}}\,(M_{\odot})}\sim 11.5$. Oversubtraction of the PSF
could suppress the measured flux and therefore the measured mass. However,
using Equation 2 and adopting $A_{V}=1$, as above, predicts a flux of 3
$\mu$Jy in F360M for $\log{M_{*}}=11$ (increasing for less extinction): we can
test this prediction by subtracting a simulated galaxy with this flux from our
PSF-subtracted images. We create mock circular Sérsic profiles (Sérsic index
$n=1$) with fluxes of 3 $\mu$Jy within a $0.4$″-radius aperture—identical to
that used to extract the galaxy flux—and a range of effective radii, and
subtract this mock galaxy profile from the PSF-subtracted quasar image. We
find that regardless of the input effective radius, a 3 $\mu$Jy total flux
invariably severely oversubtracts the galaxy emission (see Figure 6, center
panels). Conversely, a Seŕsic profile with a $\sim 1$ kpc effective radius and
total flux of $\sim 0.4\;\mu$Jy in F360M and $\sim 0.5\;\mu$Jy in F480M, when
subtracted from the galaxy images, effectively removes any evidence of the
host galaxy emission (Figure 6, right panels). This Seŕsic subtraction
predicts total galaxy fluxes $\sim 15$% higher than those measured in the
apertures in F360M and F480M, accounting for the fraction of the galaxy flux
directly under the core of the PSF.

Moreover, reducing the factor by which we scale the stellar PSF—therefore
leaving more leftover flux after PSF subtraction and pushing the mass
higher—cannot increase the mass enough to be consistent with the dynamical
mass. Above a deviation of $\sim 5\%$ from the best-fit PSF scale factor (at
which point the resulting mass is still consistent within the errors with our
original measurement) leftover PSF artifacts become obvious in the final PSF-
subtracted images. A stellar mass of $\sim 10^{10}$ rather than $10^{11}$
$M_{\odot}$ is therefore reasonable.

Like many other quasars studied at similar redshifts (e.g. Izumi et al., 2021;
Pensabene et al., 2020; Harikane et al., 2023), J2239+0207 lies above the
local $M_{\mathrm{BH}}-M_{\mathrm{bulge}}$ relation measured in elliptical
galaxies and classical bulges by Kormendy & Ho (2013) (see Figure 7). This
picture may be influenced, however, by a number of both selection and
measurement biases.

###  4.2 Potential Biases in Host Galaxy Mass Measurements

####  4.2.1 Malmquist Bias in Quasar Sample

When searching for high-redshift quasars via wide-field galaxy surveys, the
most luminous quasars—with the most active black holes—are more readily
observed; in addition, due to the rapid decline of the galaxy luminosity
function at high masses, these active, massive black holes are more often
found in more common, lower-luminosity, lower-mass galaxies (where they are
outliers) rather than rarer, bright, higher-mass galaxies (as discussed in
Lauer et al., 2007). It has therefore been difficult to say with any certainty
whether the observed tendency of high-redshift quasars to fall above the local
$M_{\mathrm{BH}}-M_{*}$ relation, as J2239+0207 does, is evidence for
evolution of $M_{\mathrm{BH}}/M_{*}$ with cosmic time, or simply a consequence
of selection bias.

To evaluate this bias, we compare our derived mass for J2239+0207 and other
measurements of $M_{\mathrm{BH}}/M_{*}$ at $z\sim 6$ to a local sample with
the same biases. For 30 extremely luminous Palomar Green (PG) quasars (with
$L_{\mathrm{AGN}}\geq 2\times 10^{12}\,L_{\odot}$, at $z\leq 0.5$), we
retrieve black hole masses obtained via single-epoch estimations (i.e. using
line widths and luminosities) from the literature (Vestergaard & Peterson,
2006; Veilleux et al., 2009; Yu et al., 2020). In addition, we adopt total
galaxy masses estimated from the quasar-galaxy image decomposition (Zhang et
al., 2016). These considerations ensure that the black hole and galaxy masses
for both this local sample and J2239+0207 are measured by similar methods. We
plot these PG quasars alongside the $z\sim 6$ samples and J2239+0207 in Figure
7. By using the most luminous PG sources, we recreate the bias towards high
$L_{\mathrm{AGN}}$ present in high-redshift quasar surveys; this sample choice
also ensures that both samples are subject to the Malmquist-type bias in
$M_{\mathrm{BH}}$ measurements derived from single-epoch spectra discussed by
Shen & Kelly (2010).

Despite both samples being similarly affected by the Lauer et al. (2007) and
Shen & Kelly (2010) biases, the high-$L_{\mathrm{AGN}}$ PG quasars and $z\sim
6$ quasars still appear to follow different distributions in Figure 7. On
average, the local quasars lie in more massive host galaxies at a given
$M_{\mathrm{BH}}$ (the PG quasars are also largely consistent with the
Kormendy & Ho (2013) local $M_{\mathrm{BH}}-M_{\mathrm{bulge}}$ relation). The
$z\sim 6$ quasars show much larger scatter, overlapping with the local
examples but extending to a factor of $\sim$ 10 above them.

####  4.2.2 Use of Dynamical Masses of Host Galaxies

J2239+0207 presents an example of a severely mismatched gas dynamical mass and
stellar mass in a high-$z$ quasar. As JWST observations create larger samples
of $z\gtrsim 6$ quasar host galaxy stellar masses, is a trend likely to
emerge? Do we expect [C II] gas dynamical masses to lie systematically above
masses measured from stellar properties?

This appears not to be an issue up to moderate redshifts. In these cases, the
dynamical masses come primarily from stellar spectroscopy, with some from
MASER observations, and agree very well with stellar masses as measured in the
$K_{S}$ band (Kormendy & Ho, 2013). Using either indicator for stellar mass,
there is a tight relation between galaxy bulge mass and SMBH mass among non-
AGN galaxies, often termed the ‘Magorrian Relation’ (Kormendy & Ho, 2013).
Many quasars have reliable stellar masses based on photometric and dynamical
data; for $z\leq 2$, their host galaxies also follow the local Magorrian
Relation (Kormendy & Ho, 2013).

However, this is a general issue at high redshift. Studies of the pre-JWST
high-$z$ quasar samples (including Izumi et al., 2021; Pensabene et al., 2020,
gray points in Figure 7) compare the galaxy dynamical mass determined mostly
from the [C II] 158 $\mu$m line measured with ALMA with the local relation
between galaxy bulge and black hole masses. As shown in Figure 7 and many
other works, at high redshift there is a large scatter relative to the local
relation, with many galaxies indicated to have SMBH masses substantially in
excess of its prediction. Is this effect an indication of a real change in the
behavior of galaxies or is it a product of using the [C II]-based dynamical
masses rather than the approaches used at lower redshift?

If the inferred dynamical masses tend to strongly underestimate stellar
masses, the observed differences between local and high-$z$ quasars in Figure
7 may simply stem from the resulting bias. On the other hand, if the dynamical
masses overestimate the stellar mass, the trend in Figure 7 is reinforced. The
usual formalism (e.g., Izumi et al., 2019) to translate line widths to
dynamical masses is summarized in equation 3:

$\frac{M_{\mathrm{dyn}}}{M_{\odot}}=1.16\times
10^{5}\left(\frac{v_{\mathrm{circ}}}{\mathrm{km~{}s}^{-1}}\right)^{2}~{}\left(\frac{D}{\mathrm{kpc}}\right)$
(3)

where $D$ is the diameter of the galaxy and

$v_{\mathrm{circ}}=0.75\frac{\mathrm{FWHM}}{\sin i}$ (4)

with $i$ the angle of inclination, which is typically deduced from the ratio
of the major and minor axes. The potential for overestimation lies primarily
in the assumption that the line FWHM is dominated by the rotation of a disk
galaxy. It is possible that non-Keplerian motions, e.g. winds, may be further
broadening the line profile and biasing the dynamical mass higher (see Izumi
et al., 2019). Another possibility is that spheroidal-dominated galaxies with
axial ratios close to unity will instead be assigned small inclinations
appropriate for face-on disks, leading to a large and incorrect upward
adjustment in $v_{\mathrm{circ}}$. The potential for underestimation lies in
assigning too small an inclination, $i$, or diameter, $D$. A detailed
theoretical simulation carried out by Lupi et al. (2019) found that the [C II]
dynamical mass underestimates the total stellar mass of a simulated
rotationally supported galaxy at $z\sim 7$ by a factor of $2-5$, and
potentially by more for a dispersion-dominated case. As they conclude,
corrections of this size in the [C II]-derived masses would substantially
reduce the difference between local and high redshift systems shown in Figure
7.

However, the dynamical mass for J2239+0207 from [CII] observations is
$2.9\times 10^{11}$ M⊙ (Izumi et al., 2019); contrary to the models of Lupi et
al. (2019), this is nearly 30 times larger than our inferred stellar mass,
despite the latter being determined in the rest-optical where the effects of
extinction and recent star formation are minimized. At least in this case, it
seems that the effects that can increase the estimated dynamical mass – winds,
incorrect values for inclination and diameter – may be dominating the [C II]
line width interpretation.

We will evaluate this issue as follows. We will compare the dynamical masses
based on submm- and mm-wave line measures with those from photometry for the
high $L_{AGN}$ PG quasars that are analogous to the high redshift quasars.
This should show if the two methods yield consistent results, as is an
implicit assumption in the comparison of the Magorrian ratios at high redshift
with the local relation.

We use CO line measurements of the PG quasar sample in place of [C II]. The
velocity widths of the [C II] line and the CO lines are very similar in local
galaxies (Röllig et al., 2012; de Blok et al., 2016), so we take the CO
measurements as a proxy for [C II] in the PG quasars. All the available CO
measurements of PG quasars that also have direct stellar mass determinations
(Zhang et al., 2016) are listed in Table 2. We use equation (3) to translate
the line widths to masses, taking the diameter of the host from Guyon et al.
(2006); Veilleux et al. (2009), with the former used for total diameters and
the latter for bulges (the distinction between total and bulge is undetermined
for high redshift host galaxies). The axial ratios for the PG quasars are
generally close to 1 - the average for all but two where they are not
available is 0.76. This would imply small inclinations and large corrections
to the mass, but to be conservative we have set $\sin i$ = 1, and consequently
our mass estimates are lower limits. Since our concern is that the masses
might be too low, this does not undermine our analysis.

As shown in Table 2, assuming bulge diameters there are 3/10 cases where the
dynamical mass lower limit is more than two times less than the directly
determined stellar mass, and 3/10 cases where it is more than two times
greater. When the dynamical masses are calculated using the total diameters of
the galaxies, all ten cases have dynamical mass lower limits greater than the
directly determined stellar masses. That is, the behavior of the PG quasars
does not support the hypothesis that the line-widths-approach results in a
significant and uniform underestimate of host galaxy masses.

This analysis may point to a genuine evolution in $M_{\mathrm{BH}}/M_{*}$ with
cosmic time, rather than an apparent evolution solely due to biases. However,
the scatter in Figure 7 suggests that any such evolution will not appear
uniformly. Instead, the $M_{\mathrm{BH}}$-$M_{*}$ ratio is likely to take on a
much larger range at high redshift than occurs locally. Many more quasar host
galaxies will be directly measured with JWST in the coming years, and will
help probe this possibility.

Figure 7: The $M_{\mathrm{BH}}-M_{*}$ relation for J2239+0207 (blue star) and a selection of other $z\sim 6$ and local quasars from the literature. $z\sim 6$ quasar dynamical masses are shown in grey (Izumi et al. (2021), diamonds; Pensabene et al. (2020), circles), while high-$z$ quasar stellar masses from recent JWST samples (Kocevski et al. (2023), squares; Harikane et al. (2023), pentagons) are shown in blue. Luminous Palomar Green (PG) quasars at $z<0.5$ with $L_{AGN}>10^{12}L_{\odot}$ are shown as green triangles. The local $M_{\mathrm{BH}}-M_{\mathrm{bulge}}$ relation in local non-AGN systems from Kormendy & Ho (2013) (grey shaded region) and the corresponding relation for the total stellar mass (grey dashed line, masses scaled in proportion to $M_{K}$) are shown. Two vertical arrows indicate the Malmquist bias we find for the luminous PG quasars. We also show the $M_{\mathrm{BH}}-M_{*}$ relation for low-mass AGN host galaxies from Reines & Volonteri (2015) as a grey solid line. PG quasars and $z\sim 6$ quasars appear to follow different distributions, with low-z quasars containing more massive host galaxies for a given black hole mass; however, a number of biases may affect this conclusion. Table 2: Comparison of Mass Determinations ID | CO FWHM | refaa1 – Shangguan et al. (2020): 2 – Evans et al. (2001): 3 – Evans et al. (2009) | log(M∗/M⊙)bbFrom stellar photometry, Zhang et al. (2016). | r1/2,bulgeccFrom Veilleux et al. (2009), Table 5 | log(Mdyn/M⊙)ddBased on computing the dynamical mass for rbulge. These values are lower limits because sin$i$ has been set to 1.  | Mdyn/M∗eeLower limits, based on masses for rbulge. | rtotalffThe larger of the bulge and disk components fitted by Guyon et al. (2006), table 6. | log(Mdyn/M⊙)ggLower limit because sin$i$ is set to 1, based on dynamical mass computed for rtotal. | Mdyn/M∗hhLower limits, based on masses for rtotal.
---|---|---|---|---|---|---|---|---|---
| (km s-1) |  |  | (kpc) | (lower limit) | (lower limit) | (kpc) | (lower limit) | (lower limit)
PG 0007+106 | 387 | 1 | 11.02 | 2.97 | 11.76 | 5.54 |  |  |
PG 0050+124 | 378 | 1 | 11.3 | 1.03 | 11.28 | 0.96 |  |  |
PG 0838+770 | 90 | 2 | 11.32 | 0.56 | 9.77 | 0.03 |  |  |
PG 0923+129 | 362 | 1 | 11.27 |  |  |  |  |  |
PG 1119+120 | 213 | 1 | 10.85 | 0.45 | 10.42 | 0.38 | 1.47 | 10.94 | 1.23
PG 1126-041 | 467 | 1 | 11.03 | 0.79 | 11.35 | 2.10 | 3.08 | 11.94 | 8.19
PG 1211+143 | 66 | 1 | 10.56 |  |  |  | 11.50 | 10.81 | 1.79
PG 1229+204 | 202 | 1 | 11.12 | 3.65 | 11.29 | 1.48 | 4.16 | 11.35 | 1.68
PG 1426+015 | 344 | 1 | 11.23 | 1.64 | 11.40 | 1.49 | 8.70 | 12.13 | 7.90
PG 1440+356 | 370 | 2 | 11.23 | 0.39 | 10.84 | 0.41 | 1.53 | 11.44 | 1.61
PG 1613+658 | 490 | 2 | 11.64 |  |  |  | 15.23 | 12.68 | 10.93
PG 1700+518 | 260 | 3 | 11.57 |  |  |  | 18.78 | 12.22 | 4.46
PG 2130+099 | 548 | 1 | 11.03 | 2.83 | 12.05 | 10.36 | 1.53 | 11.78 | 5.61
PG 2214+139 | 180 | 1 | 11.16 | 2.78 | 11.07 | 0.81 | 3.45 | 11.16 | 1.01

##  5 Summary and Conclusions

We utilize dedicated observations of a PSF star in three JWST/NIRCam medium
bands, as well as a NIRSpec prism spectrum, to remove the contribution from
the central point source of the quasar J2239+0207 and characterize its host
galaxy.

  * •

Using the broad H$\alpha$ line in the NIRSpec spectrum, we obtain a black hole
mass of $3.4\times 10^{8}$ M⊙, somewhat lower than previously estimated values
from Mg II and [C IV].

  * •

We subtract the PSF in our three bands and detect emission from the host
galaxy in the two filters longward of the Balmer break, measuring host galaxy
fluxes of $0.34\;\mu$Jy in F360M and $0.44\;\mu$Jy in F480M (host-to-quasar
flux ratios of 4.6% and 3.9% respectively). We relate the F360M flux to galaxy
mass using IRAC Band 1 observations of galaxies from $5<z<8$ in the
literature, obtaining a stellar mass of approximately $10^{10}$ M⊙, which is
significantly less than the existing dynamical mass measurement from 158
$\mu$m [C II]. J2239+0207 lies above the local $M_{\mathrm{BH}}-M_{*}$
relation, and is among the lower-mass host galaxies explored at $z\sim 6$.

  * •

To explore the $M_{\mathrm{BH}}-M_{*}$ relation of J2239+0207 and other $z\sim
6$ quasars and its deviation from the local Universe, we utilize a sample of
high-luminosity PG quasars with masses from single-epoch spectra and total
host galaxy masses derived from photometry. This selection ensures that both
of our samples are subject to the Lauer et al. (2007) bias, and any evolution
with redshift is not a result of this bias. The host of J2239+0207 is
significantly undermassive compared with typical local quasar hosts.

  * •

At a given $M_{\mathrm{BH}}$, the $z\sim 6$ quasars from the literature appear
to show a much larger range of host galaxy masses than found in local cases.
Most of the $z\sim 6$ quasars use [C II] dynamical mass as a proxy for stellar
mass. We evaluate the determination of host galaxy masses through this method
by comparing with similar measurements using CO lines from PG quasars. We find
no evidence for systematic underestimates that might result in the relatively
low host masses deduced for many high redshift systems. We hypothesize that
the $M_{\mathrm{BH}}$-$M_{*}$ ratio will possess more scatter at higher
redshift than in the local Universe, but larger samples are required to
confirm this behavior.

JWST’s wavelength range and sensitivity will for the first time allow the
detection and characterization of large numbers of high-redshift quasar hosts
in the rest-frame optical, using NIRCam photometry. The method of PSF
subtraction using observations of a dedicated PSF star is able to retrieve
galaxy signals even at the $<5$% host-to-quasar level. As stellar masses are
obtained and compared to [C II] dynamical masses for more high-redshift
quasars and local and high-redshift samples can be compared consistently, our
analysis of J2239+0207 demonstrates that JWST will allow for the robust
characterization of any evolution of the $M_{\mathrm{BH}}-M_{*}$ relation over
cosmic time.

We thank Tod Lauer for important discussions regarding selection biases and
Feige Wang, Jinyi Yang, and Xiaohui Fan for their insightful comments. The
anonymous referee’s feedback also greatly improved this work. MS, SA, JL, and
GR acknowledge support from the JWST Mid-Infrared Instrument (MIRI) grant
80NSSC18K0555, and the NIRCam science support contract NAS5-02105, both from
NASA Goddard Space Flight Center to the University of Arizona. The JWST data
presented in this paper were obtained from the Mikulski Archive for Space
Telescopes (MAST) at the Space Telescope Science Institute. The specific
observations analyzed can be accessed via http://dx.doi.org/10.17909/ajqw-7w52
(catalog DOI.)

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bradley et al. (2022) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2022, astropy/photutils: 1.6.0, 1.6.0, Zenodo, Zenodo, doi: 10.5281/zenodo.7419741
  * Coatman et al. (2017) Coatman, L., Hewett, P. C., Banerji, M., et al. 2017, MNRAS, 465, 2120, doi: 10.1093/mnras/stw2797
  * Curtis-Lake et al. (2013) Curtis-Lake, E., McLure, R. J., Dunlop, J. S., et al. 2013, MNRAS, 429, 302, doi: 10.1093/mnras/sts338
  * de Blok et al. (2016) de Blok, W. J. G., Walter, F., Smith, J. D. T., et al. 2016, AJ, 152, 51, doi: 10.3847/0004-6256/152/2/51
  * Decarli et al. (2018) Decarli, R., Walter, F., Venemans, B. P., et al. 2018, ApJ, 854, 97, doi: 10.3847/1538-4357/aaa5aa
  * Ding et al. (2022) Ding, X., Onoue, M., Silverman, J. D., et al. 2022, arXiv e-prints, arXiv:2211.14329, doi: 10.48550/arXiv.2211.14329
  * Dunlop et al. (2003) Dunlop, J. S., McLure, R. J., Kukula, M. J., et al. 2003, MNRAS, 340, 1095, doi: 10.1046/j.1365-8711.2003.06333.x
  * Elvis et al. (1994) Elvis, M., Wilkes, B. J., McDowell, J. C., et al. 1994, ApJS, 95, 1, doi: 10.1086/192093
  * Evans et al. (2001) Evans, A. S., Frayer, D. T., Surace, J. A., & Sanders, D. B. 2001, AJ, 121, 1893, doi: 10.1086/319972
  * Evans et al. (2009) Evans, A. S., Hines, D. C., Barthel, P., et al. 2009, AJ, 138, 262, doi: 10.1088/0004-6256/138/1/262
  * Farina et al. (2022) Farina, E. P., Schindler, J.-T., Walter, F., et al. 2022, ApJ, 941, 106, doi: 10.3847/1538-4357/ac9626
  * Gehren et al. (1984) Gehren, T., Fried, J., Wehinger, P. A., & Wyckoff, S. 1984, ApJ, 278, 11, doi: 10.1086/161763
  * Greene & Ho (2005) Greene, J. E., & Ho, L. C. 2005, ApJ, 630, 122, doi: 10.1086/431897
  * Guyon et al. (2006) Guyon, O., Sanders, D. B., & Stockton, A. 2006, ApJS, 166, 89, doi: 10.1086/505030
  * Harikane et al. (2023) Harikane, Y., Zhang, Y., Nakajima, K., et al. 2023, arXiv e-prints, arXiv:2303.11946, doi: 10.48550/arXiv.2303.11946
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Hutchings et al. (2002) Hutchings, J. B., Frenette, D., Hanisch, R., et al. 2002, AJ, 123, 2936, doi: 10.1086/340472
  * Izumi et al. (2019) Izumi, T., Onoue, M., Matsuoka, Y., et al. 2019, PASJ, 71, 111, doi: 10.1093/pasj/psz096
  * Izumi et al. (2021) Izumi, T., Matsuoka, Y., Fujimoto, S., et al. 2021, ApJ, 914, 36, doi: 10.3847/1538-4357/abf6dc
  * Jahnke et al. (2004) Jahnke, K., Sánchez, S. F., Wisotzki, L., et al. 2004, ApJ, 614, 568, doi: 10.1086/423233
  * Kennicutt & Evans (2012) Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531, doi: 10.1146/annurev-astro-081811-125610
  * Kim et al. (2017) Kim, M., Ho, L. C., Peng, C. Y., Barth, A. J., & Im, M. 2017, ApJS, 232, 21, doi: 10.3847/1538-4365/aa8a75
  * Kocevski et al. (2023) Kocevski, D. D., Onoue, M., Inayoshi, K., et al. 2023, arXiv e-prints, arXiv:2302.00012, doi: 10.48550/arXiv.2302.00012
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Labbé et al. (2010) Labbé, I., González, V., Bouwens, R. J., et al. 2010, ApJ, 716, L103, doi: 10.1088/2041-8205/716/2/L103
  * Lauer et al. (2007) Lauer, T. R., Tremaine, S., Richstone, D., & Faber, S. M. 2007, ApJ, 670, 249, doi: 10.1086/522083
  * Li et al. (2021) Li, J., Silverman, J. D., Ding, X., et al. 2021, ApJ, 918, 22, doi: 10.3847/1538-4357/ac06a8
  * Lupi et al. (2019) Lupi, A., Volonteri, M., Decarli, R., et al. 2019, MNRAS, 488, 4004, doi: 10.1093/mnras/stz1959
  * Lyu et al. (2016) Lyu, J., Rieke, G. H., & Alberts, S. 2016, ApJ, 816, 85, doi: 10.3847/0004-637X/816/2/85
  * Marshall et al. (2020) Marshall, M. A., Mechtley, M., Windhorst, R. A., et al. 2020, ApJ, 900, 21, doi: 10.3847/1538-4357/abaa4c
  * Marshall et al. (2023) Marshall, M. A., Perna, M., Willott, C. J., et al. 2023, arXiv e-prints, arXiv:2302.04795, doi: 10.48550/arXiv.2302.04795
  * McLeod & Bechtold (2009) McLeod, K. K., & Bechtold, J. 2009, ApJ, 704, 415, doi: 10.1088/0004-637X/704/1/415
  * McLure et al. (2011) McLure, R. J., Dunlop, J. S., de Ravel, L., et al. 2011, MNRAS, 418, 2074, doi: 10.1111/j.1365-2966.2011.19626.x
  * Mechtley et al. (2012) Mechtley, M., Windhorst, R. A., Ryan, R. E., et al. 2012, ApJ, 756, L38, doi: 10.1088/2041-8205/756/2/L38
  * Mollá et al. (2009) Mollá, M., García-Vargas, M. L., & Bressan, A. 2009, MNRAS, 398, 451, doi: 10.1111/j.1365-2966.2009.15160.x
  * Neeleman et al. (2021) Neeleman, M., Novak, M., Venemans, B. P., et al. 2021, ApJ, 911, 141, doi: 10.3847/1538-4357/abe70f
  * Onoue et al. (2019) Onoue, M., Kashikawa, N., Matsuoka, Y., et al. 2019, ApJ, 880, 77, doi: 10.3847/1538-4357/ab29e9
  * Peng et al. (2006) Peng, C. Y., Impey, C. D., Rix, H.-W., et al. 2006, ApJ, 649, 616, doi: 10.1086/506266
  * Pensabene et al. (2020) Pensabene, A., Carniani, S., Perna, M., et al. 2020, A&A, 637, A84, doi: 10.1051/0004-6361/201936634
  * Perrin et al. (2014) Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, ed. J. Oschmann, Jacobus M., M. Clampin, G. G. Fazio, & H. A. MacEwen, 91433X, doi: 10.1117/12.2056689
  * Reines & Volonteri (2015) Reines, A. E., & Volonteri, M. 2015, ApJ, 813, 82, doi: 10.1088/0004-637X/813/2/82
  * Rieke et al. (2017) Rieke, G., Alberts, S., Lyu, J., & Shivaei, I. 2017, Formation Histories and Stellar Masses of Very High-z Quasars, JWST Proposal. Cycle 1, ID. #1205
  * Rieke et al. (2005) Rieke, M. J., Kelly, D., & Horner, S. 2005, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5904, Cryogenic Optical Systems and Instruments XI, ed. J. B. Heaney & L. G. Burriesci, 1–8, doi: 10.1117/12.615554
  * Rieke et al. (2023) Rieke, M. J., Kelly, D. M., Misselt, K., et al. 2023, PASP, 135, 028001, doi: 10.1088/1538-3873/acac53
  * Röllig et al. (2012) Röllig, M., Simon, R., Güsten, R., et al. 2012, A&A, 542, L22, doi: 10.1051/0004-6361/201218935
  * Schneider et al. (2010) Schneider, D. P., Richards, G. T., Hall, P. B., et al. 2010, AJ, 139, 2360, doi: 10.1088/0004-6256/139/6/2360
  * Schramm et al. (2019) Schramm, M., Rujopakarn, W., Silverman, J. D., et al. 2019, ApJ, 881, 145, doi: 10.3847/1538-4357/ab2cd3
  * Shangguan et al. (2020) Shangguan, J., Ho, L. C., Bauer, F. E., Wang, R., & Treister, E. 2020, ApJS, 247, 15, doi: 10.3847/1538-4365/ab5db2
  * Shao et al. (2017) Shao, Y., Wang, R., Jones, G. C., et al. 2017, ApJ, 845, 138, doi: 10.3847/1538-4357/aa826c
  * Shao et al. (2022) Shao, Y., Wang, R., Weiss, A., et al. 2022, A&A, 668, A121, doi: 10.1051/0004-6361/202244610
  * Shen & Kelly (2010) Shen, Y., & Kelly, B. C. 2010, ApJ, 713, 41, doi: 10.1088/0004-637X/713/1/41
  * Stanway et al. (2016) Stanway, E. R., Eldridge, J. J., & Becker, G. D. 2016, MNRAS, 456, 485, doi: 10.1093/mnras/stv2661
  * Targett et al. (2012) Targett, T. A., Dunlop, J. S., & McLure, R. J. 2012, MNRAS, 420, 3621, doi: 10.1111/j.1365-2966.2011.20286.x
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Veilleux et al. (2009) Veilleux, S., Kim, D. C., Rupke, D. S. N., et al. 2009, ApJ, 701, 587, doi: 10.1088/0004-637X/701/1/587
  * Venemans et al. (2016) Venemans, B. P., Walter, F., Zschaechner, L., et al. 2016, ApJ, 816, 37, doi: 10.3847/0004-637X/816/1/37
  * Vestergaard & Peterson (2006) Vestergaard, M., & Peterson, B. M. 2006, ApJ, 641, 689, doi: 10.1086/500572
  * Walter et al. (2004) Walter, F., Carilli, C., Bertoldi, F., et al. 2004, ApJ, 615, L17, doi: 10.1086/426017
  * Wang et al. (2013) Wang, R., Wagg, J., Carilli, C. L., et al. 2013, ApJ, 773, 44, doi: 10.1088/0004-637X/773/1/44
  * Wolff et al. (2023) Wolff, S. G., Gáspár, A., Rieke, G. H., Ballering, N., & Ygouf, M. 2023, arXiv e-prints, arXiv:2302.10085, doi: 10.48550/arXiv.2302.10085
  * Yu et al. (2020) Yu, L.-M., Bian, W.-H., Zhang, X.-G., et al. 2020, ApJ, 901, 133, doi: 10.3847/1538-4357/abb01e
  * Yue et al. (2018) Yue, M., Jiang, L., Shen, Y., et al. 2018, ApJ, 863, 21, doi: 10.3847/1538-4357/aacf04
  * Yue et al. (2021) Yue, M., Yang, J., Fan, X., et al. 2021, ApJ, 917, 99, doi: 10.3847/1538-4357/ac0af4
  * Zhang et al. (2023) Zhang, H., Behroozi, P., Volonteri, M., et al. 2023, arXiv e-prints, arXiv:2303.08150, doi: 10.48550/arXiv.2303.08150
  * Zhang et al. (2016) Zhang, Z., Shi, Y., Rieke, G. H., et al. 2016, ApJ, 819, L27, doi: 10.3847/2041-8205/819/2/L27