# Detection of stellar light from quasar host galaxies at redshifts above 6

Xuheng Ding, Masafusa Onoue, John Silverman, Yoshiki Matsuoka, Takuma
Izumi, Michael Strauss, Knud Jahnke, Camryn Phillips, Junyao Li, Marta
Volonteri, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

Xuheng Ding, Masafusa Onoue, John Silverman, Yoshiki Matsuoka, Takuma Izumi, et al. Detection of stellar light from quasar host galaxies at redshifts above 6. Nature, 2023, 10.1038/s41586-023-06345-5. insu-04154889

HAL Id: insu-04154889

https://insu.hal.science/insu-04154889v1

Submitted on 3 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
University of Groningen

# Detection of stellar light from quasar host galaxies at redshifts above 6

Ding, Xuheng; Onoue, Masafusa; Silverman, John D.; Matsuoka, Yoshiki; Izumi, Takuma; Strauss, Michael A.; Jahnke, Knud; Phillips, Camryn L.; Li, Junyao; Volonteri, Marta

Published in:
Nature

DOI:
10.1038/s41586-023-06345-5

**IMPORTANT NOTE:** You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from it. Please check the document version below.

Document Version
Publisher's PDF, also known as Version of record

Publication date:
2023

Link to publication in University of Groningen/UMCG research database

## Citation for published version (APA):

Ding, X., Onoue, M., Silverman, J. D., Matsuoka, Y., Izumi, T., Strauss, M. A., Jahnke, K., Phillips, C. L., Li, J., Volonteri, M., Haiman, Z., Andika, I. T., Aoki, K., Baba, S., Bieri, R., Bosman, S. E. I., Bottrell, C., Eilers, A. C., Fujimoto, S., ... Yang, J. (2023). Detection of stellar light from quasar host galaxies at redshifts above 6. *Nature*, **621**(7977), 51-55. https://doi.org/10.1038/s41586-023-06345-5

## Copyright

Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).

The publication may also be distributed here under the terms of Article 25fa of the Dutch Copyright Act, indicated by the "Taverne" license. More information can be found on the University of Groningen website: https://www.rug.nl/library/open-access/self-archiving-pure/taverne-amendment.

## Take-down policy

If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.

Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the number of authors shown on this cover page is limited to 10 maximum.
# Detection of stellar light from quasar host galaxies at redshifts above 6

https://doi.org/10.1038/s41586-023-06345-5

Received: 25 November 2022

Accepted: 20 June 2023

Published online: 28 June 2023

<figure><img src="image_3.png" /></figure>

Check for updates

Xuheng Ding¹,²,³, Masafusa Onoue¹,³,⁴,⁵, John D. Silverman¹,²,⁵, Yoshiki Matsuoka⁶,
Takuma Izumi⁷,⁸, Michael A. Strauss⁹, Knud Jahnke⁴, Camryn L. Phillips⁹, Junyao Li¹⁰,
Marta Volonteri¹¹, Zoltan Haiman¹²,¹³, Irham Taufik Andika¹⁴,¹⁵, Kentaro Aoki¹⁶,
Shunsuke Baba¹⁷, Rebekka Bier¹⁸, Sarah E. I. Bosman⁴, Connor Bottrell¹,²,
Anna-Christina Eilers¹⁹, Seiji Fujimoto²⁰, Melanie Habouzit⁴,²¹, Masatoshi Imanishi⁷,²²,
Kohei Inayoshi³, Kazushi Iwasawa²³,²⁴, Nobunari Kashikawa⁵,²⁵, Toshihiro Kawaguchi²⁶,
Kotaro Kohno²⁵,²⁷, Chien-Hsiu Lee²⁸, Alessandro Lupi²⁹, Jianwei Lyu³⁰, Tohru Nagao⁶,
Roderik Overzier³¹, Jan-Torge Schindler³², Malte Schramm³³, Kazuhiro Shimasaku⁵,²⁵,
Yoshiki Toba⁷,³⁴, Benny Trakhtenbrot³⁵, Maxime Trebitsch³⁶, Tommaso Treu³⁷,
Hideki Umehata³⁸,³⁹, Bram P. Venemans³², Marianne Vestergaard³⁰,⁴⁰, Fabian Walter⁴,
Feige Wang³⁰ & Jinyi Yang³⁰

The detection of starlight from the host galaxies of quasars during the reionization epoch ($z > 6$) has been elusive, even with deep Hubble Space Telescope observations¹,². The current highest redshift quasar host detected³, at $z = 4.5$, required the magnifying effect of a foreground lensing galaxy. Low-luminosity quasars⁴⁻⁶ from the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP)⁷ mitigate the challenge of detecting their underlying, previously undetected host galaxies. Here we report rest-frame optical images and spectroscopy of two HSC-SSP quasars at $z > 6$ with the JWST. Using near-infrared camera imaging at 3.6 and 1.5 µm and subtracting the light from the unresolved quasars, we find that the host galaxies are massive (stellar masses of $13 \times 3.4 \times 10^{10} M_{\odot}$, respectively), compact and disc-like. Near-infrared spectroscopy at medium resolution shows stellar absorption lines in the more massive quasar, confirming the detection of the host. Velocity-broadened gas in the vicinity of these quasars enables measurements of their black hole masses ($1.4 \times 10^9$ and $2.0 \times 10^8 M_{\odot}$, respectively). Their location in the black hole mass–stellar mass plane is consistent with the distribution at low redshift, suggesting that the relation between black holes and their host galaxies was already in place less than a billion years after the Big Bang.

The quasars J2236+0032 and J2255+0251 (coordinates listed in Extended Data Table 1) are at redshifts of 6.40 and 6.34. They have relatively low luminosity (absolute magnitude at 1450 Å, $M_{1450} = -23.87$ and $-23.66$ mag, respectively), roughly ten times less luminous than typical bright quasars at these redshifts⁸⁻¹⁰. Using the JWST with the Near-Infrared Camera (NIRCam)¹¹, we observed J2236+0032 on 6 November 2022 and J2255+0251 on 26 October 2022. Images were taken in two broad-band filters (F356W and F150W) to detect their host galaxies at 3.56 and 1.50 µm, thus bracketing the rest-frame 4,000 Å break, a characteristic feature of galaxy spectra that constrains the age of the stellar population. Total exposure times were roughly 3,100 s for each filter per target.

In Fig. 1, we show a 24 × 24" region and zoomed-in cut-outs of the JWST images centred on the quasars for both filters. As expected, both quasars are detected at high significance. The first quasar J2236+0032 (magnitude in filters F356W and F150W, respectively, $m_{F356W} = 21.75 \pm 0.02$ mag, $m_{F150W} = 22.78 \pm 0.02$) has no evidence of close neighbours. The second quasar J2255+0251 ($m_{F356W} = 22.20 \pm 0.02$ mag, $m_{F150W} = 23.03 \pm 0.02$) has two neighbouring sources (less than 2.6")

that do not have spectroscopic redshifts at this point; however, their detection in the Hyper Suprime-Cam (HSC) *r* and *i* bands suggests that both of them are at lower redshifts than the central quasar.

The detection of an underlying host galaxy requires careful decomposition of the two-dimensional (2D) light distribution to separate the unresolved quasar from its host galaxy. We perform a two-component model fit of the rest-frame optical emission in a region of the NIRCam images centred on each target. We use the software galight¹² to identify stars from each image to model the 2D point-spread function (PSF). We then model each quasar image in each filter as a sum of a point-like quasar and an extended host. The host galaxy was modelled using a 2D Sérsic profile¹³ whose parameters include the position, size ($R_{eff}$), axis ratio, orientation, flux and the radial profile and/or shape of the light distribution (that is, Sérsic index), allowing a flexible description of the a priori unknown host galaxy shape. In Extended Data Table 1, we report the best-fit parameters and their uncertainties on the basis of the dispersion resulting from different PSF models with different modelling techniques (Methods).
<figure><img src="image_4.png" /><figcaption>**Fig. 1** | NIRCam images of the field around J2236+0032 and J2255+0251. Both JWST filters F356W (left panels) and F150W (right panels) are shown. The middle panels are zoomed-in regions of each image. One of the PSF stars used for 2D image decomposition is indicated in each quasar field.</figcaption></figure>

As shown in Fig. 2, we detect the host galaxies for both targets, which are visible after subtracting the quasar point-source model (third column). For J2236+0032, the quasar host is significantly detected in both filters ($m_{F356W}$ = 23.12 ± 0.20; $m_{F150W}$ = 25.12 ± 0.29 mag). The inferred host-to-total flux ratios are 25.5 ± 4.4 and 10.2 ± 2.8% with sizes ($R_{eff}$) of 0.7 ± 0.1 and 0.5 ± 0.2 kpc for F356W and F150W, respectively. The light distribution is elongated in the same direction (that is, NE-SW) in the two filters, suggesting a mildly inclined disc-like host. For this target, we fix the Sérsic index to 1 in both bands—corresponding to an exponential disc model—to avoid non-physical fits; this is supported by a subsequent fit to the host-galaxy emission after subtracting the quasar (Methods). For J2255+0251, the host galaxy is clearly detected in the F356W band ($m_{F356W}$ = 24.58 ± 0.30 mag) with a host-to-total flux ratio of 9.8 ± 2.6%. The host emission extends southeast (left) of the quasar’s location. The size of the host is inferred to be $R_{eff}$ = 1.5 ± 1.1 kpc with a Sérsic index of 1.5 ± 1.4, suggesting a more disc-like than spheroidal profile. The host galaxy is not detected in the F150W image, resulting in an upper limit of $m_{F150W}$ = 26.3 mag (Methods). This upper limit provides an important constraint on the galaxy spectral energy distribution (SED). A non-detection is not surprising, given the low host-to-total flux ratio in the F356W band and the fact that galaxies are typically fainter bluewards of the 4,000 Å break. We use NIRSpec¹⁴ observations to demonstrate that the contribution of extended line emission to the host photometry in F356W is minimal (Methods).

Our two photometric bands lie on either side of the 4,000 Å break, allowing us to quantify the contributions of young and relatively old stars and thus to estimate the stellar masses of the two quasar host galaxies. We fit the SED of each host with a single stellar population,

in which the stellar metallicity (the ratio of the host galaxies’ metallicity to the solar metallicity in log space) $logZ/Z_⊙$ is constrained to the range [-1, -0.3], the stellar age is within the range [0.01, 0.84] Gyr (capped by the age of the Universe at redshift 6.4) and the extinction $A_ν$ is in the range of [0, 5.0] mag. With standard assumptions on the stellar initial mass function, we find that J2236+0032 has a stellar mass of $1.3_{-0.6}^{+2.0} \times 10^{11} M_⊙$, whereas J2255+0251 is less massive with $3.4_{-1.9}^{+7.6} \times 10^{10} M_⊙$. We use simulations and three different SED fitting codes to demonstrate that the results are robust, although further observations will lessen potential systematic errors in these mass estimates (see Methods for full details on our mass estimates and their uncertainties). Our host galaxies are among the highest stellar mass galaxies known at $z \simeq 6$ (ref. 15).

A rest-frame optical spectrum of each quasar was obtained by the NIRSpec instrument on JWST (Fig. 3). These spectra cover 2.87–5.27 µm, or roughly 4,000–7,000 Å in the rest-frame, revealing strong emission lines including the Hβ + [O III] multiplet, Hα, and blended Fe II pseudo-continuum. The 0.2”-wide slit was aligned along the major axis of the host galaxy of J2236+0032 (Fig. 2). We serendipitously identify stellar absorption lines (Hγ, Hδ and He) in J2236+0032. Thus, the contribution of stellar light is confirmed both from the imaging and spectroscopy in this object. We use the velocity width of the broad Hβ emission line and the rest-frame 5,100 Å continuum luminosity to estimate the mass of the central black hole, under the assumption that the high-velocity gas is governed by the gravity of the supermassive black hole (SMBH). Using a recipe calibrated in the local Universe¹⁶ and subtracting the stellar emission, we obtain virial masses of black hole as $M_{BH}$ = $1.4 \times 10^9 M_⊙$ for J2236+0032 and $M_{BH}$ = $2.0 \times 10^8 M_⊙$ for J2255+0251, respectively (Methods).
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 2 | 2D decomposition into quasar and host-galaxy emission.** Panels are as follows from left to right: quasar image (data), best-fit model (quasar + host galaxy), host galaxy only (data – model quasar point-source component) and normalized residual image: that is, (data – model)/σ, where σ is the flux</figcaption></figure>

uncertainty of each pixel. The target name is above each row of panels and the filter is indicated on the left. The alignment of the NIRSpec slit (0.2 × 0.6") is shown in the left panel of the F356W image for each target. Pixel units are megajansky per steradian. A scale bar of 0.5" corresponds to 2.8 kpc at z ≈ 6.4.

Figure 4 shows the two $z \approx 6.4$ quasars in the black hole mass-stellar mass plane (left panel), as well as the inferred underlying relation compared to the local relation (right panel). As is observed in the local Universe, the quasar hosting the more massive black hole is hosted by the more massive galaxy. To evaluate their location in the black hole mass-stellar mass ($M_{BH}-M$) plane relative to the local mass relation, we carry out a simulation, following ref. 17, that takes into account the selection effects for these two quasars in the HSC survey. For this purpose, we combine the quasar flux limits, the estimated black hole mass and accretion rate functions at $z \approx 6$, as well as uncertainty in our stellar mass estimates to predict the expected $M_{BH}$ for these two quasars,

under the baseline assumption that the $M_{BH}-M$ relation does not evolve between $z=6.4$ and the present. In the left panel, the orange contours show the resulting bivariate probability distribution for the location of these two quasars, with the measured values overplotted in red. Our two $z>6$ quasars fall within the probability distribution (orange contours) of this simulated $z \approx 6.4$ sample and are thus consistent with no evolution of the $M_{BH}-M$ relation from $z \approx 6.4$ to $z=0$, as shown for the inferred underlying relation in the right panel of Fig. 4. With just two quasars, the probability distribution of the nominal evolution strength is broad as suggested by the width of these uncertainties (see also Extended Data Fig. 6 in the Methods). However, continuing
<figure><img src="image_7.png" /></figure>

**Fig. 3 | NIRSpec rest-frame optical spectra of J2236+0032 and J2255+0251.** For each quasar, the left panel shows the full spectrum and the top-right panels zoom in on the Hβ + [O III] λλ4960, 5008 doublet. The flux density Fλ is in units of 10−20 erg cm−2 s−1 Å−1. The data are shown in black with the errors per pixel in grey. The stellar absorption lines in J2236+0032 are identified in Hγ, Hδ and Hε,

which are also presented in the inset panel. Best-fit models of power-law continuum (blue), Fe II (green), broad Hβ (magenta), narrow Hβ + [O III] (orange thick line) and broad [O III] (orange thin line) are overplotted in the top-right panels. The sum of the continuum and emission line models is shown in red. The residuals are shown in the bottom-right panels.

JWST observations will give us a significantly larger sample, allowing us to better constrain models for the mutual evolution of the black hole and stellar populations in galaxies¹⁸⁻²⁰.

Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-023-06345-5.

1. Mechtley, M. et al. Near-infrared imaging of a z=6.42 quasar host galaxy with the Hubble Space Telescope Wide Field Camera 3. *Astrophys. J.* **756**, L38 (2012).

2. Marshall, M. A. et al. Limits to rest-frame ultraviolet emission from far-infrared-luminous z-6 quasar hosts. *Astrophys. J.* **900**, 21 (2020).

3. Peng, C. Y. et al. Probing the coevolution of supermassive black holes and galaxies using gravitationally lensed hosts. *Astrophys. J.* **649**, 616–634 (2006).

4. Matsuoka, Y. et al. Subaru high-z exploration of low-luminosity quasars (SHELLQs). I. Discovery of 15 quasars and bright galaxies at 5.7<z<6.9. *Astrophys. J.* **828**, 26 (2016).

5. Matsuoka, Y. et al. Subaru high-z exploration of low-luminosity quasars (SHELLQs). IV. Discovery of 41 quasars and luminous galaxies at 5.7≤z≤6.9. *Astrophys. J. Suppl. Ser.* **237**, 5 (2018).

6. Matsuoka, Y. et al. Subaru high-z exploration of low-luminosity quasars (SHELLQs). XVI. 69 new quasars at 5.8≤z≤7.0. *Astrophys. J. Suppl. Ser.* **259**, 18 (2022).

7. Aihara, H. et al. The Hyper Supreme-Cam SSP survey: overview and survey design. *Publ. Astron. Soc. Jpn* **70**, S4 (2018).

8. Fan, X. et al. A survey of z>5.8 quasars in the Sloan digital sky survey. I. Discovery of three new quasars and the spatial density of luminous quasars at z-6. *Astrophys. J.* **122**, 2833–2849 (2001).

9. Bañados, E. et al. An 800-million-solar-mass black hole in a significantly neutral Universe at a redshift of 7.5. *Nature* **553**, 473–476 (2018).

10. Fan, X., Banados, E. & Simcoe, R. A. Quasars and the intergalactic medium at cosmic dawn. Preprint at https://arxiv.org/abs/2212.06907 (2022).

11. Rieke, M. J., Kelly, D. & Horner, S. Overview of James Webb Space Telescope and NIRCam's role. In Proc. SPIE **5904**, 1–8 (2005).

12. Ding, X. et al. The mass relations between supermassive black holes and their host galaxies at 1<z<2 HST-WFC3. *Astrophys. J.* **888**, 37 (2020).

13. Sérsic, J. L. Boletín de la Asociación Argentina de Astronomía. *La Plata Argentina* **6**, 41–43 (1963).
<figure><img src="image_8.png" /><figcaption>**Fig. 4** | Black hole mass-stellar mass plane at $z \approx 6.4$. The left shows the observed locations of our two $z \approx 6.4$ quasars marked in red. The orange contours (1–3σ levels) show the expected distribution of a mock quasar sample at $z \approx 6.4$ for which the ratio $M_{BH}/M_\odot$ of the intrinsic population does not evolve to $z=0$. The right shows the $M_{BH}-M_\odot$ relation at $z \approx 6.4$ (red solid line) inferred</figcaption></figure>

from our two quasars after accounting for selection effects. The shaded region represents the 68% confidence interval, and the dashed line is the observed relation at $z=0$ (refs. 21,22). Note that we assume a linear relationship between black hole mass and stellar mass in log-log space. The slope is fixed to the value observed at $z=0$.

14. Jakobsen, P. et al. The Near-Infrared Spectrograph (NIRSpec) on the James Webb Space Telescope. I. Overview of the instrument and its capabilities. *Astron. Astrophys.* **661**, A80 (2022).

15. Grazian, A. et al. The galaxy stellar mass function at $3.5 \le z \le 7.5$ in the CANDELS/UDS, GOODS-South, and HUDF fields. *Astron. Astrophys.* **575**, A96 (2015).

16. Vestergaard, M. & Peterson, B. M. Determining central black hole masses in distant active galaxies and quasars. II. Improved optical and UV scaling relationships. *Astrophys. J.* **641**, 689–709 (2006).

17. Li, J. et al. On the connection between supermassive black holes and galaxy growth in the Reionization Epoch. *Astrophys. J. Lett.* **931**, L11 (2022).

18. Kormendy, J. & Ho, L. C. Coevolution (or not) of supermassive black holes and host galaxies. *Annu. Rev. Astron. Astrophys.* **51**, 511–653 (2013).

19. Habouzit, M. et al. Co-evolution of massive black holes and their host galaxies at high redshift: discrepancies from six cosmological simulations and the key role of JWST. *Mon. Not. R. Astron. Soc.* **511**, 3751–3767 (2022).

20. Inayoshi, K. et al. Rapid growth of seed black holes during early bulge formation. *Astrophys. J.* **927**, 237 (2022).

21. Häring, N. & Rix, H.-W. On the black hole mass-bulge mass relation. *Astrophys. J. Lett.* **604**, L89–L92 (2004).

22. Bennert, V. N., Auger, M. W., Treu, T., Woo, J.-H. & Malkan, M. A. A local baseline of the black hole mass scaling relations for active galaxies. I. Methodology and results of pilot study. *Astrophys. J.* **726**, 59 (2011).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

© The Author(s), under exclusive licence to Springer Nature Limited 2023

Peking University, Beijing, China. ⁴Max Planck Institute for Astronomy, Heidelberg, Germany.

⁵Department of Astronomy, School of Science, The University of Tokyo, Bunkyo-ku, Japan.

⁶Research Center for Space and Cosmic Evolution, Ehime University, Matsuyama, Japan.

⁷National Astronomical Observatory of Japan, Osawa, Mitaka, Japan. ⁸Department of Physics, Graduate School of Science, Tokyo Metropolitan University, Hachioji, Japan. ⁹Department of Astrophysical Sciences, Princeton University, Princeton, NJ, USA. ¹⁰Department of Astronomy, University of Illinois at Urbana-Champaign, Urbana, IL, USA. ¹¹Institute of Astrophysics of Paris, CNRS, Sorbonne Université, Paris, France. ¹²Department of Astronomy, Columbia University, New York, NY, USA. ¹³Department of Physics, Columbia University, New York, NY, USA.

¹⁴Physics Department, Technical University of München, Garching bei München, Germany.

¹⁵Max Planck Institute for Astrophysics, Garching bei München, Germany. ¹⁶Subaru Telescope, National Astronomical Observatory of Japan, Hilo, HI, USA. ¹⁷Graduate School of Science and Engineering, Kagoshima University, Kagoshima, Japan. ¹⁸Institute for Computational Science, University of Zurich, Zürich, Switzerland. ¹⁹MIT Kavli Institute for Astrophysics and Space Research, Cambridge, MA, USA. ²⁰Department of Astronomy, The University of Texas at Austin, Austin, TX, USA. ²¹Centre for Astronomy at the University of Heidelberg (ITA), Heidelberg, Germany. ²²Department of Astronomy, School of Science, Graduate University for Advanced Studies (SOKENDAI), Mitaka, Japan. ²³Institute of Sciences of the Cosmos (ICCUB), University of Barcelona (IEEC-UB), Barcelona, Spain. ²⁴ICREA, Barcelona, Spain. ²⁵Research Center for the Early Universe, Graduate School of Science, The University of Tokyo, Bunkyo-ku, Japan. ²⁶Department of Economics, Management and Information Science, Onomichi City University, Onomichi, Japan. ²⁷Institute of Astronomy, Graduate School of Science, The University of Tokyo, Mitaka, Tokyo, Japan. ²⁸W. M. Keck Observatory, Kamuela, HI, USA. ²⁹G. Occhialini' Physics Department, University of Studi di Milano-Bicocca, Milano, Italy.

³⁰Steward Observatory, University of Arizona, Tucson, AZ, USA. ³¹Observatoryl/MCTI, Rua General José Cristino, Rio de Janeiro, Brazil. ³²Leiden Observatory, Leiden University, Leiden, The Netherlands. ³³University of Potsdam, Potsdam, Germany. ³⁴Academia Sinica Institute of Astronomy and Astrophysics, Taipei, Taiwan. ³⁵School of Physics and Astronomy, Tel Aviv University, Tel Aviv, Israel. ³⁶Kapteyn Astronomical Institute, University of Groningen, Groningen, The Netherlands. ³⁷Department of Physics and Astronomy, University of California, Los Angeles, CA, USA. ³⁸Institute for Advanced Research, Nagoya University, Nagoya, Japan.

³⁹Department of Physics, Graduate School of Science, Nagoya University, Nagoya, Japan.
⁴⁰DARK, Niels Bohr Institute, Copenhagen N, Denmark. ⁵e-mail: xuheng.ding@ipmu.jp; onoue@pku.edu.cn

<footnote n="1">¹Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI), The University of Tokyo, Chiba, Japan. ²Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Japan. ³Kavli Institute for Astronomy and Astrophysics,</footnote>
Cosmological model

A standard concordance cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_m = 0.30$ and $\Omega_\Lambda = 0.70$ is adopted, where $H_0$ is the Hubble constant and $\Omega_m$ and $\Omega_\Lambda$ are the matter density and dark energy density, respectively, at the present time, which gives a scale of 5.63 kpc per "at $z = 6.4$. All magnitudes are presented in the AB system. Milky Way dust extinction is negligible at the near-infrared wavelengths of interest and is not corrected for.

Target selection

The two objects presented here are part of our Cycle 1 JWST programme (Observation ID 1967; principal investigator, M. Onoue), which targets 12 of the lowest-luminosity quasars at redshift 6.0 to 6.4 (ref. 23) discovered by the HSC-SSP<sup>7</sup>, an optical wide-field survey using the 8.2-m Subaru Telescope. These quasars were spectroscopically confirmed by ground-based follow-up observations<sup>4–6</sup>. The quasars show broad Ly$\alpha$ emission and strong intergalactic medium absorption bluewards of Ly$\alpha$. This is a deep survey, which is sensitive to quasars a factor of ten less luminous than those discovered by shallower surveys such as the Sloan Digital Sky Survey (Extended Data Fig. 1).

NIRCam observations and data reduction

The data presented in this paper were taken with Module B of the NIRCam instrument, which has a field of view of $2.2 \times 2.2 \text{ arcmin}^2$. Total exposures of 3,100 s in the two filters (F356W and F150W) were obtained simultaneously. A $4 \times 4$ primary and subpixel dithering pattern was used to mitigate cosmic ray hits and bad pixels in the detector and to ensure subpixel resampling during the stacking step. We used the INTRAMODULEBOX and STANDARD dithering patterns for the primary and subpixel dithers, respectively. We used the BRIGHT1 readout mode.

The data were processed using the standard procedures in JWST pipeline v.1.7.2. The precalibrated 'Stage 2' image frames were downloaded from the MAST archive. These images have the pipeline parameter reference files jwst_1009.pmap for J2255+0251 and jwst_1011.pmap for J2236+0032, as registered in the JWST Calibration Reference Data System (https://jwst-crds.stsci.edu). For individual frames, global background light was first subtracted using the Background2D function of Photutils<sup>24</sup>. The archived images clearly have horizontal and vertical stripe noise patterns, known as '1/f noise'. This 1/f noise was subtracted by first masking bright objects, then collapsing the 2D images along each axis of the detectors and estimating the noise amplitudes by measuring sigma-clipped median values. These amplitudes were then subtracted from each row and column. The horizontal stripes were measured for each of the four detector amplifiers separately. These postprocessed stage 2 image frames were then aligned and stacked with inverse-variance weighting using the stage 3 pipeline, keeping the original position angle of the detector for the purpose of building the PSF library. Both F356W and F150W images were resampled with a pixel scale a factor of two smaller than that of the detector, using the drizzling algorithm implemented in the Resample step of the pipeline. The final pixel scales for F356W and F150W are 0.0315 and 0.0153", respectively.

2D image decomposition of quasar and host-galaxy emission using galight

Accurate decomposition of the image of a quasar into a central point source and an extended host galaxy requires a high-quality model for the PSF based on stars detected in the same image as the quasars. Space-based telescopes have a much sharper and more stable PSF than ground-based telescopes, and the Hubble Space Telescope (HST) has been used to measure quasar host galaxies to redshift up to $z \approx 2$ (refs. 12,25–27). However, HST's roughly 90-min orbit means that it is continually passing between Earth's shadow and direct sunlight, causing the telescope to expand and contract ('orbital breathing') and

giving rise to a time-dependent PSF. The difficulty of modelling the PSF has not allowed quasar host galaxies to be detected significantly beyond $z \ge 3$ (refs. 1,2,28,29), where host starlight can be only measured through the SED fitting method<sup>30</sup>.

Here, we follow our tested strategy<sup>12</sup>, including the analysis of JWST images<sup>31</sup>, and build a PSF library by identifying all isolated, unsaturated stars of sufficient signal-to noise ratio in our images. We identified 12/16 PSF stars in filter F150W|F356W for J2236+0032 and 9/5 PSF stars in filter F150W|F356W in the images for J2255+0251. We use our 2D modelling software galight<sup>12</sup> to fit the quasar images with a model of a scaled PSF (the spatially unresolved point-like quasar) and a PSF-convolved 2D Sérsic profile (the host galaxy). galight is a python-based open-source package (https://galight.readthedocs.io/en/latest/) that uses the image modelling capabilities of lenstronomy<sup>32</sup>. It has been adopted for quasar decomposition of images by HST<sup>12</sup>, JWST<sup>31</sup> and the HSC survey<sup>33</sup>. We adopt uniform priors for the effective radius $R_{\text{eff}} \in [0.03", 2.0"]$ and the Sérsic index ($n) \in [0.3, 9]$ of the host to avoid unphysical parameter inference. For J2255+0251, there is one nearby object that we fit simultaneously with a Sérsic profile, allowing us to correct for its effect on the host-galaxy properties.

We obtain a weighted inference for the decomposition result as briefly described here. After subtracting the remaining local background, we use each PSF in our library in turn to fit the image. The performance of each PSF is determined by its best-fit $\chi^2$ value. We select groups of two, three and five PSFs from the library that have the top level $\chi^2$ performance and then average them using psfr (S. Birrer et al., manuscript in preparation). To optimize our modelling of the unresolved quasar emission, we consider both best-fit models using individual stars and average models based on the combined PSF stars described above. Thus, we add the three averaged PSFs as new members to the PSF library. We take the results from the five top-performing PSF models in the updated library. We determine our final result parameters by weighting their $\chi^2$ values, as defined by ref. 12, equations (3)–(6). The uncertainties we quote are based on the dispersion in the host properties from these five different PSF models. Errors based on variation in the PSF models are often significantly larger than the formal fitting errors, as has been recently demonstrated using simulations and observations of quasars and their host galaxies<sup>34</sup>. The final inferred host filter magnitude, size, Sérsic index and other fit parameters are presented in Extended Data Table 1. The images in Fig. 2 are based on the PSF with the best performance. The sizes of our host galaxies are defined using the Sérsic effective radius $R_{\text{eff}}$ along the semimajor axis measured by galight.

The quasar host of J2236+0032 seems to be a compact galaxy. As a result, we find that the Sérsic $n$ is poorly constrained, and the resulting host residual is point-like with an inferred $R_{\text{eff}}$ of 0.03" (the lower limit of the prior), indicating that the central PSF has not been fully removed. This model implies a stellar mass with an implausibly high value above $10^{11.8} M_{\odot}$ given the SED fitting described below. Thus, we refit this object, fixing $n = 1$ in both F356W and F150W. This gives a much-improved model. The host is elongated in F356W, suggesting that it is a disc-like edge-on galaxy. We also find a consistent position angle between F356W and F150W, as shown in Fig. 2. Thus, for J2236+0032, we adopt the results with Sérsic $n$ fixed to 1. We further directly fit the 'data - quasar' residual image (that is, host) with the Sérsic model and allow the Sérsic index to vary. When we do so, the inferred $n$ is 1.55 with values of $R_{\text{eff}}$ and the host magnitude consistent with those with the $n = 1$ fit.

As J2255+0251 has no clear host detection in the F150W band, we remodel the data and fix the host-galaxy parameters (host Sérsic index, $R_{\text{eff}}$, ellipticity and position angle) to the values inferred from the F356W band, while allowing the central position and amplitude to vary. We still see no evidence for a host residual, but the model gives a formal host magnitude of $26.3 \pm 0.2$ mag. We quote these inferences in Extended Data Table 1 and adopt a lower limit to the magnitude of 26.3 mag.
## Confirmation of our 2D decomposition using galfit

In the fitting routine described above, the results are sensitive to the setting of certain parameters. For instance, a subpixel shift of the PSF is typically needed to optimally align with the position of the quasar. This interpolation within a pixel can result in modest variations in the allocation of flux between the quasar and its host galaxy. In galight, the parameter point_source_supersampling_factor controls this interpolation. For guidance here on the optimal setting, we compared our results with the decomposition from galfit³⁵ while implementing different supersampling factor values, that is, 1 (no interpolation) and 2. The choice of supersampling factor had a greater impact on the F150W decomposition than that on the F356W image. We found that the two codes are in close agreement when using a supersampling_factor of 2 as shown in Extended Data Fig. 2 and reported in Extended Data Table 1. As a result, to account for the uncertainty associated with this setting, we report the averaged results obtained using both supersampling_factor values of 1 and 2, with uncertainties covered by these two settings.

## Tests of the robustness of the host-galaxy properties

We perform several sanity checks for our two quasars to confirm that the detections of their host galaxies are real. First, our final reduced data are co-added using 16 dithers. To check whether the apparent host is caused by a random ghost from any particular dither frame, we reanalyse the data using the first eight and the second eight dithers separately. We are able to detect the host from both halves of the data. To ensure that the apparent host-galaxy flux is not dominated by a mismatched PSF core, we mask the quasar centre (using a 0.12 and 0.06" radius aperture for F356W and F150W, respectively) and redo the fit. Again, the host is well-detected.

To rule out the possibility that the residual emission (host galaxy) is an artefact of the chosen PSF stars, we selected the five top-performing PSF in the library, and identified 20 pairs from these. We then fit one PSF with the other without assuming any host component for all pairs. We performed this test for stars in the fields of J2236+0032 and J2255+0251 in both bands; residual images show no extended features in any case. This is illustrated for J2255+0251 in the F356W filter in the top panels of Extended Data Fig. 3. We also fit galaxy + host models to pure PSF images, in which case the inferred host magnitudes were roughly two magnitudes fainter than the hosts detected around the quasars. By contrast, we clearly detect the quasar host when we use any of the five library PSFs and, using these stellar residual maps to measure the standard deviation, we show the maps of the host signal-to-noise ratio for all images of our two quasars in the bottom panels of Extended Data Fig. 3.

We carried out three sets of simulations to further validate our host detections, including J2255+0251 in the F356W filter and J2236+0032 in the F356W and F150W filters. In one realization of the simulation, we randomly select a PSF from the five top-ranked PSFs in the library, and add a Sérsic profile using the best-fit parameters from our model inference. We sample the resulting mock quasar + host image at the observed resolution and insert it into a random empty position in the JWST NIRCam data. We then model the images with the other four PSF stars, and select the result with the best performance. We repeated this simulation 100 times to quantify the bias and standard deviation of the inferred values around the true values. We found that the key Sérsic parameters can be obtained accurately, with biases well below $1\sigma$ in all quantities. For J2255+0251, the bias (obtained – truth) and scatter of the inferred host magnitude are $-0.02 \pm 0.11$. For J2236+0032, the magnitude bias and scatter are $0.01 \pm 0.14$ and $-0.13 \pm 0.16$ for F356W and F150W, respectively. We present the distributions of the inferred host magnitudes in the simulations of J2255+0251 (F356W) in Extended Data Fig. 4. We note that these simulation tests were conducted with point_source_supersampling_factor set to 1. The similar scatter distribution can be found when setting this parameter to 2.

Finally, we tested the fidelity of the host magnitude inference using a joint fit across the two bands. We fix the host-galaxy parameters for J2236+0032 (where the host is detected in F150W), setting the Sérsic parameters to those inferred from F356W; we find a very similar host magnitude to the value when all parameters are allowed to vary: the changes in the inferred host magnitude are within 0.1 mag.

## Stellar mass of the host galaxy

We use SED fitting to estimate the stellar masses of our two host galaxies at $z \simeq 6.4$ using the photometry and its uncertainty based on NIRCam imaging with two filters (F356W and F150W) that straddle the rest-frame 4,000 Å break. The host galaxy of J2255+0251 is not detected in the F150W filter, thus we take 26.3 mag as a lower limit on the magnitude. We adopt a Chabrier initial mass function in our analysis to allow direct comparison with the local black hole–host-galaxy mass relations (for example, refs. 21,22). Contributions of nebular emission lines from the host galaxies are also allowed with a uniform log ionization parameter $\log U$ over the range [-3, -1]. Three key parameters define the shape of the SED template: age, metallicity and dust attenuation ($A_\nu$). We adopt a uniform prior on the age over the range [0.01, 0.84] (in Gyr), extending to the age of the Universe at $z = 6.4$. We adopt a uniform prior in metallicity $\log Z/Z_\odot$ over the range [-1, -0.3] following refs. 36,37. We assume an $A_\nu$ range of [0, 5.0] mag for both targets. For J2255+0251, the F150W band has an upper limit on the flux, thus the constraint on $A_\nu$ is weaker; hence, we assume $A_\nu$ follows a log-normal distribution with a median value at 0.85 as supported by recent JWST observations of high-redshift galaxies³⁸. A measurement of the Balmer decrement (that is, the ratio of narrow Hα and Hβ in the spectrum) implies an attenuation $A_\nu = 0.67 \pm 0.13$ mag for the narrow-line region (M. Onoue et al., manuscript in preparation), consistent with the median extinction value of our prior. We also note that for J2236+0032, the contribution of the host in the spectrum (M. Onoue et al., manuscript in preparation) has a magnitude in close agreement with the NIRCam photometry, supporting our stellar mass estimate.

We use gsf³⁹ to perform the SED model fitting. This software generates a set of templates with a range of ages and metallicities according to the prior to fit a composite stellar population-like star formation history. A random parameter sampling is performed through Markov chain Monte Carlo to infer the probability distribution of the SED parameters. For J2236+0032 and J2255+0251, the inferred values of $\log M_\odot$ are $11.12_{-0.27}^{+0.40}$ and $10.53_{-0.37}^{+0.51}$ (in units of $M_\odot$), respectively, as shown in Extended Data Table 1. In Extended Data Fig. 5, we present the best-fit SEDs of the two quasar hosts, in which the Markov chain Monte Carlo inference indicates that the stellar mass is mildly sensitive to age and metallicity but is more sensitive to $A_\nu$. We also adopt the same configuration for the independent code Bagpipes⁴⁰ and CIGALE⁴¹ to perform SED fits; we find values of $M_\odot$ that change by 0.1 dex or less. The robust measurement of stellar mass is due to having photometry on both sides of the 4,000 Å break, accurate spectroscopic redshifts and a firm upper limit of the stellar age 0.84 Gyr, given the age of the Universe at $z = 6.4$. Even so, we note that extra systematic uncertainties may be unaccounted for thus leading to greater uncertainty in the mass estimates than we have quoted.

## NIRSpec observations and data reduction

The NIRSpec fixed-slit data for the two quasars were obtained on 30 October 2022 for J2236+0032 and 28 October 2022 for J2255+0251. The targets were aligned onto the $0.2 \times 3.2''$ S200A2 slit and observed with the medium-resolution ($R \simeq 1000$) G395M grating, which covers $2.87-5.27$ µm in the observed frame. Total exposures of 1,970 s were obtained for each target.

The spectroscopic data were processed with the JWST pipeline v1.8.5. Raw (uncal) data were downloaded from the MAST archive, which has the parameter reference files jwst_1027.pmap. The stage 1 calibration was performed with the default parameters but with the function to
detect large cosmic ray hits on the detector turned on during the Jump step. The 1/f/noise along the vertical direction was subtracted from the stage 1 products using the same method as that for NIRCam images. Pixels of the science frames flagged as bad or saturated were masked on the basis of the associated data quality frames. During the stage 2 data processing, including flat-fielding and flux calibration, we use in-flight measurements of spectrophotometric standard stars (jwst\_nirspec\_fflat\_0065.fits) to convert the detector signals into physical units. The signal loss due to the optical system and the finite width of the S200A2 slit were corrected in the Pathloss step by applying built-in correction factors for a point source as a function of wavelength. The 2D spectra at each dither position were then stacked using the stage 3 pipeline. The pixel resampling scale was kept to the original pixel scale of the detector. A 6-pixel-wide box-car aperture (0.6" wide) was used to extract the one-dimensional spectrum from the stacked spectrum.

## Spectrum fitting

We use a public code QSOFitMORE⁴² (v1.2.0) (https://doi.org/10.5281/ zenodo.4893646) with custom modification to fit the NIRSpec rest-optical spectra of the two quasars. The continuum emission is fitted with two components: power-law emission from the accretion disc and a pseudo-continuum from the forest of singly ionized iron emission. QSOFitMORE uses an empirical iron template of ref. 43. The rest-frame 5,100 Å monochromatic luminosity ($L_{5,100}$) is derived from the best-fit power-law continuum model with a power-law index of $\alpha_{\lambda}$ (that is, flux density $F_{\lambda} \propto \lambda^{\alpha_{\lambda}}$).

Emission lines are modelled with a combination of one or more Gaussian profiles. J2236+0032 needs one broad component for $H\beta$, whereas the adjacent forbidden emission of [O III] $\lambda\lambda$4960, 5008 requires one narrow component and another blueshifted broad component. The narrow [O III] doublet is forced to have the same line width without velocity offsets between the two. The [O III] redshifts of the two targets are derived from the peaks of the narrow [O III] components. The [O III] redshift of J2236+0032 agrees well with the location of the absorption lines we detect in Hγ, Hδ and Hε (Fig. 3). The $H\beta$ + [O III] line profile as well as the continuum properties of J2236+0032 are reported in Extended Data Table 2. A full description of the fitting procedure and other emission line properties will be presented in M. Onoue et al. (manuscript in preparation).

## Black hole mass estimates

The central black hole masses $M_{BH}$ of the two quasars are measured by the single-epoch method using broad $H\beta$ emission lines. Assuming that the broad $H\beta$ emission traces the kinematics of the gas orbiting in the vicinity of the SMBHs (the so-called broad-line region), one can use the line width of $H\beta$ as a proxy of the rotation velocity of the broad-line-region gas. There is also an empirical correlation between the rest-frame 5,100 Å continuum luminosity ($L_{5,100}$) and the radius from the SMBH to the broad-line-region gas, which is observed by the reverberation mapping technique⁴⁴. In this paper, we use the recipe provided by ref. 16 to estimate the black hole mass:

$$M_{BH} = 10^{6.91} \left( \frac{L_{5100}}{10^{44} \text{ erg s}^{-1}} \right)^{0.50} \left( \frac{\text{FWHM}_{H\beta}}{10^3 \text{ km s}^{-1}} \right)^2 M_{\odot} \quad (1)$$

where FWHM means the full-width at half-maximum. Based on the spectrum fitting results described above, we derive virial black hole masses of $M_{BH} = (1.54 \pm 0.27) \times 10^9 M_{\odot}$ for J2236+0032 and $M_{BH} = (2.02 \pm 0.17) \times 10^8 M_{\odot}$ for J2255+0251. Note that these virial SMBH masses have intrinsic uncertainties of 0.4 dex (ref. 16) that are not included in the error budget above. We report the derived BH properties in Extended Data Table 2.

We estimate the host-galaxy contribution to the rest-optical continuum of the NIRSpec spectra, using the decomposed host-galaxy images in the NIRCam's F356W filter. Referring to the relative position

angle of the NIRCam and NIRSpec observations, we find that the stellar emission falling inside the $0.2 \times 0.6''$ extraction aperture is 21.8% (J2236+0032) and 4.8% (J2255+0251) with respect to the total emission. We interpret this as the fraction of the host starlight in the observed continuum emission in the NIRSpec G395M spectrum for each source, which allows us to correct the 5,100 Å continuum luminosity to that from the quasar emission alone. Accordingly, the virial BH masses become $M_{BH} = (1.36 \pm 0.15) \times 10^9 M_{\odot}$ for J2236+0032 and $M_{BH} = (1.97 \pm 0.17) \times 10^8 M_{\odot}$ for J2255+0251, respectively. These corrected BH properties are reported in Extended Data Table 2, and are adopted in Fig. 4. This correction is more significant for J2236+0032, in which stellar absorption lines are detected (Fig. 3).

## Emission from $H\beta$ + [O III] in the hosts

Our SED model allows for the inclusion of the $H\beta$ + [O III] line emission that falls within the F365W filter (Extended Data Fig. 5). However, if those lines were of high equivalent width in the hosts, they could affect the broad-band photometry used in the stellar mass inference, especially if the gas is ionized by the quasar⁴⁵,⁴⁶. To ensure that our NIRCam detection of extended host emission is not dominated by these lines, we directly compare the total flux of the $H\beta$ + [O III] narrow lines seen in our spectra with the contribution from the host in the NIRCam image within the aperture of the spectral slit. For J2236+0032, narrow $H\beta$ emission is not detected and the [O III] emission in NIRSpec contributes only roughly 3% of the host-galaxy light.

For J2255+0251, the narrow $H\beta$ + [O III] emission is comparable to the host-galaxy brightness measured in the F356W filter. Much of this emission may be coming from the narrow-line region of the quasar itself. To quantify the impact on the detection of the stellar continuum, we use the 2D NIRSpec spectroscopic data and fit a Gaussian to the unresolved spatial profile at each wavelength, allowing the centroid to vary linearly with wavelength, setting the width to be the same as the spatial PSF as measured from a publicly available reference star in the same spectroscopic configuration as our observations and letting the amplitude be a free parameter. The best-fit Gaussian profile at each wavelength is then subtracted from the 2D spectrogram to show any spatially extended [O III] emission in the spectrum. A background region of equivalent size is subtracted from the extended emission region. We calculate the flux density of the extended emission in the F356W band and compare it to the flux density of the NIRCam host galaxy within the aperture of the NIRSpec slit on the sky. For J2255+0251, the [O III] extended emission has a magnitude of 27.6, which is 11% of the flux of the observed host galaxy within the slit aperture. Note that the measurement of the extended emission is based on the observational fact that the apparent minor axis of J2255+0251 is wider than the scale of the PSF in F356W (Fig. 2). We thus conclude that most of the narrow $H\beta$ + [O III] emission in J2255+0251 comes from the unresolved quasar, which has been removed by our quasar decomposition. For J2236+0032, we use the same approach and find that the extended [O III] has a magnitude of 27.5 mag. The flux ratio between the extended [O III] and the host galaxy within the slit is roughly 3%. No extended $H\beta$ emission was detected in either target.

## Construction of mock quasar and host-galaxy sample at z = 6.4

We use the method described in ref. 17 to generate mock samples to assess the expected location of our two quasars in the black hole–stellar mass plane. We start with a large galaxy sample representative of the population at z ≈ 6 based on the stellar mass function¹⁵. We then assume that the ratio of black hole to host stellar mass and its dispersion are the same as the local values (that is, no evolution). As an update on this procedure, we now use the black hole Eddington rate distribution of ref. 47 to assign quasar luminosities. Observational uncertainty on M. and intrinsic uncertainty on $M_{BH}$ are incorporated with a selection on magnitude (−23.9 < $M_{1450}$ < −23.8) is applied. The orange contours in Fig. 4 show the resulting bivariate probability distribution in the stellar
mass-black hole mass plane. We further generate samples for different evolution rates ($M_{BH}/M. \propto (1+z)^\gamma$) in steps of 0.1 in $\gamma$ and under the assumption that the dispersion is the same as the local value (refs. 17,48 demonstrates that the evolution rate and dispersion are degenerate). The probability distribution of $\gamma$ is determined given our measurements of both $M_{BH}$ and $M$ of our two quasars (Extended Data Fig. 6); $\gamma$ lies in the range of -0.40 to +0.84 with 68% confidence.

## Data availability

The JWST data are available in the MAST and can be accessed through https://archive.stsci.edu/doi/resolve/resolve.html?doi=10.17909/2z59-3v32 (JWST GO 1967) after a 12-month exclusive access period. Reduced quasar NIRCam image cut-outs are also available from the corresponding author upon request or at https://github.com/dartoon/publication/tree/main/SHELLQs_J2236_J2255_data.

## Code availability

The JWST data were processed with the JWST calibration pipeline (https://jwst-pipeline.readthedocs.io). Public tools were used for data analysis: galight¹², galfit³⁵ and PyQSOFit⁴².

23. Matsuoka, Y. et al. Subaru high-z exploration of low-luminosity quasars (SHELLQs). V. Quasar luminosity function and contribution to cosmic reionization at z=6. *Astrophys. J.* **869**, 150 (2018).

24. Bradley, L. et al. astropy/photutils 1.4.0. Zenodo https://zenodo.org/record/6385735 (2022).

25. Jahnke, K. et al. Ultraviolet light from young stars in GEMS Quasar host galaxies at 1.8<z<2.75. *Astrophys. J.* **614**, 568–585 (2004).

26. Mechtley, M. et al. Do the most massive black holes at z=2 grow via major mergers? *Astrophys. J.* **830**, 156 (2016).

27. Marian, V. et al. Major mergers are not the dominant trigger for high-accretion AGNs at z-2. *Astrophys. J.* **882**, 141 (2019).

28. Glikman, E. et al. Major mergers host the most-luminous red quasars at z-2: a Hubble Space Telescope WFC3/IR study. *Astrophys. J.* **806**, 218 (2015).

29. Chen, Y.-C. et al. A close quasar pair in a disk-disk galaxy merger at z=2.17. *Nature* **616**, 45–49 (2023).

30. Dong, X. Y. & Wu, X.-B. Herschel observed Stripe 82 quasars and their host galaxies: connections between AGN activity and host galaxy star formation. *Astrophys. J.* **824**, 70 (2016).

31. Ding, X., Silverman, J. D. & Onoue, M. Opening the era of quasar-host studies at high redshift with JWST. *Astrophys. J. Lett.* **939**, L28 (2022).

32. Birrer, S. et al. Lenstronomy II: a gravitational lensing software ecosystem. *J. Open Source Softw.* **6**, 3283 (2021).

33. Li, J. et al. The sizes of quasar host galaxies in the Hyper Suprime-Cam Subaru Strategic Program. *Astrophys. J.* **918**, 22 (2021).

34. Zhuang, M.-Y. & Shen, Y. Characterization of JWST NIRCam PSFs and implications for AGN+host image decomposition. Preprint at https://arxiv.org/abs/2304.13776 (2023).

35. Peng, C. Y., Ho, L. C., Impey, C. D. & Rix, H.-W. Detailed structural decomposition of galaxy images. *Astron. J.* **124**, 266–293 (2002).

36. Bouwens, R. J. et al. UV-continuum slopes of >4000 z-4-8 galaxies from the HUDF/XDF, HUDF09, ERS, CANDELS-South, and CANDELS-North fields. *Astrophys. J.* **793**, 115 (2014).

37. Ucci, G. et al. Astraeus - II. Quantifying the impact of cosmic variance during the Epoch of Reionization. *Mon. Not. R. Astron. Soc.* **506**, 202–214 (2021).

38. Shapley, A. E., Sanders, R. L., Reddy, N. A., Topping, M. W. & Brammer, G. B. JWST/NIRSpec Balmer-line measurements of star formation and dust attenuation at z-3-6. Preprint at https://arxiv.org/abs/2301.03241 (2023).

39. Morishita, T. et al. Massive dead galaxies at z-2 with HST Grism spectroscopy. I. Star formation histories and metallicity enrichment. *Astrophys. J.* **877**, 141 (2019).

40. Carnall, A. C., McLure, R. J., Dunlop, J. S. & Davé, R. Inferring the star formation histories of massive quiescent galaxies with BAGPIPES: evidence for multiple quenching mechanisms. *Mon. Not. R. Astron. Soc.* **480**, 4379–4401 (2018).

41. Boquien, M. et al. CIGALE: a python code investigating galaxy emission. *Astron. Astrophys.* **622**, A103 (2019).

42. Fu, Y. QSOFITMORE (version 1.2.0): a python package for fitting UV-optical spectra of quasars. Zenodo https://zenodo.org/record/5810042 (2021).

43. Boroson, T. A. & Green, R. F. The emission-line properties of low-redshift quasi-stellar objects. *Astrophys. J. Suppl. Ser.* **80**, 109 (1992).

44. Kaspi, S. et al. Reverberation measurements for 17 quasars and the size-mass-luminosity relations in active galactic nuclei. *Astrophys. J.* **533**, 631–649 (2000).

45. Kashino, D. et al. EIGER I. A large sample of [O iii]-emitting galaxies at 5.3<z<6.9 and direct evidence for local reionization by galaxies. *Astrophys. J.* **950**, 66 (2023).

46. Mathee, J. et al. EIGER II. First spectroscopic characterisation of the young stars and ionised gas associated with strong Hβ and [O iii] line-emission in galaxies at z=5-7 with JWST. *Astrophys. J.* **950**, 67 (2023).

47. Wu, J. et al. Demographics of z-6 quasars in the black hole mass-luminosity plane. *Mon. Not. R. Astron. Soc.* **517**, 2659–2676 (2022).

48. Li, J. et al. Synchronized coevolution between supermassive black holes and galaxies over the last seven billion years as revealed by Hyper Suprime-Cam. *Astrophys. J.* **922**, 142 (2021).

**Acknowledgements** We thank T. Morishita and L. Yang for help with the gsf and Bagpipes packages. We thank Y. Fu for his help on the use of QSOFitMORE. We thank X. Fan, L. Ho and V. Bennert for useful discussions. This work is based on observations made with the NASA/ESA/CSA JWST. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract no. NAS 5-03127 for JWST. These observations are associated with programme no. 1967. Support for programme no. 1967 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract no. NAS 5-03127. This paper is based on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by the Subaru Telescope and Astronomy Data Center at NAOJ. Data analysis was in part carried out with the cooperation of the Center for Computational Astrophysics (CfCA), NAOJ. We are honoured and grateful for the opportunity of observing the Universe from Maunakea, which has cultural, historical and natural significance in Hawaii. This work was supported by World Premier International Research Center Initiative (WPI), MEXT, Japan. X.D., J.S., Y.M., T.I., K.K. and H.U. are supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI grant nos. JP22K14071, JP18H01251, JP22H01262, JP17H04830, JP21H04494, JP20K14531, JP17H06130 and JP20H01953. F.W., S.E.I.B., B.T., J.T.S. and M.O. acknowledge support from the European Research Council (ERC) grant nos. 740246, 885301 and 950533. M.O. and K.I. acknowledge support from the National Natural Science Foundation of China grant nos. 12150410307 and 12073003. Y.M. acknowledges support from the Mitsubishi Foundation grant no. 30140. S.F. acknowledges support from NASA through the NASA Hubble Fellowship grant no. HST-HF2-51505.001-A awarded by the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract no. NAS 5-26555. Z.H. acknowledges support from the National Science Foundation grant no. AST-2006176. K.I. acknowledges support by the Spanish MCIN under grant no. PID2019-105510GB-C33/AEI/10.13039/501100011033. A.L. acknowledges funding from MIUR under the grant no. PRIN 2017-MB8AEZ. B.T. acknowledges support from the Israel Science Foundation (grant no. 1849/19) and from the ERC under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 950533). M.T. acknowledges support from the NWO grant no. 0.16. VIDI.189.162 ('ODIN'). M.V. acknowledges support from the Independent Research Fund Denmark by grant no. DFF 8021-00130.

**Author contributions** X.D. led the analysis of the quasar decomposition, simulation test, SED fitting and preparation of the manuscript. M.O. reduced the NIRCam and NIRSpec data. He also led the spectroscopic data analysis and wrote relevant sections. These first two authors equally contributed to this work. J.D.S. contributed to the NIRCam data analysis and preparation of the manuscript. C.L.P. performed the NIRSpec 2D spectroscopic analysis and evaluated the strength of the extended [O iii] emission. J.L. constructed the mock quasar sample at z=6.4, which was used to assess the redshift evolution of the black hole-stellar mass relation. M.O., X.D., J.D.S., Y.M., T.I., M.A.S. and K.J. led the project design and management, also developing the main interpretation of the results. All coauthors contributed to the discussion of the presented results and the preparation of the manuscript.

**Competing interests** The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41586-023-06345-5.

**Correspondence and requests for materials** should be addressed to Xuheng Ding or Masafusa Onoue.

**Peer review information** Nature thanks the anonymous reviewers for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permissions information** is available at http://www.nature.com/reprints.
<figure><img src="image_9.png" /></figure>

**Extended Data Fig. 1** | Distribution of known quasars at redshift $z > 5.8$. The y-axis indicates the absolute magnitudes at rest-frame 1450 Å. The two target quasars in this work are shown in red (J2236+0032 at $z = 6.40$, and J2255+0251

at $z = 6.34$), while other low-luminosity quasars from the HSC-SSP are shown in blue. The J/WST 12 Cycle 1 targets in GO #1967 are highlighted with open circles. Other known quasars are shown in black.
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

Extended Data Fig. 2 | Quasar-subtracted images of J2236+0032 based on 2D decomposition by two independent codes. The inferred host magnitude from each code is indicated. For this comparison, the point_source_

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

supersampling_factor is set to 2 for galight. The same PSF is used for this comparison. Coordinates are in pixel units; the pixel scales for F356W and F150W are 0.0315 and 0.0153, respectively.
<figure><img src="image_18.png" /><figcaption>**Extended Data Fig. 3** | Testing the effect of PSF star mismatch on the inferred presence of a host galaxy. Top: The residuals using the top five PSF star pairs for J2255+0251 in F356W (10 of 20 are shown) – we see no evidence for extended residual emission. Using these 20 PSF residual maps, the PSF</figcaption></figure>

mismatch *standard deviation* is calculated for each pixel. *Bottom*: Maps of host galaxy signal (i.e., data – point-source model) divided by the standard deviation (as indicated by the color scale) for both targets in two filters. The host is clearly detected in all cases except for J2255+0251 in F150W.
<figure><img src="image_19.png" /><figcaption>**Extended Data Fig. 4** | Distribution of the inferred host magnitude from simulations with and without a host present. When the host of J2255+0251 is added into the mock image in the simulations, the derived host magnitude distribution based on 100 realizations (dark orange) is consistent with the</figcaption></figure>

input true value, with an r.m.s. scatter of 0.12 mag. When the host light is not included in the mock image and we still fit the host light in our model, we infer host magnitudes with a broad distribution (tan histogram), centered two magnitudes fainter.
<figure><img src="image_20.png" /></figure>

**Extended Data Fig. 5** | The SED inference using the host galaxy two-band photometry based on gsf. **Top:** The red data points with errors indicate the inferred host fluxes. The gray region represents the 1σ range of variation in SED templates obtained through MCMC, while the black line represents the median

SED template. The blue diamonds indicate the predictions using this median template. **Bottom:** The distribution of the host galaxy SED properties using MCMC inference.
<figure><img src="image_21.png" /><figcaption>**Extended Data Fig. 6** | Probability distribution of the evolution parameter based on observations of our two quasars. The evolution parameter *y* is inferred as M<sub>BH</sub>/M<sub>☉</sub> ∝ (1+z)<sup>y</sup>. The 68% confidence interval in *y* lies between −0.40 and +0.84. The black dashed line indicates the 50 percentile position (+0.21).</figcaption></figure>
Extended Data Table 2 | Spectroscopic properties and virial BH masses

<table><thead><tr><th></th><th>J2236+0032</th><th>(corrected)</th><th>J2255+0251</th><th>(corrected)</th></tr></thead><tbody><tr><td>[O III] redshift</td><td>6.4039 ± 0.0009</td><td></td><td>6.333 ± 0.001</td><td></td></tr><tr><td>F<sub>Hβ,broad</sub></td><td>2.13 ± 0.10</td><td></td><td>7.01 ± 0.09</td><td></td></tr><tr><td>F<sub>Hβ,narrow</sub></td><td></td><td></td><td>0.36 ± 0.07</td><td></td></tr><tr><td>F<sub>[O III]λ4960,core</sub></td><td>0.07 ± 0.10</td><td></td><td>0.26 ± 0.04</td><td></td></tr><tr><td>F<sub>[O III]λ4960,wing</sub></td><td>0.49 ± 0.14</td><td></td><td>0.31 ± 0.02</td><td></td></tr><tr><td>F<sub>[O III]λ5008,core</sub></td><td>0.34 ± 0.13</td><td></td><td>0.66 ± 0.09</td><td></td></tr><tr><td>F<sub>[O III]λ5008,wing</sub></td><td>0.96 ± 0.18</td><td></td><td>1.36 ± 0.08</td><td></td></tr><tr><td>α<sub>λ</sub></td><td>-1.752 ± 0.003</td><td></td><td>-1.75 ± 0.02</td><td></td></tr><tr><td>L<sub>5100</sub> (10<sup>45</sup> erg s<sup>-1</sup>)</td><td>2.29 ± 0.02</td><td>1.79 ± 0.01</td><td>1.531 ± 0.008</td><td>1.458 ± 0.008</td></tr><tr><td>Hβ FWHM<sub>broad</sub> (km s<sup>-1</sup>)</td><td>6290 ± 560</td><td></td><td>2520 ± 110</td><td></td></tr><tr><td>M<sub>BH</sub> (10<sup>8</sup> M<sub>☉</sub>)</td><td>15.4 ± 2.7</td><td>13.6 ± 1.5</td><td>2.02 ± 0.17</td><td>1.97 ± 0.17</td></tr></tbody></table>

The corrected values are based on the quasar continuum emission after subtracting the host galaxy light falling onto the NIRSpec fixed slit. Flux values are in units of 10<sup>-17</sup> erg cm<sup>-2</sup> s<sup>-1</sup>.