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

# The $M_{BH}-M_*$ Relation up to $z \sim 2$ through Decomposition of COSMOS-Web NIRCam Images

Takumi S. Tanaka¹,²,³, John D. Silverman¹,²,³,⁴, Xuheng Ding¹,²,³,⁵, Knud Jahnke⁶, Benny Trakhtenbrot⁷, Erini Lambrides³,²,⁸, Masafusa Onoue²,⁹, Irham Taufik Andika¹⁰,¹¹, Angela Bongiorno¹², Andreas L. Faisst¹³, Steven Gillman¹⁴,¹⁵, Christopher C. Hayward¹⁶, Michaela Hirschmann¹⁷, Anton Koekemoer¹⁸, Vasily Kokorev¹⁹, Zhaoxuan Liu¹,²,³, Georgios E. Magdis¹⁴,²⁰,²¹, Alvio Renzini²², Caitlin Casey¹⁴,²³, Nicole E. Drakos²⁴, Maximilien Franco²³, Ghassem Gozaliasl²⁵,²⁶, Jeyhan Kartaltepe²⁷, Daizhong Liu²⁸, Henry Joy McCracken²⁹, Jason Rhodes³⁰, Brant Robertson³¹, and Sune Toft¹⁴,²¹

¹Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan; takumi.tanaka@ipmu.jp

²Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

³Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

⁴Center for Astrophysical Sciences, Department of Physics & Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA

⁵School of Physics and Technology, Wuhan University, Wuhan 430072, People's Republic of China

⁶Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

⁷School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel

⁸NASA-Goddard Space Flight Center, Code 662, Greenbelt, MD, 20771, USA

⁹Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, People's Republic of China

¹⁰Technical University of Munich, TUM School of Natural Sciences, Department of Physics, James-Franck-Str. 1, D-85748 Garching, Germany

¹¹Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany

¹²INAF-Osservatorio Astronomico di Roma, Via Frascati 33, I-00078, Monte Porzio Catone, Italy

¹³Caltech/IPAC, 1200 E. California Boulevard, Pasadena, CA 91125, USA

¹⁴Cosmic Dawn Center (DAWN), Denmark

¹⁵DTU Space, Technical University of Denmark, Elektrovej, Building 328, 2800, Kgs. Lyngby, Denmark

¹⁶Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA

¹⁷Institute for Physics, Laboratory for Galaxy Evolution and Spectral Modelling, EPFL, Observatoire de Sauverny, Chemin Pegasi 51, 1290 Versoix, Switzerland

¹⁸Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

¹⁹Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, NL-9700AV Groningen, The Netherlands

²⁰DTU-Space, Technical University of Denmark, Elektrovej 327, DK2800 Kgs. Lyngby, Denmark

²¹Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N, Denmark

²²Osservatorio Astronomico di Padova, Vicole dell'Osservatorio 5, Padova, I-35122, Italy

²³The University of Texas at Austin, 2515 Speedway Boulevard Stop C1400, Austin, TX 78712, USA

²⁴Department of Physics and Astronomy, University of Hawaii, Hilo, 200 W Kawili Street, Hilo, HI 96720, USA

²⁵Department of Computer Science, Aalto University, PO Box 15400, Espoo, FI-000 076, Finland

²⁶Department of Physics, Faculty of Science, University of Helsinki, 00014-Helsinki, Finland

²⁷Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA

²⁸Max-Planck-Institut für extraterrestrische Physik, Gießenbachstraße 1, 85748 Garching b. München, Germany

²⁹Institut d'Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98 bis boulevard Arago, F-75014 Paris, France

³⁰Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91001, USA

³¹Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA

Received 2024 January 25; revised 2024 December 3; accepted 2024 December 3; published 2025 January 28

## Abstract

Our knowledge of relations between supermassive black holes and their host galaxies at $z \gtrsim 1$ is still limited, even though being actively sought out to $z \sim 6$. Here, we use the high resolution and sensitivity of JWST to measure the host galaxy properties for 107 X-ray-selected type-I active galactic nuclei (AGNs) at $0.68 < z < 2.5$ with rest-frame optical/near-infrared imaging from COSMOS-Web and PRIMER. Black hole masses ($\log(M_{BH}/M_{\odot}) \sim 6.9-9.6$) are available from previous spectroscopic campaigns. We extract the host galaxy components from four NIRCam broadband images and the Hubble Space Telescope/Advanced Camera for Surveys F814W image by applying a 2D image decomposition technique. We detect the host galaxy for $\sim 90\%$ of the sample after subtracting the unresolved AGN emission. With host photometry free of AGN emission, we determine the stellar mass of the host galaxies to be $\log(M_*/M_{\odot}) \sim 9.5-11.6$ through spectral energy distribution fitting and measure the evolution of the mass relation between SMBHs and their host galaxies. Considering selection biases and measurement uncertainties, we find that the $M_{BH}/M_*$ ratio evolves as $(1+z)^{0.48_{-0.62}^{+0.31}}$ thus remains essentially constant or exhibits mild evolution up to $z \sim 2.5$. We also see an amount of scatter ($\sigma_{\mu} = 0.30_{-0.13}^{+0.14}$), similar to the local relation and consistent with low-$z$ studies, and a noncausal cosmic assembly history where mergers contribute to the statistical averaging toward the local relation is still feasible. We highlight

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

<footnote n="32">³² NPP Fellow.</footnote>
improvements to come with larger samples from JWST and, particularly, Euclid, which will exceed the statistical power of current wide and deep surveys.

*Unified Astronomy Thesaurus concepts: Supermassive black holes (1663); Active galactic nuclei (16); AGN host galaxies (2017); Galaxy evolution (594)*

*Materials only available in the online version of record: machine-readable table*

**1. Introduction**

With our understanding that galaxies grow by increasing their stellar mass through mergers and in situ star formation from gas accretion, there are still many unresolved questions in galaxy evolution. One of the most important challenges in galaxy formation is understanding the physical processes that relate the growth of supermassive black holes (SMBHs) alongside the growth of the galaxies that harbor them. Observational studies, mainly in the local Universe, have unveiled tight correlations between the mass of SMBHs ($M_{BH}$) and the physical properties of their host galaxies, such as stellar velocity dispersion $\sigma_*$ and stellar mass $M_*$ (J. Magorrian et al. 1998; L. Ferrarese & D. Merritt 2000; A. Marconi & L. K. Hunt 2003; N. Häring & H.-W. Rix 2004; K. Gültekin et al. 2009; A. W. Graham et al. 2011; A. Beifiori et al. 2012; J. Kormendy & L. C. Ho 2013; A. E. Reines & M. Volonteri 2015). How and through what physical processes such a tight relation is formed are still unclear, and the origin of the mass relation can shed light on the evolution of not only SMBHs but also galaxies.

A widely considered scenario, in answer to this question, is a coevolution scheme, where galaxies and black holes mutually increase their mass at a correlated pace. As a potential physical cause for a coevolution scenario, some studies implement active galactic nuclei (AGNs) feedback, where the energy released from AGNs heats the gas and controls star formation or gas accretion through radio jets or AGN winds (e.g., V. Springel et al. 2005; P. F. Hopkins et al. 2008; T. Di Matteo et al. 2008; A. C. Fabian 2012; C. DeGraf et al. 2015; C. M. Harrison 2017). Additionally, studies support a common gas supply simultaneously fueling both SMBHs and their host galaxies by increasing the black hole (BH) accretion rate and star formation rate (SFR; R. Cen 2015; N. Menci et al. 2016). On the other hand, others have shown that, even in the absence of a close physical connection between SMBHs and host galaxies, the mass relation can be achieved through a noncausal connection; major mergers have averaged the mass relation statistically (cosmic averaging scenario; C. Y. Peng 2007; M. Hirschmann et al. 2010; K. Jahnke & A. V. Macciò 2011).

To unravel the cause of the mass relation, an effective approach is observing the relation between $M_{BH}$ and $M_*$ throughout cosmic history. With such observations, we can directly determine whether the relation in the local Universe, both its ratio and dispersion, evolves with redshift. Then, comparisons of the observational results with simulations (e.g., X. Ding et al. 2020a; M. Habouzit et al. 2021; X. Ding et al. 2022b) based on various physical models can allow us to discuss the physical processes that establish the galaxy–BH relations and further constrain the physics of black hole formation and galaxy evolution.

Before the advent of James Webb Space Telescope (JWST), statistical studies using 2D image decomposition analyses were conducted using images obtained by Hubble Space Telescope (HST; e.g., C. Y. Peng et al. 2006; K. Jahnke et al. 2009; V. N. Bennert et al. 2011a; M. Cisternas et al. 2011;

B. D. Simmons 2011; B. D. Simmons et al. 2012; M. Schramm & J. D. Silverman 2013; M. Mechtley et al. 2016; X. Ding et al. 2020b; V. N. Bennert et al. 2021; J. I. H. Li et al. 2023b) and ground-based telescopes including Subaru's Hyper Suprime-Cam (HSC; T. Ishino et al. 2020; J. Li et al. 2021), and Pan-STARRS1 (M.-Y. Zhuang & L. C. Ho 2023). In summary, these studies have concluded that the relation between $M_{BH}$ and $M_*$ does not evolve with redshift at $z \lesssim 2$. However, studies using a large statistical and universal sample have yet to be achieved at $z > 1$, and the possibility of the redshift evolution is still left (e.g., D. Park et al. 2015). At these redshifts, we can obtain information longer than the 4000 Å break to constrain $M_*$ from observations at near-infrared wavelengths. However, there are a limited number of previous studies using near-infrared data at $z \gtrsim 1$ (e.g., K. Jahnke et al. 2009; B. D. Simmons et al. 2012; M. Mechtley et al. 2016; X. Ding et al. 2020b, utilizing HST/NICMOS or HST/WFC3). Other studies carry out statistical AGN samples using a spectral energy distribution (SED) fitting based 1D decomposition method (e.g., A. Merloni et al. 2010; M. Sun et al. 2015; H. Suh et al. 2020).

JWST is now revolutionizing the field of AGN-host galaxy relations up to $z \sim 6$ and beyond based on its high spatial resolution and unprecedented sensitivity. For instance, X. Ding et al. (2022a) applied a 2D decomposition analysis on early JWST NIRCam data from CEERS that successfully detected the host galaxies of five quasars at $z \sim 1.6-3.5$ from the SDSS DR17Q catalog (B. W. Lyke et al. 2020). They also succeeded in detecting clear substructure and performed pixel-by-pixel SED fitting for one of the five targets, SDSSJ1420+5300A, at $z \sim 1.6$. Other studies have also performed 2D decompositions of AGN host galaxies using JWST imaging; for instance, J. Li et al. (2023a) have analyzed a galaxy that is one of the most promising candidates for having a recoiling SMBH ($z \sim 0.36$) while D. D. Kocevski et al. (2023) present the host properties of five X-ray-luminous AGNs ($3 < z < 5$) in CEERS. Y. Zhang et al. (2023) also utilized JWST NIRCam data to assess the validity of $M^*$ estimated from 1D decomposition (spectrum-based) method for the HETDEX type-I AGNs ($2 < z < 2.5$).

At $z > 6$, X. Ding et al. (2023) conducted decomposition analysis of two low-luminosity quasars, thus representing the highest-redshift record for detection of host stellar emission. They suggest that $z \sim 6$ low-luminosity quasars have a mass relation consistent with the local relation after considering selection biases and measurement uncertainties, albeit with a small sample. Equally remarkable, JWST studies of high-$z$ AGNs are revealing a higher abundance of lower mass black holes that are actively accreting within very dusty and compact objects (e.g., M. Onoue et al. 2023; D. D. Kocevski et al. 2023; J. Matthee et al. 2023; Y. Harikane et al. 2023; R. Maiolino et al. 2023; J. E. Greene et al. 2023; H. B. Akins et al. 2024; D. D. Kocevski et al. 2024; T. S. Tanaka et al. 2024a, 2024b; T. S. Tanaka et al. 2024; B. Wang et al. 2024). However, in general, the sample sizes of these high-$z$ AGNs and the accuracy of $M_*$ estimation are still limited. Therefore, the redshift evolution of the mass relation remains highly uncertain.

Therefore, as an important next step to investigate the evolution of the mass relation, we perform a 2D decomposition
analysis with JWST/NIRCam data for a sample whose redshift range significantly improves upon what was statistically analyzed before the JWST era, further bridging the gap between low-z statistical and limited high-z studies. We use a sample larger than previous studies using JWST that reaches up to $z \sim 2-3$ when AGN and star formation activities peaked in cosmic history. With NIRCam images of $N = 107$ AGNs in COSMOS-Web (C. M. Casey et al. 2023) and Public Release IMaging for Extragalactic Research (PRIMER)-COSMOS, we conduct 2D decomposition analyses and then statistically discuss the evolution of the $M_{BH}-M_{*}$ relation with consideration of the selection bias and measurement uncertainty (T. R. Lauer et al. 2007; Y. Shen & B. C. Kelly 2010; A. Schulze & L. Wisotzki 2011, 2014). In addition, we report on the ability to accurately model the JWST point-spread function (PSF) in each band and the impact on the derived host galaxy properties.

This paper is organized as follows. Section 2 describes the JWST/NIRCam data and the sample selection. Section 3 describes the detailed analysis method including 2D image decomposition and careful PSF modeling, SED fitting, and construction of mock data. In Section 4, we present our fitting results and discuss the PSF effect on the results. Then, we show the evolution of $M_{BH}-M_{*}$ relation in Section 5 with considering the selection bias. Also, we discuss the possibility of scatter evolution in the mass relation and summarize the challenges with 2D decomposition methods in Section 6. We present the conclusion and prospects for future studies in Section 7. In this paper, all magnitudes are AB magnitude (J. B. Oke 1974), and we assume a standard cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_m = 0.30$, and $\Omega_\Lambda = 0.70$.

## 2. Data

### 2.1. COSMOS-Web

COSMOS-Web (PI: Jeyhan Kartaltepe and Caitlin Casey, GO1727, see C. M. Casey et al. 2023, for the overview) is a 270 hr treasury survey program in JWST Cycle 1, covering $0.54 \text{ deg}^2$ with NIRCam (M. J. Rieke et al. 2023) in four filters (F115W, F150W, F277W, F444W) and $0.19 \text{ deg}^2$ with MIRI (P. Bouchet et al. 2015) using F770W. Due to the large field, the COSMOS-Web field was split into 20 tiles.

The data are reduced with the JWST Calibration Pipeline<fnref n="33" /> (H. Bushouse et al. 2023) version 1.14.0 and the calibration Reference Data System version 1223. The $5\sigma$ depth in an aperture with a radius of $0.0'15$ ranges from 26.7 to 27.5 mag in F115W and 27.5 to 28.2 mag in F444W, depending on the number of integrations (also see Section 2.1 of C. M. Casey et al. 2023). The mosaic images have a resolution of $0.0'030 \text{ pixel}^{-1}$. Details of the image reduction process will be described in M. Franco et al. (2025, in preparation). In addition to NIRCam four-band data, we use HST/ Advanced Camera for Surveys (ACS) F814W data (A. M. Koekemoer et al. 2007). In this study, we do not use MIRI data because it is challenging to apply the 2D decomposition method (Section. 3.2) due to its lower spatial resolution and larger PSF.

### 2.2. PRIMER

PRIMER (PI: James Dunlop, GO1837) is a 195 hr treasury program of JWST Cycle 1, targeting two equatorial HST

CANDELS Legacy Fields: COSMOS and Ultra Deep Survey (UDS). PRIMER-COSMOS covers 144 arcmin² with eight NIRCam (M. J. Rieke et al. 2023) filters (F090W, F115W, F150W, F200W, F277W, F356W, F410M, F444W) and 112 arcmin² with two MIRI (P. Bouchet et al. 2015) filters (F770W and F1800W) in the COSMOS field. The processed COSMOS-PRIMER data consist of one mosaic image. In this analysis, we use NIRCam eight-band data and HST/ACS F814W data. The PRIMER-COSMOS data are reduced with the JWST Calibration Pipeline (H. Bushouse et al. 2023) version 1.8.3 and the calibration Reference Data System version 1017. The $5\sigma$ depth in an aperture with a radius of $0.0'15$ has a wide range from 27.9 to 28.3 mag in F090W and 28.4 to 28.9 mag in F444W, depending on the number of integrations, ~1 mag deeper than the COSMOS-Web data. The mosaic images have a resolution of $0.0'030 \text{ pixel}^{-1}$.

### 2.3. Broad-line AGN Sample

To evaluate the relation between $M_{*}$ and $M_{BH}$, we use the type-I AGN sample with $M_{BH}$ estimates available in A. Schulze et al. (2015, 2018). They fitted broad emission lines from the data of large spectroscopic surveys and estimated $M_{BH}$ with virial mass estimators. A. Schulze et al. (2015) presents the redshift evolution of AGN population based on spectroscopically observed type-I AGNs from zCOSMOS (S. J. Lilly et al. 2007; S. J. Lilly et al. 2009), VVDS (O. Le Fèvre et al. 2005, 2013; B. Garilli et al. 2008), and SDSS (D. P. Schneider et al. 2010; Y. Shen & B. C. Kelly 2012). A. Schulze et al. (2018) provides the properties of X-ray selected and spectroscopically confirmed type-I AGNs in the FMOS-COSMOS survey (D. Kashino et al. 2013; J. D. Silverman et al. 2015). Here, we select the targets in A. Schulze et al. (2015) and A. Schulze et al. (2018) that are also detected by Chandra (Chandra-COSMOS Survey; M. Elvis et al. 2009; F. Civano et al. 2012) or XMM-Newton (XMM-COSMOS; N. Cappelluti et al. 2009; M. Brusa et al. 2010). The 2-10 keV flux sensitivity is $7.3 \times 10^{-16} \text{ erg cm}^{-2} \text{ s}^{-1}$ for Chandra and $3 \times 10^{-15} \text{ erg cm}^{-2} \text{ s}^{-1}$ for XMM-Newton.

By combining A. Schulze et al. (2015) and A. Schulze et al. (2018), we have $M_{BH}$ estimated from spectral data that are acquired by the FMOS-COSMOS and zCOSMOS surveys with some AGNs having measurements from both. Considering the quality of the spectroscopy, the error on $M_{BH}$ estimation, and the fact that $H\beta$ line is used for calibrating virial mass estimators, we use the $M_{BH}$ from $H\beta$ (FMOS-COSMOS), $H\alpha$ (FMOS-COSMOS), and MgII (zCOSMOS-Bright, Deep) in order of preference; e.g., FMOS $H\beta$ is used for an object with both FMOS $H\beta$ and zCOSMOS MgII estimation. As shown in Figure 10 of A. Schulze et al. (2018), there is a very good agreement between the FMOS $H\alpha$- and $H\beta$-based $M_{BH}$ compared to those using MgII. A. Schulze et al. (2018) also compared full width at half-maximum (FWHM) measurements taken at different times for each object (see Figure 7 of A. Schulze et al. 2018) and confirmed that the FWHM values are consistent with each other. Note that 19 AGNs in our sample are listed in the SDSS DR16 quasar catalog (B. W. Lyke et al. 2020); we do not use these SDSS MgII-based $M_{BH}$ measurements given the benefits of the FMOS and deeper zCOSMOS spectroscopy. The number and redshift range of each measurement are summarized in Table 1.

For $H\beta$, A. Schulze et al. (2018) used the virial mass estimation relation by M. Vestergaard & B. M. Peterson

<footnote n="33">³³ https://github.com/spacetelescope/jwst</footnote>
<table><caption>**Table 1**
Sample Size for Each Single-epoch $M_{BH}$ Estimation</caption><thead><tr><th>Survey</th><th>Line</th><th>COSMOS-Web</th><th>PRIMER</th><th>z Range</th></tr></thead><tbody><tr><td>FMOS-COSMOS</td><td>Hα</td><td>52</td><td>3</td><td>0.68–1.7</td></tr><tr><td>...</td><td>Hβ</td><td>21</td><td>2</td><td>1.2–2.5</td></tr><tr><td>zCOSMOS-bright</td><td>Mg II</td><td>26</td><td>0</td><td>1.1–2.1</td></tr><tr><td>zCOSMOS-deep</td><td>Mg II</td><td>3</td><td>0</td><td>0.94–1.2</td></tr><tr><td>Total</td><td>...</td><td>102</td><td>5</td><td>0.68–2.5</td></tr></tbody></table>

(2006),

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

where $L_{5100}$ is continuum luminosity at 5100 Å, and $FWHM_{H\beta}$ is FWHM of Hβ broad line. Then, the Hα-based masses are calculated as given in A. Schulze et al. (2017, 2018),

$$M_{BH}(H\alpha) = 10^{6.71} \left( \frac{L_{H\alpha}}{10^{42} \text{ erg s}^{-1}} \right)^{0.48} \times \left( \frac{\text{FWHM}_{H\alpha}}{1000 \text{ km s}^{-1}} \right)^{2.12} M_{\odot}, \quad (2)$$

where $L_{H\alpha}$ is the Hα luminosity, and $FWHM_{H\alpha}$ is the FWHM of the broad Hα emission line. For MgII-based $M_{BH}$ estimation (zCOSMOS-bright MgII, and zCOSMOS-deep MgII), the calibration by Y. Shen et al. (2011) is used,

$$M_{BH}(MgII) = 10^{6.74} \left( \frac{L_{3000}}{10^{44} \text{ erg s}^{-1}} \right)^{0.62} \times \left( \frac{\text{FWHM}_{MgII}}{1000 \text{ km s}^{-1}} \right)^2 M_{\odot}, \quad (3)$$

where $L_{3000}$ is continuum luminosity at rest-frame 3000 Å, and $FWHM_{MgII}$ is FWHM of MgII broad emission line. These single-epoch virial mass estimations have an uncertainty due to possible variability and uncertainties in the modeling of broad-line regions (see., Y. Shen 2013). In this paper, we use the $M_{BH}$ uncertainties that also consider uncertainties from the single-epoch virial mass estimation, typically ~0.4 dex. We also consider this uncertainty in generating mock data (Section 3.6).

From the parent catalog, we select broad-line AGNs that reside in the COSMOS-Web and PRIMER fields. The final sample size has 107 AGNs with black hole mass estimation summarized in Table 1. Figure 1 shows the spatial location of the AGN sample within the COSMOS-Web footprint. We use five broadband images from HST/ACS (F814W) and JWST/NIRCam (F115W, F150W, F277W, and F444W) for the sample residing in COSMOS-Web. Four more broadband (F090W, F200W, F356W) and medium-band (F410M) images are available for five galaxies in the PRIMER field.

We also compare the optical color index $g-i$ calculated based on the COSMOS2020 photometry (J. R. Weaver et al. 2023) with SDSS quasars and hard-X-ray-detected AGNs at the same redshift range (Figure 6 in J. D. Silverman et al. 2005). Our sample has a wide $g-i$ distribution of $g-i = 0-2.4$. Thus, our sample includes both unobscured and dust-obscured AGNs and is not significantly biased to either sample.

Figure 2 shows the distribution of 2–10 keV X-ray luminosity $L_{[2–10\text{keV}]}$ (panel (a)) and $M_{BH}$ (panel (b)) as a function of redshift, where $L_{[2–10\text{keV}]}$ is calculated with an

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Location of the sample compared to the footprints of COSMOS-Web and PRIMER-COSMOS. The background is the mosaic image of HST/ACS F814W (A. M. Koekemoer et al. 2007) for the COSMOS field. Red stars show the position of each AGN. The region enclosed by the solid cyan line shows the PRIMER-COSMOS field, and gray shaded region corresponds to the COSMOS-Web field.</figcaption></figure>

X-ray spectral index $\Gamma = 1.8$ (e.g., M. Brightman et al. 2013). Since the sample consists of X-ray-selected objects and is flux limited, there is a tendency for higher-$z$ objects to have larger $L_{[2–10\text{keV}]}$ over the sensitivity limit. We can also see the trend of $M_{BH}$ increasing with redshift. The sample biases likely influence this trend from the flux sensitivity and the availability of broad-line FWHM measurements from spectroscopic data. Figure 2(c) displays the relation between the Eddington ratio $L_{bol}/L_{edd}$ and $M_{BH}$ with Eddington ratio decreasing as $M_{BH}$ increases. This trend is due to the observational flux limitation and that $L_{edd}$ is proportional to $M_{BH}$. While there is a correlation between $M_{BH}$ and $L_{bol}$, dividing $L_{bol}$ by $L_{edd}$ to calculate the Eddington ratio cancels out this weaker correlation between $L_{bol}$ and $M_{BH}$.

Figure 3 shows the original F277W images, i.e., before the decomposition analysis, of representative AGNs in our sample. In some targets, we can recognize the extended components of the host galaxies far from the central PSF-like feature. However, a central AGN component, especially those with spiky diffraction features in the outer part, dominates the system and buries the host galaxy component. These dominant PSF components prevent us from obtaining host galaxy information directly and make the 2D decomposition analysis necessary.

### 3. Method

To extract host galaxy components from the original AGN + host galaxy NIRCam images, we apply a 2D image analysis tool galight (X. Ding et al. 2020b). With galight, we perform forward modeling of each image as a superposition of a PSF component and PSF-convolved Sérsic components corresponding to the light from an AGN and its host galaxy, respectively. We then obtain images of the host galaxy, free
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Characteristics of the type-I AGN sample with JWST imaging from COSMOS-Web and PRIMER: (a) Distribution of 2–10 keV X-ray luminosity $L_{[2–10keV]}$ as a function of z. Black dashed lines indicate the $F_{[2–10keV]}$ sensitivity in the sample, (b) black hole mass $M_{BH}$ as a function of z, and (c) relation between the Eddington ratio $L_{bol}/L_{edd}$ and $M_{BH}$. The color and shape of the points indicate the JWST survey field and the source of $M_{BH}$ estimation, respectively, as shown.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Original images in F277W of some targets in the order of redshift. The target IDs and the redshifts are shown in the left corner of each image. The white bars indicate a 1" scale. Depending on the host-to-total flux ratio ($H/T$) that varies with the sample, a central AGN component can dominate an entire system, thus burying a host galactic component.</figcaption></figure>

of the AGN, by subtracting the fitted PSF component from the original image.

### 3.1. Comparative Analysis of Model PSF Construction

Considering that the AGN can account for up to ~95% of the total flux (e.g., X. Ding et al. 2020b; X. Ding et al. 2022a), the results significantly depend on the accuracy of reconstructing the PSF. There are different strategies to reconstruct PSF images based on either using theoretical PSFs (e.g., K. A. Suess et al. 2022) or stellar images (e.g., D. Nardiello et al. 2022; X. Ding et al. 2022a; X. Ding et al. 2023; M.-Y. Zhuang & Y. Shen 2023; W. M. Baker et al. 2023). The former uses the theoretical PSF model such as WebbPSF (M. D. Perrin & R. Soummer 2012; M. D. Perrin et al. 2014), and the latter uses natural stellar images as the PSF directly or modeled PSF with tools such as PSFEx (E. Bertin 2011). Many previous studies concluded that the synthetic PSF simulated by WebbPSF is narrower than the PSFs reconstructed with the natural stars (Y. Ono et al. 2022; X. Ding et al. 2022a; M. Onoue et al. 2023). Note that K. Ito et al. (2023) used the intermediate

method between the former and the latter; they used theoretical PSFs with WebbPSF and smoothed it by comparing the surface brightness profile of natural stellar images. In this paper, we reconstruct the PSF using three methods based on natural star images for each region and filter. Then, we compare the results with the different PSF reconstruction methods and discuss the host galaxy characteristics with the method dependence.

#### 3.1.1. $\chi^2$-based Methods

First, we follow the strategy of X. Ding et al. (2020b) and X. Ding et al. (2022a). They first construct an empirical PSF library for which each PSF is represented by the image of a single star. Then, the 2D decomposition analysis is run with all single PSFs in the library, separately. Then, they sort the results in the order of reduced chi-square $\chi^2$ and stack the PSFs with the top 3, 5, and 8 $\chi^2$ values. Using the single PSFs and the stacked PSFs, they select the PSF with the smallest $\chi^2$ as the final PSF. This method is based on the $\chi^2$; i.e., they assumed that the lower $\chi^2$ is indicative of a better (the closer to the more accurate) PSF.

Following their strategy, we apply the find_PSF function in galight to list PSF candidates, then select PSF candidates manually for each mosaic image and filter. In this manual selection process, only obvious PSF candidates with the PSF-like complex hexagonal and spiky diffraction features and without a galaxy-like broad component are selected. We cropped the images of the selected PSF candidates for the short-wavelength-channel filters (F090W, F115W, F150W, and F200W) and long-wavelength-channel filters (F277W, F356W, F410M, and F444W) as squares with 150 and 240 pixels per side, corresponding to 4"/5 and 7"/2, respectively. After removing neighbor objects using clean_PSF function in galight, the PSF libraries contain ~30–50 PSF candidates depending on the filter and region. Then, we fit each AGN target with a superposition of a PSF-convolved Sérsic profile and each single PSF candidate in the PSF library. With the fitting results of each single PSF, we select PSFs with the top-5 $\chi^2$ and top-75% $\chi^2$ and stack them to generate an averaged PSF image. These top-5 stacked and top-75% stacked PSFs are finally used to estimate the parameters in this method. Note that each target has its own top-5 and top-75% PSFs generated from single PSFs with the lowest $\chi^2$ selected for each target.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison of FWHM (semimajor axis) for PSFs from different PSF reconstruction methods for all individual AGNs. The x-axis represents the ratio of the FWHMs for the top-5 PSFs to those for the PSFEx PSFs, and the y-axis represents the ratio of the FWHMs for the top-75% PSFs to those for the PSFEx PSFs. Blue, green, orange, and red colors indicate the different filters: F115W, F150W, F277W, and F444W. Median values and $1\sigma$ confidence range are denoted by star symbols and error bars. The gray dashed line indicates $y=x$, i.e., the same FWHMs for the top-5 and top-75% PSFs. Notably, F277W and F444W exhibit an FWHM bias among different PSF reconstruction methods.</figcaption></figure>

### 3.1.2. Modeling Method

M.-Y. Zhuang & Y. Shen (2023) compare JWST/NIRCam PSFs modeled with different methods (Swarm, photutils, and PSFEx), and concluded that PSFEx reconstructed PSFs provide the best performance. From the 2D decomposition of simulated broad-line AGNs, M.-Y. Zhuang & Y. Shen (2023) also suggested that smaller $\chi^2$ values do not necessarily provide a means to distinguish which PSFs are more likely to characterize the AGN with higher accuracy. Following the conclusion by M.-Y. Zhuang & Y. Shen (2023), we use PSFEx and compare the results with $\chi^2$-based selected PSFs (Section 3.1.1).

PSFEx constructs an empirical PSF model based on the output catalog of SExtractor (E. Bertin & S. Arnouts 1996). We first run SExtractor for source detection, and then run PSFEx for modeling the PSF for each mosaic image and filter. PSFEx can also reconstruct local PSFs as a function of positions on the detector. We do not use local PSFs because M.-Y. Zhuang & Y. Shen (2023) also concluded that a universal or global PSF usually shows “satisfactory” fitting results, and the sample region (COSMOS-Web and PRIMER-COSMOS) has a much smaller number of stars than the south continuous viewing zone, which M.-Y. Zhuang & Y. Shen (2023) tested local PSF reconstruction.

### 3.1.3. Comparing the Final PSFs

Now, we have three final PSFs for comparison; the top-5 $\chi^2$ stacked PSF, the top-75% $\chi^2$ stacked PSF, and the PSFEx PSF (Sections 3.1.1 and 3.1.2). To compare PSFs, we perform, a 2D Gaussian fitting for each PSF image and measure the FWHMs along the semimajor axis. Figure 4 compares the FWHMs of PSFs obtained with each method, target, and filter. The x- and y-axes of Figure 4 show the ratio $FWHM_{\text{Top-5}}/\text{FWHM}_{\text{PSFEx}}$ and $FWHM_{\text{Top-75\%}}/\text{FWHM}_{\text{PSFEx}}$, where FWHM is the value along the semimajor axis.

First, regardless of the filters, we can see that the distribution extends farther in the x-axis direction than the y-axis. This can be attributed to greater variation in the FWHMs for the top-5 stacked PSFs. We use the same PSFEx PSF for each target in

the same field, and the top-75% stacking in the same field uses mostly the same single PSFs in the field. In contrast, top-5 stacking employs only the best-fit single PSFs with the lowest $\chi^2$. As a result, the FWHM variation for each galaxy is largest for the top-5 PSF followed by the top-75% PSF and PSFEx PSF, and $FWHM_{\text{Top-5}}/\text{FWHM}_{\text{PSFEx}}$ have a larger scatter than $FWHM_{\text{Top-75\%}}/\text{FWHM}_{\text{PSFEx}}$.

Second, we focus on the FWHM bias between the methods for each filter. As suggested by M.-Y. Zhuang & Y. Shen (2023), for short-wavelength filters (F115W and F150W), there is a significant scatter in the FWHM ratio. On the other hand, for long-wavelength filters (F277W and F444W), the scatter is smaller than the short-wavelength side. These results imply that, in the long-wavelength filters, PSFEx PSFs are sharper than the top-5 and Top-75% PSFs. The impact of these trends on the 2D decomposition analysis is discussed in Section 4.2.

Note that these trends can depend on the visual inspection performed when constructing the PSF library (Section 3.1.1) and the settings used for Sextractor and PSFEx (Section 3.1.2). For example, visual inspection can be biased by the hexagonal diffraction features of the JWST PSF, an appropriate FWHM range, and the absence of extended structures originating from host galaxies. If this selection process is strongly biased by the hexagonal features, it might lead to a selective choice of brighter PSFs. As a result, the parameter distributions presented here may not necessarily match the distribution of the actual PSF. Nonetheless, even different PSF reconstruction methods can result in different FWHMs. Therefore, when performing a 2D decomposition analysis with only one PSF reconstruction method and not considering the possibilities of other PSFs, 2D decomposition results can be biased by a specific PSF. Because determining the PSF shape perfectly is challenging, it is also important to discuss uncertainties by considering the results obtained with possible different PSFs. We discuss how different PSF reconstruction methods affect the results of the 2D decomposition and the final $M_*$ estimation in Section 4.2, and we perform a detailed comparison of the obtained final PSFs in Appendix B.

### 3.2. Decomposition

Using galight, we fit the AGN + host galaxy images with the composite model of a PSF component and a single Sérsic profile (J. L. Sérsic 1963; J. L. Srsic 1968) convolved by a PSF. Note that we do not assume the same morphology in every band, i.e., the fitting is performed in each filter independently. In cases where there are nearby galaxies that can affect the fitting, these galaxies are also modeled as Sérsic components and fitted simultaneously. Our Sérsic model has seven free parameters: amplitude, Sérsic index *n*, effective radius $r_e$, coordinates of the center $x_c$, $y_c$, and ellipticity $e_1, e_2$. The PSF model corresponding to the AGN component has three free parameters: amplitude and coordinates of the center $x_c$, $y_c$. Thus, the number of free parameters for PSF + single Sérsic component is 10 in total. Note that the actual number of free parameters in the fitting changes depending on the number of nearby objects also fitted with a Sérsic profile.

To avoid unphysical results, the range of *n* and $r_e$ is constrained to [0.3, 7] and [0.06, 2.0], respectively. Note that some galaxies show clear substructures that do not suit a Sérsic profile, such as bars and spiral arms. Thus, using a Sérsic profile is a first-order approximation to model the global component of AGN host galaxies.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Example fits for CID-62 ($z \sim 1.92$) using F444W. Each row shows the result for a different PSF (top-5 stacked, top-75% stacked, and PSFEx, from top to bottom). In each row, the images are as follows: original, model, data-model point source (host galaxy only), and residuals normalized by the error map, from left to right. $\chi^2$ values are shown in the panels of data-model point source. The right panel shows the 1D surface brightness profile where dashed lines indicate half-width at half-maximum (HWHM) of each PSF. In the data point source image, we reveal a disklike host galaxy, which is buried under the PSF before subtraction, regardless of the PSF reconstruction method.</figcaption></figure>

In the fitting process, we cut the image into square regions centered on the target with a radius 7 times the standard deviation along the semimajor axis of the 2D Gaussian fitting with photutils (L. Bradley et al. 2023). Then, the above model is optimized with a particle swarm optimizer (J. Kennedy & R. Eberhart 1995). galight also supports Markov Chain Monte Carlo (MCMC) to estimate the posterior parameter distributions. As suggested by X. Ding et al. (2023), we also confirm that the uncertainty estimated with MCMC is much smaller than the uncertainties from different PSFs. Thus, we do not use MCMC in the decomposition analysis, and we estimate the uncertainty from the results with different single PSFs (Section 3.4). We set the supersampling factor relative to pixel resolution to 3, which controls interpolation within a pixel to perform a subpixel shift of the PSF (see., X. Ding et al. 2023).

As described above, using the PSF library constructed with galight, we fit with every single PSF in the library and sort them in order of $\chi^2$. Then, we fit with the three final PSFs: top-5 $\chi^2$ stacked, top 75% $\chi^2$ stacked, and PSFEx. Figure 5 compares the fitting results with the final PSFs for the example galaxy CID-62 in the F444W. We can find that the host component is more prominent than the central PSF component at larger radii, and 2D decomposition makes it possible to detect the host galaxy initially buried under the PSF component with all three model PSFs.

For targets that have $n > 6.5$ in any band other than F277W, we rerun the fit while fixing $n$ to the value found for the F277W band. F277W has the lowest number of values hitting the upper limit on $n$, falls above the rest-frame 4000 Å break, and is somewhat central to JWST wavelength coverage. This pertains to 42, 30, and 26 sources detected in F115W, F150W, and F444W, respectively. From mock tests, we confirm that galight can return $n \sim 7$ even if the actual value is much smaller (Appendix C). This is likely due to fitting where some of the AGN emission is attributed to a central stellar concentration of the host thus overestimating the host galaxy flux as well.

### 3.3. Detection of Host Galaxy

Figure 5 shows an example where the host galaxy is clearly detected. However, for some galaxies, the strong PSF component dominates the total flux maybe due to not only a low host-to-total flux ratio $H/T$ but also compact morphology, making it challenging to distinguish the host signal from their PSF component. To gauge, in a quantitative manner, which AGNs have accurate host galaxy information, we exclude cases where the host galaxy is undetectable, following three strategies below.

1. *Bayesian Information Criteria.* In addition to the PSF + Sérsic model (PS+SE model) described above, we also fit with a model containing only a PSF component
(PS model). Then, we calculate the Bayesian Information Criteria (BIC; G. Schwarz 1978) for the two models, PS +SE and PS model, as

$$BIC = \chi^2 + k \ln(n), \quad (4)$$

where k is the number of free parameters, and n is the number of data points. We regard that the PS+SE model provides a better description of the data than the PS model when BIC<sub>PS+SE</sub> is much smaller than BIC<sub>PS</sub>, as

$$BIC_{PS+SE} < BIC_{PS} - 10. \quad (5)$$

We decide the threshold value of the BIC difference of 10 based on R. E. Kass & A. E. Raftery (1995).

2. *S/N of the host galaxy.* To estimate the significance of the detection of the host galaxy, we calculate the signal-to-noise ratio (S/N) of the host galaxy as done in X. Ding et al. (2023). We construct an error map of the PSF-subtracted images considering two sources of error: noise from the observed images and the uncertainty propagated from different PSF reconstructions. Due to the intrinsic variations of the PSF image even in the same field of view (M. Yue et al. 2023; M.-Y. Zhuang & Y. Shen 2023) and errors in the observed image used in PSF reconstruction, the reconstructed PSFs contain uncertainties. Thus, considering the uncertainty from PSF reconstruction is indispensable to calculate the S/N of detected host galaxies. For the PSF uncertainty, we use the pixel-by-pixel standard deviation of the fitted PSF components when fitting the PS+SE model with each single PSF in the PSF library. Then, we calculate the final noise map as a composite of the observed noise map and the PSF uncertainty map. With the final noise map, we calculate the S/N of host galaxy S/N<sub>host</sub> within a radius of 2r<sub>e</sub>. Then, we define the detection as cases with high S/N<sub>host</sub>, as

$$S/N_{\text{host}} > 5. \quad (6)$$

3. *Manual inspection and removal.* We find that some objects have invalid central values (e.g., masked due to saturation) in their F444W image and shallower surface brightness profiles than any PSFs. The fitting of these galaxies fails even with applying a mask in the central region. We find four cases (CID-50, CID-208, CID-668, and CID-112) with such features and label them as nondetections.

CID-142 is located near the edge of the image, and the host galaxy is partially cut off. It is also the possible that CID-142 has a mismatched PSF from PSFs in other fields. Therefore, obtaining accurate photometry of CID-142 is challenging; thus, we manually exclude CID-142 in the following discussion.

We also find some obvious false detections in F814W, where the image shows a dominant PSF feature and no extended host-like feature, and galight fit the PSF-like features as a host galaxy with H/T ~ 1. We confirm that the decomposition of the JWST images clearly detects a host-like extended feature. Thus, in addition to the above strategies, we recognize 31 obvious false detections in F814W as nondetection cases.

With the above three strategies, we confidently report the detection of an AGN host galaxy for those that fulfill both the conditions of Equations (5) and (6) for each of the final PSFs; i.e., we determine whether the host is detected or not for each top-5, top-75%, and PSFEx PSF, separately. Also note that this

decision is made for each sample and each filter, i.e., a galaxy detected in one filter may not be detected in another filter. The number of objects detected over two filters of NIRCam is 102, occupying ~95% out of the entire N = 107 sample with the top-5 PSF (see Appendix A for the number of detection in each filter).

### 3.4. Photometry of Host Galaxies

We calculate the flux of the host galaxies, using the Sérsic fits, considering Galactic dust extinction (D. J. Schlegel et al. 1998) for the detected cases. For the photometric accuracy, we set an error of 0.2 mag, which represents likely systematic uncertainties (e.g., X. Ding et al. 2022a; M.-Y. Zhuang & Y. Shen 2023; M.-Y. Zhuang et al. 2023) and errors discussed in Section 3.3, considering both observational errors and PSF uncertainty in a radius of 2r<sub>e</sub>. We use the 3σ value for the undetected filters as the upper limit.

Here, we also fit with a model containing only a Sérsic component (SE model). For targets with BIC<sub>PS+SE</sub> > BIC<sub>SE</sub> - 10, the SE model is better or flexible enough to describe the data than the PS+SE model, and we use the Sérsic photometry calculated with the SE model instead of the PS+SE model. Such cases are observed only in F444W and are very limited (two or three objects depending on the PSF).

### 3.5. SED Fitting

We fit the photometry of the host galaxies with CIGALE (v2022.1, M. Boquien et al. 2019; G. Yang et al. 2022) SED fitting library. For a stellar population, we use the single stellar population model by G. Bruzual & S. Charlot (2003; bc03 module) and the Chabrier initial mass function (IMF; G. Chabrier 2003) with the M<sub>*</sub> cutoff of 0.1M<sub>☉</sub> and 100M<sub>☉</sub> (G. Bruzual & S. Charlot 2003). We assume a delayed-τ model for a star formation history (SFH), where SFR at each look-back time t is modeled as

$$SFR(t) \propto \begin{cases} (t - t_{\text{age}}) \exp\left(-\frac{t - t_{\text{age}}}{\tau}\right) & (t > t_{\text{age}}), \\ 0 & (t < t_{\text{age}}), \end{cases} \quad (7)$$

where t<sub>age</sub> and τ indicate the starting time of star formation activity and the declining timescale of SFR. We also consider a nebular emission with nebular module and a dust attenuation with dustatt_modified_starburst module that assumes the modified D. Calzetti et al. (2000) law. The free parameters in the CIGALE fitting are E(B - V), M<sub>*</sub>, t<sub>age</sub>, and τ. Their grid values are decided basically following M.-Y. Zhuang et al. (2023) and summarized in Table 2. To avoid unphysical solutions, we set the upper limit of t<sub>age</sub> to 0.95t<sub>H</sub>, where t<sub>H</sub> indicates the cosmic age at each redshift. Otherwise, stellar metallicity and the ionization parameter U are fixed at solar metallicity and log U = -2.

We also apply the Bayesian-based SED fitting code, Prospector (J. Leja et al. 2017; B. D. Johnson et al. 2021), to assess the uncertainty M<sub>*</sub> with different SED fitting codes. Prospector is based on the flexible stellar population synthesis (FSPS; C. Conroy et al. 2009; C. Conroy & J. E. Gunn 2010) to generate model SEDs of galaxies. For comparison, we use almost the same settings with CIGALE: a Chabrier IMF, a delayed-τ SFH model, solar metallicity, D. Calzetti et al. (2000) dust attenuation law, and nebular emission with
<table><caption>**Table 2**

Major Parameter Settings in SED Fitting</caption><thead><tr><th>Parameter</th><th>Description</th><th>Values (CIGALE)</th><th>Prior (Prospector)</th></tr></thead><tbody><tr><td>log (M<sub>*</sub>/M<sub>☉</sub>)</td><td>total stellar mass</td><td>scaled with the data</td><td>uniform: min = 9, max = 13</td></tr><tr><td>log (Z<sub>*</sub>/Z<sub>☉</sub>)</td><td>stellar metallicity</td><td>fixed at 0</td><td>fixed at 0</td></tr><tr><td>log U</td><td>ionization parameter</td><td>fixed at -2</td><td>fixed at -2</td></tr><tr><td>E(B - V)/mag</td><td>Color excess for the nebular lines</td><td>0, 0.001, 0.005, 0.01, 0.03, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7</td><td></td></tr><tr><td>τ<sub>dust,2</sub></td><td>optical depth of the diffuse dust attenuation</td><td></td><td>uniform: min=0, max=3</td></tr><tr><td>τ/Gyr</td><td>timescale of delayed-τ SFH</td><td>0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10, 12, 14, 16, 18, 20</td><td>uniform in log (τ/Gyr): min=-2, max=1.5</td></tr><tr><td>t<sub>age</sub>/Gyr</td><td>starting time of delayed-τ SFH</td><td>equally sampled in [1 Gyr, 1.95t<sub>H</sub>] with the separation of 0.1 Gyr</td><td>uniform: min=1, max=0.95t<sub>H</sub></td></tr></tbody></table>

log U = -2. Prospector can fit with the nonparametric SFH, which separates galaxy formation history into several age bins and assumes a constant SFR in each bin (e.g., J. Leja et al. 2019). S. Lower et al. (2020) input cosmological hydrodynamic simulation data into Prospector and concluded that nonparametric SFH tends to reconstruct M* more accurately than parametric SFHs. S. Lower et al. (2020) also suggested that parametric SFH tends to underestimate M*. However, in this study, the available photometry is limited in the number and wavelength range (only five/nine bands in the near-infrared wavelength range for the COSMOS-Web/PRIMER-COSMOS field). Furthermore, the photometry derived in Section 3.4 contain uncertainties from the 2D decomposition analysis. Therefore, we choose to use the parametric SFH (delayed-τ model) instead of the nonparametric assessment. The parameter prior settings in the MCMC run are summarized in Table 2. We compare the results with CIGALE and Prospector in Section 4.4.

### 3.6. Generating Mock Data to Consider Selection Effects

As mentioned, our sample is X-ray-flux limited (Section 2.3), raising the possibility of bias toward larger M<sub>BH</sub> or higher Eddington ratio (T. R. Lauer et al. 2007; A. Schulze & L. Wisotzki 2011, 2014). Due to this selection effect, a direct comparison of the observational results with the local relation is not appropriate. Thus, in this study, we generate a mock AGN-galaxy catalog based on the procedure in J. Li et al. (2021) and apply the mock observation (adding selection biases and observational effects) to discuss the intrinsic evolution of the mass relation. The procedure for generating the mock catalogs is described below.

First, we generate the mock redshift $z_{\text{mock}}$ and mock true stellar mass $M_{*,\text{mock}}$ based on the COSMOS2020 stellar mass function (SMF) by J. R. Weaver et al. (2023). Next, we use the $M_{*,\text{mock}}$ to generate the mock true BH mass $M_{\text{BH, mock}}$ for the mock sample. Here, we assumed the local relation as

$$ \log \left( \frac{M_{\text{BH}}}{M_{\odot}} \right) = \alpha \log \left( \frac{M_{*}}{M_{\odot}} \right) + \beta, \quad (8) $$

where $\alpha$ and $\beta$ indicate the slope and the intercept of the local $M_{\text{BH}}-M*$ plane. Then, assuming a normal distribution, we calculate $M_{\text{BH, mock}}$ as

$$ \log M_{\text{BH, mock}} = \mathcal{N} \left( \alpha \log \left( \frac{M_{*}}{M_{\odot}} \right) + \beta + \gamma \log(1+z), \sigma_{\mu}^{2} \right). \quad (9) $$

Here, the parameter $\gamma$ indicates the strength of redshift evolution of the mass relation, and $\sigma_{\mu}$ is the intrinsic scatter of the mass relation. Figure 6(a) shows $M_{\text{BH, mock}}$ versus $M_{*,\text{mock}}$ distribution with the example parameter set of $\gamma = 0$ and $\sigma_{\mu} = 0.3$.

Then, based on $z_{\text{mock}}$ and $M_{\text{BH, mock}}$, we assign mock Eddington ratio $\lambda_{\text{Edd, mock}}$ by sampling the Eddington ratio distribution function by A. Schulze et al. (2015). Using the $\lambda_{\text{Edd, mock}}$ and $M_{\text{BH, mock}}$, we calculate the bolometric luminosity $L_{\text{bol, mock}}$.

From $L_{\text{bol, mock}}$, we obtain $L_{[2-10\text{keV}, \text{mock}]}$ using the bolometric correction by F. Duras et al. (2020). With the calculated $L_{[2-10\text{keV}, \text{mock}]}$, the mock X-ray flux $F_{[2-10\text{keV}, \text{mock}]}$ is determined with the assumption of $\Gamma = 1.8$, the same assumption as in Section 2.3. Figure 6(b) shows $F_{[2-10\text{keV}, \text{mock}]}$ distribution with the example parameter set of $\gamma = 0$ and $\sigma_{\mu} = 0.3$.

We also generate mock virial BH mass $M_{\text{vir, mock}}$ with the $L_{\text{bol, mock}}$, $M_{\text{BH, mock}}$, and the assumed FWHM distribution. We first calculate mock continuum and line luminosity $L_{5100, \text{mock}}$, $L_{3000, \text{mock}}$, and $L_{H\alpha, \text{mock}}$. For $L_{5100, \text{mock}}$ and $L_{3000, \text{mock}}$, we used the bolometric correction by H. Netzer & B. Trakhtenbrot (2007) and B. Trakhtenbrot & H. Netzer (2012), respectively. For $L_{H\alpha, \text{mock}}$, we use the $L_{5100, \text{mock}}$ and scaling relation between $L_{H\alpha}$ and $L_{5100, \text{mock}}$ by H. D. Jun et al. (2015). These bolometric correction and scaling relations are the same as those in A. Schulze et al. (2018). Then, we assume that FWHM of the emission line follows the log-normal distribution with the scatter of 0.17 dex (Y. Shen et al. 2008), and we generate the observed FWHM of $H\alpha$, $H\beta$, and MgII emission lines following Equations (1), (2), and (3). To consider the bias in single-epoch virial mass estimation, we add the bias with $\beta_{\text{bias}} = 0.6$ in calculating FWHM. $\beta_{\text{bias}}$ represents the proportion by which $\Delta L$, the variation of luminosity from the mean luminosity $\bar{L}$, affects the variation in FWHM. Thus, the FWHM is generated by a log-normal distribution with the scatter of 0.17 dex and the mean value corresponding to the luminosity of $\bar{L} + \beta\Delta L$ (see Y. Shen 2013; J. Li et al. 2021). Finally, the mock virial BH mass $M_{\text{vir, mock}}$ is calculated using the mock observed FWHM and the luminosity following Equations (1), (2), and (3). Figures 6(c) and (d) show the FWHM $H\alpha, \text{mock}$ histogram and $M_{\text{vir, mock}} - M_{\text{BH, mock}}$ distribution with the example parameter set of $\gamma = 0$ and $\sigma_{\mu} = 0.3$.

Considering possible systematic uncertainties from 2D decomposition analysis, we added an error based on a normal distribution with 0.2 dex to the $M_{*,\text{mock}}$ to consider the uncertainty of $M*$ derived from observations, resulting in the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Procedure for generating mock galaxies at $z \sim 1.5$ with an example parameter set of $\alpha = \alpha_{\text{local}}$, $\beta = \beta_{\text{local}}$, $\gamma = 0$, and $\sigma_{\mu} = 0.3$. Orange and blue contours show the distribution of the entire sample (1σ–5σ) and that above the detection limit of the observed sample (1σ–3σ, see Sections 5.1 for the definition of the detection). The individual panels are as follows: (a) true black hole mass $M_{\text{BH,mock}}$ vs. true stellar mass $M_{\text{BH,mock}}$, before adding observational uncertainties, (b) 2–10 keV flux ($F_{[2–10\text{keV}],\text{mock}}$), (c) Hα FWHM (FWHM$_{H\alpha,\text{mock}}$), (d) pseudo-observed virial mass $M_{\text{vir,mock}}$ vs. true black hole mass $M_{\text{BH,mock}}$ where $M_{\text{vir,mock}}$ is the mock black hole mass after considering the observational bias. (e) Pseudo-observed stellar mass $\tilde{M}_{\text{star,mock}}$ vs. stellar mass $M_{\text{star,mock}}$ where $\tilde{M}_{\text{star,mock}}$ is the stellar mass after considering the observational bias, and (f) pseudo-observed virial mass $M_{\text{vir,mock}}$ vs. pseudo-observed stellar mass $\tilde{M}_{\text{star,mock}}$. The latter is to be compared with the observed $M_{\text{BH}}-M_{\star}$ plane shown in the left panel of Figure 12. Black dashed lines in panels (a) and (f) indicate $M_{\text{BH}}-M_{\star}$ local relation obtained from the fitting of the local galaxies (N. Häring & H.-W. Rix 2004; V. N. Bennert et al. 2011b; see Section 5), and gray dashed lines in panels (d) and (e) indicate $y = x$ line.</figcaption></figure>

mock observed stellar mass $\tilde{M}_{\star,\text{mock}}$. The resulted $\tilde{M}_{\star,\text{mock}} - M_{\star,\text{mock}}$ distribution is shown in Figure 6(e).

Finally, we get the mock data set of the redshift $z_{\text{mock}}$, observed stellar mass $\tilde{M}_{\star,\text{mock}}$, and virial BH mass $M_{\text{vir,mock}}$. Figure 6(f) shows the final mock observed mass distribution $(M_{\text{vir,mock}} - \tilde{M}_{\star,\text{mock}})$ with the example parameter set of $\gamma = 0$ and $\sigma_{\mu} = 0.3$. In Section 5, we discuss the mass relation, comparing our results with this mock catalog.

# 4. Results of AGN–Host Decomposition

## 4.1. Morphological Parameters

First, we examine the results obtained from the 2D decomposition. In Figures 7(a)–(c), we show the distribution of $n$ and $r_e$ (in the unit of pixel and kiloparsecs) for each filter using the top-5 stacked PSF. The distribution of $n$ has a peak in the distribution at $n \sim 1–2$, clearly seen in the F277W and F444W filters. This is similar to studies of AGN hosts at high redshift, which show hosts characterized by disklike morphology.

In Figure 7(d), we compare $H/T$ with $S/N_{\text{host}}$. We can see a strong correlation between the reconstructed $H/T$ and $S/N_{\text{host}}$. High $H/T$ means host galaxies dominate the AGN + host galaxy composite images, and we can easily detect host galaxies with high $S/N_{\text{host}}$; thus, our results indicate the validity of our analyses for the majority of the sample. Estimated morphological parameters and basic information for each host galaxy are reported in Table 3. In Appendix C, we

confirm that galight can reconstruct $r_e$ and $H/T$ correctly by running galight on mock galaxy images.

## 4.2. Impact of Different PSF Models

In Figures 8(a)–(c), we compare morphological parameters, $n$, $r_e$, and the host-to-total flux ratio $H/T$, in F444W obtained using each final PSFs. These parameters fall mostly along the $y = x$ line, indicating a strong correlation. Thus, we can say that different PSF reconstruction methods do not significantly affect the results.

In Figure 8(a), the size ($r_e$) comparison shows more consistent results than $n$, indicating that it is minimally affected by different PSFs. On the other hand, in comparing $n$ (Figure 8(b)), regardless of the PSF, we can find a consistent estimation on the low $n$ side ($n \lesssim 3$). However, at larger $n$ ($\gtrsim 4$), the scatter increases. This tendency is likely because larger $n$ implies a compact Sérsic profile that resembles the PSF, making it challenging to distinguish from the PSF. We also find a tendency that PSFEx PSFs result in slightly larger $n$ than the top-5 or top-75% stacked PSFs. Even when considering the other filters, we find highly consistent results, with a trend of increased scatter at higher $n$. However, we note that the individual values of $n$ have considerable uncertainty as shown in Appendix C.

Regarding $H/T$, we primarily see consistent results with strong correlations in Figure 8(c). While $H/T$ estimated using the top-5 and top-75% PSFs is largely consistent with each other, the PSFEx PSFs tend to result in slightly higher $H/T$.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Summary of the decomposition results in each filter. Columns (a), (b), and (c) respectively show the distribution of the estimated *n*, log *r*<sub>e</sub> in the unit of pixel and kiloparsecs. In column (d), the relation between *H/T* and *S/N*<sub>host</sub> is shown. In each panel, red (gray) indicates detection (nondetection) in each filter.</figcaption></figure>

than the other two PSFs. M.-Y. Zhuang & Y. Shen (2023)
suggested that using a narrower PSF than an exact PSF could
overestimate the host flux and *n*. Thus, the above biases in *n*
and *H/T* could be explained by the fact that PSFEx PSFs for
F444W have a slightly narrower PSF than the other PSFs, as
shown in Figure 4. A detailed comparison of the estimated *H/T*
in other filters is summarized in Appendix B.

In conclusion, the decomposition results using different
PSFs are generally highly consistent with each other. A
comprehensive discussion of the technical and practical
differences between PSF reconstruction methods will also be
provided in Section 6.3. Nonetheless, these discussions are
based on the comparisons between estimated values, and here,
we cannot definitively determine the true exact value. Related
to this, we provide the result of mock tests using different final
PSFs in Appendix C. Also, note that all final PSFs are not
single stellar images but stacked or modeled PSFs based on
multiple stellar images. As mentioned in Section 3.1.3, single
PSFs often have a wide range of sizes and shapes. Thus, using
a single stellar image without testing other stellar images can
risk misinterpreting the PSF image and giving different results.

4.3. *Host Images*

As seen in these high-quality host galaxy images, 2D decomposition analyses of JWST images open up the potential for a more detailed image-based galaxy analysis, such as studies traditionally conducted on inactive galaxies. Figure 9 shows three-color images (F277W, F150W, and F115W for red, green, and blue, hereafter RGB) of each host galaxy created by subtracting the PSF and the nearby Sérsic component.

First, thanks to the high spatial resolution, deep observations,
and meticulous decomposition analysis, we can access the
highest-quality AGN-host galaxy images up to z ~ 2.5,
allowing us to identify substructures. Particularly, in the case of
CID-273 (z = 1.85) and CID-307 (z = 2.05), despite their
redshifts ~2, we can clearly identify blue spiral arms with an
overall diffuse red broad component.

Additionally, galaxies such as CID-54 (z = 0.97), CID-510
(z = 1.12), CID-361 (z = 1.18), CID-445 (z = 1.26), and CID-
452 (z = 1.41) show more reddish colors at their centers than in
the outer region, indicating the possibility of having a bulge-like
structure or highly dust-obscured region (e.g., K. Ito et al. 2024).
<table><caption>**Table 3**
Data Summary of Our Sample<sup>a</sup></caption><thead><tr><th>CID<sup>b</sup></th><th>XID<sup>c</sup></th><th>z</th><th>R.A.<sup>d</sup><br>(deg)</th><th>Decl.<sup>d</sup><br>(deg)</th><th>log M<sub>BH</sub><sup>e</sup><br>(M<sub>☉</sub>)</th><th>line<sup>f</sup></th><th>log M<sub>*</sub><sup>g</sup><br>(M<sub>☉</sub>)</th><th>r<sup>h</sup><sub>e,F115W</sub><br>(arcsec)</th><th>n<sub>F115W</sub><sup>h</sup></th><th>log H/T<sub>F115W</sub><sup>h</sup></th><th>b/a<sub>F115W</sub><sup>h</sup></th><th>PA<sub>F115W</sub><sup>h,i</sup><br>(deg)</th><th>BIC<sub>F115W</sub><sup>h</sup></th></tr></thead><tbody><tr><td>1:157</td><td>5404</td><td>1.33</td><td>149.6751</td><td>1.9828</td><td>8.73</td><td>2</td><td>10.43 ± 0.24</td><td>0.27</td><td>5.2</td><td>-0.38</td><td>0.76</td><td>48.43</td><td>10,197.79</td></tr><tr><td>2:162</td><td>169</td><td>2.46</td><td>149.7359</td><td>2.0276</td><td>9.30</td><td>1</td><td>11.28 ± 0.17</td><td>0.06</td><td>7.0</td><td>-0.57</td><td>0.33</td><td>72.45</td><td>4176.85</td></tr><tr><td>4:416</td><td>404</td><td>1.80</td><td>149.8129</td><td>2.3455</td><td>8.41</td><td>3</td><td>10.98 ± 0.13</td><td>0.22</td><td>1.9</td><td>-0.36</td><td>0.78</td><td>-31.62</td><td>8593.73</td></tr><tr><td>5:203</td><td>341</td><td>1.36</td><td>149.8142</td><td>2.0164</td><td>8.34</td><td>3</td><td>10.18 ± 0.25</td><td>0.18</td><td>1.5</td><td>-0.37</td><td>0.42</td><td>-54.74</td><td>3194.34</td></tr><tr><td>6:307</td><td>285</td><td>2.05</td><td>149.8227</td><td>2.0897</td><td>8.93</td><td>3</td><td>11.08 ± 0.21</td><td>1.5</td><td>7.0</td><td>-0.39</td><td>0.63</td><td>51.64</td><td>6704.09</td></tr></tbody></table>

**Notes.**

<sup>a</sup> The full version of this catalog is to be available, which will include results based on all bands.

<sup>b</sup> ID in Chandra-COSMOS Survey (M. Elvis et al. 2009; F. Civano et al. 2012)

<sup>c</sup> ID in XMM-COSMOS (N. Cappelluti et al. 2009; M. Brusa et al. 2010)

<sup>d</sup> J2000.

<sup>e</sup> $M_{BH}$ is estimated with single-epoch virial mass estimation using each line. We assume the 0.4 dex uncertainty in this study.

<sup>f</sup> The emission line used to estimate $M_{BH}$. 1:Hβ (FMOS-COSMOS), 2:Hα (FMOS-COSMOS), 3:Mg II (zCOSMOS-bright), 4:Mg II (zCOSMOS-deep).

<sup>g</sup> $M_{*}$ is estimated using CIGALE with the top-5 PSF results.

<sup>h</sup> Table 3 only shows morphological parameters in F115W measured with the top-5 PSF. Parameters in the other filters are to be available in the full version.

<sup>i</sup> Position angle measured from north to west.

(This table is available in its entirety in machine-readable form in the online article.)

Furthermore, there are cases with extended red structures (e.g., CID-668; $z = 0.97$), which may indicate the presence of dust lanes as seen in the X-ray obscured (type-II) AGNs (J. D. Silverman et al. 2023).

## 4.4. Estimated $M_{*}$ and the Comparison of Different SED Fitting Methods

As described in Section 3.5, we perform SED fitting of the host galaxy photometry to estimate $M_{*}$. For five galaxies detected in less than two bands, we estimate $M_{*}$ upper limits by assuming $H/T = 0.2$ and the average F277W photometry-to-$M_{*}$ ratio in the sample. The assumed $H/T$ of 0.2 is because galaxies with $H/T \gtrsim 0.2$ are almost detected in our methods (Figure 7(d)). The inferred $M_{*}$ for each host galaxy is reported in Table 3. Figure 10 shows four examples of the SED fitting results with the residuals. Generally, for the objects with $z \gtrsim 1$, the photometry or upper limit from F814W and F115W fall at a rest-frame wavelength shorter than 4000 Å break and are important in constraining stellar population parameters of host galaxies.

We independently employ two distinct SED fitting codes, CIGALE and Prospector, as explained in Section 3.5. Both codes are run having as similar parameter settings as possible. In Figure 11, we compare the parameters obtained from both codes. As shown in Figure 11(a), the results exhibit a significantly high positive correlation. However, we find an offset of approximately $\Delta \log M_{*} = +0.13$ dex (corresponding to $1\sigma$). This offset is not far from common systematic $M_{*}$ uncertainties among SED fitting methods reported in C. Pacifici et al. (2023). It remains challenging to determine whether $M_{*}$ from CIGALE or Prospector is more accurate. In this study, we primarily use the CIGALE $M_{*}$ in the main discussion to maintain consistency with previous studies for AGN-host galaxies (e.g., F. Zou et al. 2019; L. Shen et al. 2020; T. Ishino et al. 2020; J. Li et al. 2021; L. Koutoulidis et al. 2022; M.-Y. Zhuang et al. 2023; J. Li et al. 2024).

Finally, in Figure 8(d), we compare the estimated $M_{*}$ with each final PSF. The estimated $M_{*}$ are on the 1:1 line and show a strong correlation, suggesting that the estimated $M_{*}$ with different final PSFs remain highly consistent. This consistency

can be attributed to the fact that, as shown in Figure 8(d), $H/T$ or host flux from different final PSFs is also consistent.

We compare our CIGALE $M_{*}$ with M.-Y. Zhuang et al. (2023), which also utilized CIGALE on the COSMOS-Web data, and find that these two measurements well agree with each other with a scatter of $\Delta \log M_{*} = +0.08_{-0.18}^{+0.20}$. This is very encouraging since their decomposition analysis is independent of our effort.

In Figures 11(b)–(d), we compare results for the other output SED model parameters: $A_{V}$, $t_{age}$, and $\tau$. While $A_{V}$ and $t_{age}$ exhibit large uncertainties, they show a positive correlation, with the median offset being close to zero within the range of uncertainties. Regarding $\tau$, it is evident that significantly inconsistent values are observed around $\log(\tau/Gyr) \sim 0.5$. Considering that $t_{age}$ is generally $\log(t_{age}/Gyr) \lesssim 0.6$ in our sample, this discrepancy can be attributed to the challenge of accurately determining SFH when $\tau \gg t_{age}$.

Additionally, we confirmed a strong negative correlation ($\rho = -0.60, p \ll 0.05$) between $\Delta \log \tau$ and $\Delta \log M_{*}$. This relation is because larger $\tau$ indicates a longer-lasting SFH, i.e., star formation has been persisting more recently. Consequently, galaxies with larger $\tau$ tend to host more young stellar populations, resulting in a smaller mass-to-light ratio and smaller $M_{*}$. While these differences may be attributed to differences in stellar models or fitting strategies (i.e., Bayesian or non-Bayesian), further investigation is omitted in this paper since it does not impact the main results of this study.

## 5. $M_{BH}-M_{*}$ Relation

In the left panel of Figure 12, we plot our measurements of $M_{BH}$ as a function of $M_{*}$. Based on the large sample covering a broad range in both parameters, we find a weak positive correlation with a weighted Spearman's correlation coefficient<sup>34</sup> of $\rho = 0.25$. As done in previous local studies, we attempt to fit the observational data with a linear function of

$$ \log M_{BH} = \alpha \log M_{*} + \beta. \quad (10) $$

<footnote n="34"><sup>34</sup> The weighted Spearman's correlation function is calculated with R/CRAN wCorr package (P. Bailey & A. Emad 2023).</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 8. Comparison of the decomposition results ($r_e$ in the unit of arcseconds, $n$, and $H/T$) based on the F444W filter using different PSFs (columns (a)–(c)). Data points in red have satisfied our stringent criteria for host detection (Section 3.3). Column (d) compares $M_*$ estimated by the SED fitting, with colors corresponding to whether each object was detected in more than two bands. Spearman's correlation coefficient $\rho$ for each prior is shown in the lower right corner of each panel. High correlation coefficients and the distribution around $y = x$ (black dashed line) suggest consistent results among the fitting with different PSFs.</figcaption></figure>

The red line represents the results with $\alpha \simeq 0.42_{-0.08}^{+0.12}$ and $\beta \simeq 4.08_{-1.26}^{+0.86}$, and the orange-shaded region shows the $1\sigma$ confidence interval. The coefficient of determination $R^2$ for this linear fitting is 0.049, which indicates that this linear model does not sufficiently explain the data. We test the validity of this linear model by fitting the data with a model with the fixed $\alpha$ of 0 and compared the BIC values. Although the BIC for the $\alpha$-free model is lower than the BIC for the $\alpha = 0$ model, the difference in BIC is less than 10. This means the $\alpha$-free model is not significantly better than the $\alpha = 0$ model. We also run the Shapiro-Wilk test and got a $p$-value of 0.37. Therefore, at a 5% significance level, the null hypothesis cannot be rejected, and we cannot say that the observed data do not follow a Gaussian distribution. However, it is important to note that this result does not consider the selection bias and does not necessarily represent the intrinsic distribution on the $M_{BH}-M_*$ plane at high redshift.

For comparison, we also provide the linear fit to the local sample consisting of 30 inactive galaxies (N. Häring & H.-W. Rix 2004) and 25 active galaxies (V. N. Bennert et al.

2011b). The fit to these 55 galaxies results in $\alpha_{local} \simeq 0.97_{-0.11}^{+0.10}$ and $\beta_{local} \simeq -2.48_{-1.11}^{+1.15}$. $R^2$ for this fitting is 0.63, and this suggests that the local $M_{BH}-M_*$ relation is well described by the linear model. Note that the local sample and our high-$z$ sample have different selection effects. Thus, we cannot directly compare $\alpha$ and $\beta$ with $\alpha_{local}$ and $\beta_{local}$ (see Sections 5.1 and 5.2 for the discussion with consideration of selection bias).

For investigating the redshift dependence of the mass relation, we calculate $\Delta \log (M_{BH}/M_{\odot})$, the relative offset of the black hole mass at given $M_*$ from the local relation:

$$ \Delta \log (M_{BH}/M_{\odot}) = M_{BH} - \alpha_{local} \log (M_{*}/M_{\odot}) - \beta_{local}. \quad (11) $$

For $\alpha_{local}$ and $\beta_{local}$, we use the above values from the local samples (N. Häring & H.-W. Rix 2004; V. N. Bennert et al. 2011b), i.e., 0.97 and $-2.48$, respectively. We plot $\Delta \log (M_{BH}/M_{\odot})$ as a function of $z$ in the right panel of Figure 12. We then parameterize $\Delta \log (M_{BH}/M_{\odot})$ to evolve with $z$ as

$$ \Delta \log (M_{BH}/M_{\odot}) = \gamma \log (1 + z). \quad (12) $$
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Three-color cutout PSF-subtracted images (the F277W, F150W, and F115W for RGB) of all host galaxies detected in all F115W, F150W, and F277W in the order of redshift. We have not performed deconvolutions based on the PSF FWHM differences between each filter, and consequently, each filter image possesses a different-sized PSF. The target IDs and the redshifts are shown in the up left corner of each image. The white bars in the lower left part are 1" in length. Thanks to the high spatial resolved deep observation of JWST and the careful decomposition analysis, we can clearly identify substructures, such as bulges, spirals, bars, and dust lanes.</figcaption></figure>

Here, we assume there is no redshift change in $\alpha$ and $\beta$, and $\Delta \log M_{BH}$ can be described by only the evolution from the local relation ($\log M_{BH} = \alpha_{\text{local}} \log M_{*} + \beta_{\text{local}}$). Fitting our data with Equation (12) without considering the selection bias

results in $\gamma = 1.33_{-0.14}^{+0.13}$, suggesting positive evolution. Although $R^2$ for this fitting is 0.050, the BIC for this fitting is significantly lower than the $\gamma = 0$ model (i.e., $\Delta \log (M_{BH}/M_{\odot}) = 0$ with a scatter); thus, this model describes
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Representative examples of SED fits to PSF-subtracted host galaxies. Source ID and their estimated $M_*$ with $1\sigma$ confidence range are shown in the top right corner of each plot. The inverted triangle indicates the $3\sigma$ upper limit due to the nondetection in each filter. The bottom colored shades show the transmission curves of F814W, F115W, F150W, F277W, and F444W from left to right. The lower plot in each panel shows the difference between the observational and best-fit model photometry scaled by the error of the data.</figcaption></figure>

the data better than the $\gamma = 0$ model. However, because our sample is (X-ray) flux limited, it is essential to consider the impact of selection bias when determining the mass relation (e.g., X. Ding et al. 2020b; J. Li et al. 2021, J. Li et al. 2022). Thus, for the rest of this section, we measure the intrinsic slope of the mass relation ($\alpha$) at $z \sim 1.5$ and the redshift evolution parameter ($\gamma$) by comparing our results with mock catalogs as described in Section 3.6. This approach allows us to account for selection bias and measurement uncertainties to determine the intrinsic redshift evolution and dispersion of the mass relation.

### 5.1. Intrinsic Slope ($\alpha$) of the Mass Relation at $z \sim 1.5$

Past efforts to establish the evolution of the ratio between black hole and galaxy mass assume that there is a linear relation at higher redshifts, and the slope of the relation matches the local relation (e.g., X. Ding et al. 2020b; J. Li et al. 2021). This may not necessarily be the case. Here, we assess how well the parameters of a linear relation can be constrained.

In particular, we initially assumed that $\alpha$ and $\beta$, for constructing the mock samples (Equation (9)), have constant values independent of the redshift, and the evolution is expressed simply as $\gamma \log(1+z)$. However, the results in the previous Section 5 may indicate a different $\alpha$ from the local value without considering selection biases. In this section, we examine whether the constant $\alpha$ assumption is valid by estimating the intrinsic value of $\alpha$ while considering selection bias and provide evidence for an intrinsic relation between $M_{BH}$ and $M_*$ at $z > 1$ for the first time.

Since the parameters $\alpha$, $\beta$, and $\gamma$ exhibit degeneracies, in this subsection, we constrain the $z$ range to $z = 1–2$ and perform fitting for $\alpha$ and $\beta$ with an assumption of $\gamma = 0$ over this redshift range. Consequently, we cannot compare the estimated $\beta$ with the values in the local relation directly. For the intrinsic dispersion of the mass relation $\sigma_\mu$, we assume $0.3$ dex in this subsection, thus matching local studies.

With the mock observed data described in Section 3.6, we apply the selection criteria corresponding to our sample (Section 2.3) and compare them with the observed results to
<figure><img src="image_13.png" /><figcaption>Figure 11. Comparison of SED fitting results with two different SED fitting codes. Panels (a) to (d) compare the estimated log M*, A*, tage, and log τ. X- and Y-axes correspond to the result with CIGALE and Prospector. The black dashed line indicates the x = y line. The median and 1σ confidence level of the difference between the two results (defined as Δ = Prospector − CIGALE) and Spearman's correlation coefficient are shown in the lower right part in each panel.</figcaption></figure>

constrain the evolution parameters. We assume the observation thresholds as

$$F_{[2-10 \text{ keV}]} > F_{[2-10 \text{ keV}],\text{lim}} \qquad (13)$$

$$FWHM_{\text{line}} > 1000 \text{ km s}^{-1}, \qquad (14)$$

$$M_{*} > 10^{10} M_{\odot} \qquad (15)$$

where the first and second conditions correspond to the detection limit of the X-ray observation and broad lines for single-epoch M<sub>BH</sub> estimation. The third condition accounts for the detection limit of the host galaxies, and we confirm that galaxies with masses greater than 10<sup>10</sup>M<sub>☉</sub> are successfully detected across all redshift ranges in our decomposition analysis. Thus, in comparing mock and real observations, galaxies with masses below 10<sup>10</sup>M<sub>☉</sub> are excluded to maintain consistency with this mock observation. As shown in Figure 2(a), excluding the eight objects, all sources in our sample have F<sub>[2-10keV]</sub> above the XMM-Newton F<sub>[2-10keV],lim</sub>. Note that all of the eight targets with F<sub>[2-10keV]</sub> smaller than XMM-Newton F<sub>[2-10 keV],lim</sub> are targets observed only with Chandra. We changed F<sub>[2-10 keV],lim</sub> depending on in which survey each target was detected. Because M<sub>BH</sub> of our real targets are based on single-epoch estimation with three different lines (Hα, Hβ, and MgII), we have three different selection thresholds on the mock galaxy using FWHM<sub>Hα</sub>, FWHM<sub>Hβ</sub>, and FWHM<sub>Mg ii</sub>.

For each of our targets, we first select the mock galaxies with the corresponding selection bias, i.e., the selection condition using the same line information used to estimate M<sub>BH</sub>. Furthermore, we select the mock galaxies with a similar redshift; |z − z<sub>mock</sub>| < 0.1. Then, we calculate the probability that the mock galaxies with the similar z<sub>mock</sub> as the real galaxy would have the same M&#770;<sub>*,mock</sub> with |ΔM*| < 0.1 and M<sub>vir,mock</sub> with |ΔM<sub>BH</sub>| < 0.1, i.e., calculate the probability p following

$$p = \frac{N_{[\Delta z|<0.1, \quad |\Delta M_{*}|<0.1, \quad |\Delta M_{BH}|<0.1]}}{N_{[\Delta z|<0.1]}} \qquad (16)$$

We then calculate the likelihood of our sample being observed for each parameter combination of α and β. Thereby, we estimate the probability distribution of these parameters using MCMC. In the sampling, we assume a uniform prior between 0.1 to 3 for α, and -25 to 8 for β.

The obtained α − β distribution is shown on the left panel of Figure 13, indicating a strong anticorrelated degeneracy between α and β. Because the 1σ contour includes values

(α<sub>local</sub>, β<sub>local</sub>) corresponding to the local relation, our results do not definitively reject the scenario where α and β do not evolve compared to the local relation up to z ~ 2.5. In the right panels of Figure 13, we show the M<sub>BH</sub> − M* distribution for mock observed galaxies (similar to Figure 6(f)) generated with manually sampled parameters on the ridge of the α − β degeneracy and the real observed galaxies (α ~ 0.2–1.8). The plots suggest that the mock data exhibit a similar distribution to our sample within the comparable mass range. To break this degeneracy and improve the precision of determining α and β, a larger sample with a wider M<sub>BH</sub> range in future studies is essential. Even so, we demonstrate that a relation between M<sub>BH</sub> and M* at high-z is realized based on having a statistical sample afforded by the COSMOS-Web data set. As indicated by Figure 8(d), the difference in the PSF reconstruction methods does not significantly affect the M* estimation. Therefore, the posterior distribution of α and β that is estimated based on M* also shows no significant PSF dependency.

## 5.2. Intrinsic Evolution (γ) of the M<sub>BH</sub>/M* Relation

To determine the evolution of mass relation with consideration of the selection bias and measurement uncertainties, we generate mock observed catalogs with free parameters of γ (evolution rate) and σ<sub>μ</sub> (intrinsic dispersion of the mass relation) and constrain them by comparing the mock catalogs with observational data in a similar manner to Section 5.1. Note that the assumption of γ = 0 and σ<sub>μ</sub> = 0.3 used in section 5.1 is not applied in the fitting performed in this subsection; both parameters are treated as free parameters during the fitting process here. In contrast to Section 5.1, we fixed α and β in Equation (9) to the values in the local relation (α<sub>local</sub>, β<sub>local</sub>) to consider the redshift evolution. As illustrated in Section 5.1, we cannot rule out the possibility for evolution of α. However, as described in Section 5.1, α, β, and γ are strongly degenerate, and obtaining physically meaningful results is challenging when all three parameters are left free for the sample being considered here.

There is still a possibility that σ<sub>μ</sub> depends on redshift. Nevertheless, as discussed later, imposing strong constraints on σ<sub>μ</sub> in our results is challenging due to the sample size and its uncertainties. Therefore, we set σ<sub>μ</sub> as a constant independent of redshift in the fitting. It means σ<sub>μ</sub> obtained through this method is considered to be an averaged value over z ~ 0.68–2.5. Even so, this σ<sub>μ</sub> estimation has the highest statistical significance for such a study at z ≳ 1. We discuss the redshift evolution of σ<sub>μ</sub> in
<figure><img src="image_14.png" /><figcaption>**Figure 12.** (Left) Observed $M_{BH}$ and $M_*$ distribution as shown by the red circles. A linear fit is indicated by the red line with the shaded region indicating the $1\sigma$ confidence range. The best-fit relation and Spearman's correlation coefficient are reported at the top of the panel. (Right) $\Delta \log (M_{BH}/M_{\odot})$ as a function of redshift where zero corresponds to the local relation as marked by the gray band having a width representing the local dispersion. In both panels, samples from previous studies using HST are plotted: gray for low-$z$ (N. Häring & H.-W. Rix 2004; V. N. Bennert et al. 2011b), green for intermediate-$z$ (K. Jahnke et al. 2009; V. N. Bennert et al. 2011a; M. Cisternas et al. 2011; M. Schramm & J. D. Silverman 2013), and orange indicating those from X. Ding et al. (2020b) at $z \sim 1.5$. Red triangles indicate the $M_*$ upper limit and $\Delta \log (M_{BH}/M_{\odot})$ lower limit for five undetected targets.</figcaption></figure>

Section 6.2. Finally, in this analysis, we assume no redshift-dependent parameters among the free parameters. Therefore, there is no need to restrict the redshift range within the data, as in Section 5.1.

In the fitting process, we assume a uniform prior distribution for $\gamma$ between $-1 < \gamma < 1$. Then, we have three different prior settings for $\sigma_{\mu}$: a uniform distribution between $0.01 < \sigma_{\mu} < 1.0$ (basic), a Gaussian distribution with a mean of 0.3 dex and a standard deviation of 0.1 dex (Gaussian), and a uniform distribution between $0.25 < \sigma_{\mu} < 1.0$ with a prohibition of $\sigma_{\mu} < 0.25$ (realistic).

Figure 14 shows the estimated posterior distributions of $\gamma$ and $\sigma_{\mu}$ using each prior setting with the top-5 stacked PSF results. As evident in all panels by the orange contours, the intrinsic dispersion $\sigma_{\mu}$ is strongly degenerate with the evolution rate $\gamma$ where a smaller $\gamma$ results in a larger $\sigma_{\mu}$ as demonstrated in previous studies (X. Ding et al. 2020b; J. Li et al. 2021). To reiterate, a smaller value for $\gamma$ biases the mass relation toward relatively lower $M_{BH}$ values; thus, a larger $\sigma_{\mu}$ is required to reproduce a certain set of observation data.

Considering the likelihood distribution for $\gamma$, the “basic” prior setting shows a slight positive-to-no evolution with $\gamma = 0.46_{-0.63}^{+0.52}$. Similarly, the “Gaussian” prior setting results in $\gamma = 0.48_{-0.62}^{+0.31}$. If we assume that the intrinsic dispersion should not be significantly smaller than the local dispersion, we limit the allowed range for $\sigma_{\mu}$ to be above 0.25 ("Realistic" case). In this case, we find $\gamma = 0.22_{-0.58}^{+0.39}$, closer to the case for no evolution ($\gamma = 0$) than the results with the other priors. For the latter, the intrinsic dispersion is slightly higher at $0.38_{-0.09}^{+0.12}$. In all cases, our results are consistent with very mild or essentially a lack of evolution with respect to the local relation.

Then, the left panel of Figure 15 compares mock observations using the median parameters ($\gamma = 0.48$ and $\sigma_{\mu} = 0.30$) under the assumption of “Gaussian” prior with the actual observational $M_{BH}-M_*$ distribution and relation. We can see

that the mock data can explain the observed data well. The right panel of Figure 15 compares the intrinsic relationship, i.e., the relation corrected for the selection bias, with the local relation. Again, the resulting intrinsic relation is consistent with the local relation within the range of errors.

## 6. Discussion

In Section 5.2, our findings suggest a mild or lack of evolution of the mass relation from the local relation when considering selection biases and measurement uncertainties. In this section, we first compare the derived values of $\gamma$ and $\sigma_{\mu}$ from Section 5.2 with previous studies. Then, we also discuss the cosmic averaging scenario (C. Y. Peng 2007; M. Hirschmann et al. 2010; K. Jahnke & A. V. Macciò 2011).

### 6.1. Comparison to Other Studies

First, the conclusion of no or mild evolution with $\gamma = 0.48_{-0.62}^{+0.31}$ is consistent with studies based on 2D decomposition analysis (e.g., X. Ding et al. 2020b; J. Li et al. 2021) and studies using an SED-fitting-based decomposition method (e.g., M. Sun et al. 2015; H. Suh et al. 2020).

In Figure 14, we compare our results on the estimated $\gamma - \sigma_{\mu}$ distribution to those of J. Li et al. (2021). Our sample has higher redshift range than J. Li et al. (2021), and the change of $\Delta \log (M_{BH}/M_{\odot})$ is also proportional to $\log(1+z)$ in our model. Therefore, to reproduce the observational results, when increasing $\gamma$, we need to decrease $\sigma_{\mu}$ more than that in J. Li et al. (2021). In other words, the slope of the $\gamma-\sigma_{\mu}$ degenerate relation is steeper in our study. As a result, while the sample size of this study is approximately 6 times smaller than that in J. Li et al. (2021), the uncertainty of $\gamma$ is only ~2 times larger than that in J. Li et al. (2021). On the other hand, due to the steep slope, imposing constraints on $\sigma_{\mu}$ becomes challenging, and the uncertainty becomes ~4 times larger than that in J. Li et al. (2021). Even so, our estimated value of $\sigma_{\mu}$ is $0.30_{-0.13}^{+0.14}$,
<figure><img src="image_15.png" /><figcaption>**Figure 13.** (Left) Posterior distribution in the $\alpha - \beta$ plane obtained through a comparison between mock observation and our sample. The posterior distributions for $\alpha$ and $\beta$ are shown in the right and upper 1D histogram. Each contour indicates $1\sigma$, $2\sigma$, and $3\sigma$ levels. The median and $1\sigma$ confidence levels of $\alpha$ and $\beta$ are shown in the top right corner of each panel. (Right) Comparison of mock observed masses generated with each $\alpha$ and $\beta$ at $z = 1.5$ with the real observed data at $z = 1-2$ shown in red circles. Contours show the distribution of mock observation, indicating $1\sigma$, $2\sigma$, and $3\sigma$ levels. Each parameter set, $\alpha$ and $\beta$, is shown in the upper left corner of each panel. As shown, we cannot reject the scenario that $\alpha$ does not evolve up to $z \sim 2.5$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Posterior distribution on the $\sigma_{\mu} - \gamma$ plane obtained through a comparison between mock observation and our sample (the result with the top-5 stacked PSFs). Panels (a) to (c) correspond to the result using basic, Gaussian, and realistic prior settings, respectively. The posterior distribution for $\sigma_{\mu}$ and $\gamma$ are shown in the right and upper 1D histogram. Each contour (red to orange) indicates $1\sigma$, $2\sigma$, and $3\sigma$ from inner to outer. The median and $1\sigma$ confidence levels of $\gamma$ and $\sigma_{\mu}$ are shown in the top right corner of each panel. Only in panel (b), we plot the contour from J. Li et al. (2021) in gray, and their $\gamma$ and $\sigma_{\mu}$ estimations are shown in the lower left corner based on SDSS quasars at $0.3 < z < 0.8$ with Subaru HSC imaging. Regardless of the prior, the results suggest no or mild evolution at $z < 2.5$.</figcaption></figure>

which is remarkably similar to that in J. Li et al. (2021) with $\sigma_{\mu}=0.25_{-0.04}^{+0.03}$.

It is worth highlighting that the inference on the value of $\gamma$ is very close to zero for the “Realistic” case (Figure 14(c)) where we assume that the intrinsic scatter ($\sigma_{\mu}$) cannot be lower than the local dispersion. Interestingly, if $\sigma_{\mu}$ ($\sim 0.4-0.5$) is actually higher than the local value, this would push the evolution parameter to negative values ($\gamma \sim -0.5$), thus presenting a scenario where the black holes have to catch up to their host galaxies by a bit.

## 6.2. Scatter ($\sigma_{\mu}$) Evolution and Cosmic Averaging

When assuming a noncausal cosmic averaging scenario (K. Jahnke & A. V. Macciò 2011), major mergers average and

equalize the mass ratio $M_{BH}/M_*$ through cosmic history. Thus, $\sigma_{\mu}$ should increase toward high redshift.

To test the cosmic averaging scenario with our data, we generate a mock sample (20,000 parameter sets of $M_{BH}$ and $M_*$) at $z \sim 1.4$, the median redshift of our observational sample. We generate $M_*$ based on the SMF by J. R. Weaver et al. (2023), and calculate the $M_{BH}$ assuming the $\gamma$ and $\sigma$ sampled in the MCMC run with “Gaussian” prior setting (Section 5.2). Each mock galaxy is assumed to undergo major mergers following the major merger rate from V. Rodriguez-Gomez et al. (2015) covering $z \sim 0-1.4$. We simulate the redshift evolution of $M_{BH}$ and $M_*$ by summing them with those of merging partners. Then, we calculate $\sigma_{\mu}$ at each redshift to trace the expected scatter evolution from the assumed conditions. In this simulation, we do not consider the accretion
<figure><img src="image_17.png" /><figcaption>**Figure 15.** (Left) Comparison of the observed sample and the mock observation of $z \sim 1.5$ AGN on the $M_{BH}-M_{*}$ plane. The mock data at $z \sim 1.5$ is constructed with the best-fit $\gamma$ and $\sigma_{\mu}$ (the result with the top-5 stacked PSFs assuming the “Gaussian” prior) while assuming $\alpha_{local}$ and $\beta_{local}$. Small red circles with error bars and the red dashed lines indicate the observed data and the linear fit results (same with Figure 12 left). Contours indicate a distribution of mock observed samples, showing $1\sigma$, $2\sigma$, and $3\sigma$ from inner to outer. (Right) Comparison of the intrinsic $M_{BH}-M_{*}$ relation at $z \sim 1.5$ (the red line) to the local relation (the gray line).</figcaption></figure>

onto the black hole and star formation, i.e., both $M_{BH}$ and $M_{*}$ are assumed to grow only through major mergers. We limit the mergers to those with mass ratios within $\pm 0.5$ dex.

Figure 16 compares the simulated redshift evolution of $\sigma_{\mu}$ with the results from this study and J. Li et al. (2021). When assuming the evolution only through major mergers, the growth within the $1\sigma$ uncertainty range significantly encompasses the results of J. Li et al. (2021). Moreover, our median redshift evolution is consistent with the results of J. Li et al. (2021). Therefore, the $\sigma_{\mu}$ difference between the results from this study and J. Li et al. (2021) could be interpreted as the major merger-based scatter evolution.

However, due to the large uncertainty in our results, we cannot draw any definitive conclusions regarding the redshift evolution of $\sigma_{\mu}$. Furthermore, as evident from Figure 16, our sample has a wide redshift range compared to that of J. Li et al. (2021). If $\sigma_{\mu}$ varies with redshift, the sample should be binned in a narrower redshift range to trace the redshift evolution. Nevertheless, as mentioned in Sections 5.1 and 5.2, the degeneracy relation on the $\gamma - \sigma_{\mu}$ plane tends to steepen toward higher redshifts, making it relatively challenging to impose constraints on $\sigma_{\mu}$. Future high-$z$ statistical studies will likely require samples of a similar size compared to those of J. Li et al. (2021) with $N = 584$, or even larger to address the redshift evolution of $\sigma_{\mu}$. Thus, it will be necessary to conduct comprehensive surveys of high-$z$ AGNs using next-generation survey data such as Euclid and Roman.

### 6.3. General Notes on PSF Reconstruction Methods

So far, various studies have performed decomposition of JWST images (e.g., X. Ding et al. 2022a; X. Ding et al. 2023; Y. Harikane et al. 2023; M. A. Stone et al. 2023; M. Yue et al. 2023; M.-Y. Zhuang & Y. Shen 2023; M.-Y. Zhuang et al. 2023; M. A. Stone et al. 2023). As discussed above or in the previous studies, the results of the AGN+host galaxy 2D decomposition depend significantly on the PSF reconstruction.

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Evolution of $\sigma_{\mu}$ at $z \lesssim 2$ based on a simple simulation of the cosmic averaging scenario. The vertical red error bars represent the scatter in $\sigma_{\mu}$ within the $z$ range of our sample using the basic prior. Faint gray lines correspond to each sampled parameter set while the red line and the orange filled region represent the median and $1\sigma$ confidence level of $\sigma_{\mu}$ redshift evolution. The green data denotes the constraint from J. Li et al. (2021), which is consistent within the $1\sigma$ range of our results.</figcaption></figure>

In particular, M.-Y. Zhuang & Y. Shen (2023) discussed the effect of different PSF on decomposition results. M.-Y. Zhuang & Y. Shen (2023) compared three PSF modeling methods (SWarp, photutils, and PSFEx), but they did not compare them with $\chi^2$ based methods directly. Thus, in this paper, we summarize the comparison when using different PSF reconstruction methods.

In this study, we compare three final PSFs: two obtained by $\chi^2$-based methods (a top-5 stacked PSF and a top-75% stacked
PSF) and one from an empirically modeling method, PSFEx. As demonstrated in Figure 4, we find offsets in FWHM among different PSF reconstruction methods. These PSF variations, as indicated by M.-Y. Zhuang & Y. Shen (2023), could potentially introduce biases in estimating *n*, *r*<sub>e</sub>, or *H/T*. This is due to a broader (narrower) PSF than in reality, which tends to result in smaller (larger) *n*, larger (smaller) *r*<sub>e*,*c**, and overestimation of host fluxes. However, as shown in Figure 8, morphological parameters such as *n*, *r*<sub>e*,*c**, *H/T* generally exhibit consistent relations. Besides, as shown in Figure 8(d), different PSFs have less impact on *M*<sub>*</sub> estimation than on *H/T*. This could be due to the fact that *M*<sub>*</sub> is estimated from SED fitting (Section 3.5) using multiband photometry that averages the uncertainty in each band. If so, SED fitting with a smaller number of photometric bands (e.g., X. Ding et al. 2023; M. Yue et al. 2023) may lead to more severe effects from inaccurate PSF reconstruction. We also find each method has advantages and disadvantages from a technical aspect. Lastly, we summarize below the technical comparison between each method.

*Modeling method.* The approach of constructing an empirical PSF model from numerous stars, such as PSFEx in this study, has the advantage of being less influenced by noise compared to the χ²-based methods, as depicted in Figure 17. This method also allows flexible modeling, considering PSF as a function of position or brightness.

However, a drawback is the requirement for many PSF
candidates to model a local or flux-dependent PSF, which can
be considered a trade-off. Furthermore, we also confirm that
FWHM values of PSFEx PSFs depend on configuration
parameters. When reconstructing PSFs, it is challenging to
determine the best configuration parameters because the exact
PSFs are not known.

*χ²-based selection method.* Selecting PSFs based on χ²_ν
from a substantial number of stars, such as Top-5 or Top-75%
PSFs, allows easy analysis considering the PSF uncertainties.
Also, by fitting with various different single PSFs, the
possibility for PSF mismatches is minimized. Notably, M. Yue
et al. (2023) discussed the differences in broadband PSFs
attributed to variations in SED shapes between stars and AGNs.
Our method using all single PSFs in the PSF library may
consider this PSF uncertainty as a result of selecting low-χ²_ν
single PSFs with a matched shape.

However, using fewer PSF candidates, like the top-5 stacked
PSF, might increase the noise of the final PSF, as observed in
Figure 17. Note that this noise is generally smaller than the
central main component; thus, it should not affect significantly
except in the case with small *H/T*. Additionally, our approach
involves visual inspection in PSF candidate selection, which
might cause a bias. Also, as mentioned by M.-Y. Zhuang &
Y. Shen (2023), a lower χ²_ν does not necessarily mean more
correct PSF. Moreover, as a practical demerit, this method
needs more time to create a PSF library with visual inspection
in each filter and region and more computational cost for SED
fittings with all single PSFs.

Given that we do not know the correct answers in this study,
it is challenging to discuss which method produces the most
accurate results. For some targets, decomposition is clearly
successful with one final PSF, which can then be evaluated
from the residual emission based on other PSFs. We also
confirm that these failures in fitting can occur with all final
PSFs. Therefore, we conclude that it is best to assess the impact

on derived properties (e.g., *n*, *r*<sub>e</sub>, *H/T*) by varying the method of PSF construction (Section 3.1) and place equal weight on assessing the uncertainties based on varying PSFs (Section 3.3) to obtain solid 2D decomposition results.

7. Conclusions

We performed a 2D decomposition analysis of high-z
(z ~ 0.68–2.5) type-I AGNs using the COSMOS-Web
(C. M. Casey et al. 2023) and PRIMER-COSMOS surveys to
measure the black hole–stellar mass relation at high-z. Our
sample contains 107 targets that are X-ray-selected, broad-line
AGNs with single-epoch black hole mass estimates
(log (M<sub>BH</sub>/M<sub>☉</sub>) ~ 6.9 – 9.6) based on Hβ, Hα, and MgII
from previous spectroscopic surveys (e.g., A. Schulze et al.
2018).

By utilizing HST/ACS + JWST/NIRCam imaging that
covers the rest-frame optical to near-infrared, we obtained
multiband information of the AGN host galaxies with
unprecedented spatial resolution in which we can clearly
identify substructures such as bars and spirals arms (Figure 9).
Since AGN–host 2D decomposition is known to be sensitive to
the PSF reconstruction methods, we compared the results with
three final PSFs reconstructed using the modeling method
PSFEx and a χ²-based selection method. Through a meticulous
decomposition analysis using various PSFs, we successfully
detected the host galaxies in more than two filters for over 90%
of the entire sample. Then, we confirmed that host morpho-
logical parameters such as *n*, *r*<sub>e</sub>, and *H/T* remain relatively
consistent regardless of the PSF reconstruction method used
(Figure 8). Furthermore, given the high quality of the host
galaxy images, this study is expected to serve as a crucial
stepping stone for image-based spatially resolved investiga-
tions of AGN host galaxies, such as double-Sérsic model fitting
(decomposition fitting with AGN, bulge, and disk components)
or parametric/nonparametric substructure analysis.

With AGN-subtracted photometry of the host galaxy in
multiple bands, we estimate $M_*$ by performing SED fitting and
present the $M_{BH}-M_*$ relation at $z \sim 0.68-2.5$ (Figure 12). There
is a weak positive correlation between $M_{BH}$ and $M_*$ with the
correlation coefficient of $\rho = 0.25$ ($p = 0.010$). We fit the mass
relation by a simple (log-)linear relation of $\log(M_{BH}/M_{\odot}) =
\alpha \log(M_*/M_{\odot}) + \beta$ with consideration of selection biases and
measurement uncertainties. Our results show that the slope of the
mass relation at $z \sim 2$ ($\alpha = 0.89^{+0.61}_{-0.41}$) is consistent with the local
relation ($\alpha_{\text{local}} = 0.97^{+0.10}_{-0.11}$) (Figure 13).

Assuming the redshift evolution term of the mass relation to
be γ log (1 + z), we further determine the evolution factor γ
and the intrinsic scatter of the mass relation σ<sub>μ</sub> while
considering selection biases and uncertainties based compar-
isons to mock catalogs (Figure 14). Even though the estimated
probability distribution shows strong degeneracy between γ
and σ<sub>μ</sub>, we find no or mild evolution with γ = 0.48<sup>+0.31</sup><sub>-0.62</sub>. If we
assume that σ<sub>μ</sub> is not smaller than the local value, we obtain
γ = 0.22<sup>+0.39</sup><sub>-0.58</sub>, which is more consistent with the no- or mild-
evolution scenario. The estimated γ - σ<sub>μ</sub> distribution is largely
consistent with J. Li et al. (2021) based on the HSC imaging of
SDSS quasars at z < 0.8. Given the higher redshift range of our
sample, the slope of the degeneracy relation between γ and σ<sub>μ</sub>
is steeper than J. Li et al. (2021). Therefore, despite the sample
being approximately 6 times smaller, the estimated γ
uncertainty is just slightly larger than that from J. Li et al.
(2021).
Furthermore, the estimated value of the intrinsic scatter is $\sigma_{\mu} = 0.30_{-0.13}^{+0.14}$, which is consistent with the local relation and the recent estimate by J. Li et al. (2021). We show that this value $\sigma_{\mu}$ at high-$z$ may not be in contradiction to a cosmic averaging scenario (Figure 16) as recently put forward by J. Li et al. (2021) and X. Ding et al. (2022b) where AGN feedback is invoked to explain the constant level of dispersion with redshift. However, due to the small sample size, high redshift, and wide redshift range, our constraints on the redshift evolution of $\sigma_{\mu}$ are weak. Thus, a larger sample size at $z \sim 1-3$ is needed, especially at high-$z$.

Future large-scale surveys such as Euclid and Roman will significantly augment the sample size, along with deeper observations by JWST, to provide stronger constraints on SMBH and galaxy evolution. For future large imaging data sets, visual inspection for all multicomponent fits and manually exclusion of anomalous results will not be feasible; thus, improvements in 2D decomposition techniques or the imposition of more sophisticated conditions to confirm the robustness of host detection is needed. Additionally, by leveraging the high spatial resolution of JWST images, it is important to compare the morphology, substructures (shortly introduced in Section 4.3), and stellar populations of AGN-host galaxies with non-AGN galaxies, as well as to discuss the presence or absence of a bulge component and the $M_{BH}-M_{bulge}$ relation. In order to address these challenges, it is imperative to enhance the 2D decomposition analysis by mitigating the uncertainties related to PSF reconstruction. This involves conducting a meticulous analysis of AGN PSFs, determining the validity of applying stellar PSFs to AGN with different SEDs than stars, identifying the most effective methods for accurate PSF reconstruction, and establishing a framework for evaluating the uncertainties in reconstructed PSFs.

## Acknowledgments

This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. 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 NAS 5-03127 for JWST. These observations are associated with program IDs 1727 and 1837, and the specific observations analyzed can be accessed via DOI: 10.17909/g0bh-7854. Numerical computations were in part carried out on the iDark cluster, Kavli IPMU. This work was made possible by utilizing the CANDIDE cluster at the Institut d'Astrophysique de Paris. The cluster was funded through grants from the PNCG, CNES, DIM-ACAV, the Euclid Consortium, and the Danish National Research Foundation Cosmic Dawn Center (DNRF140). It is maintained by Stephane Rouberol. We thank Marko Shuntov, Kei Ito, and Mingyang Zhuang for giving us useful advice for reconstructing the PSFs. We thank Mingyang Zhuang, Yue Shen, and Junyao Li for giving us their measurements for comparison. We thank the anonymous referee for helpful feedback. Kavli IPMU is supported by World Premier International Research Center Initiative (WPI), MEXT, Japan. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant DNRF140. T.T. is supported by Forefront Physics and Mathematics Program to Drive Transformation (FoPM), a World-leading Innovative Graduate Study (WINGS) Program, the University of Tokyo. J.S. is supported by JSPS KAKENHI (JP22H01262) and the World Premier International Research Center Initiative (WPI), MEXT, Japan. This work was supported

by JSPS Core-to-Core Program (grant No. JPSCCA20210003). B.T. acknowledges support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 950533) and from the Israel Science Foundation (grant No. 1849/19). J.R. is supported by JPL, which is run by Caltech for NASA. G.E.M. and S.G. acknowledge the Villum Fonden research grant 13160 "Gas to stars, stars to dust: tracing star formation across cosmic time," grant 37440, "The Hidden Cosmos," and the Cosmic Dawn Center of Excellence funded by the Danish National Research Foundation under grant No. 140.

**Facilities:** JWST (NIRCam), HST (ACS).

**Software:** astropy (T. P. Robitaille et al. 2013; A. M. Price-Whelan et al. 2018; Astropy Collaboration et al. 2022), CIGALE (M. Boquien et al. 2019; G. Yang et al. 2022), FSPS (C. Conroy et al. 2009; C. Conroy & J. E. Gunn 2010), galight (X. Ding et al. 2020b), lenstronomy (S. Birrer & A. Amara 2018; S. Birrer et al. 2021), matplotlib (J. D. Hunter 2007), wCorr (P. Bailey & A. Emad 2023), numpy (C. R. Harris et al. 2020), photutils (L. Bradley et al. 2023), psfr (S. Birrer et al. 2025, in preparation), Prospector (J. Leja et al. 2017; B. D. Johnson et al. 2021), jwst (H. Bushouse et al. 2023).

## Appendix A
### Detected Numbers in Each Filter

Table 4 summarizes the number of undetected host galaxies based on the conditions described in Section 3.3 for each PSF and filter, while Table 5 then summarizes the number of detected hosts for each final PSF. Regarding the JWST filters, we have less galaxies that are classified as nondetection due to the BIC condition than due to the $S/N_{host}$ condition. The number of undetected hosts is lowest in F277W, followed by F150W and F444W. On the other hand, F115W and F814W show a larger number of undetected hosts, especially due to $S/N_{host}$. This trend is because F814W is the HST observation with lower resolution and shallower depth than the JWST observations, and both F814W and F115W are in the shorter-wavelength side of the Balmer break, leading to a tendency for a smaller intrinsic $H/T$.

<table><caption>**Table 4**
The Number of the Undetected Host Galaxies for Each Condition, PSF, and Filter</caption><thead><tr><th>Filter</th><th>PSF</th><th>(1)<br>BIC</th><th>(2)<br>S/N</th><th>(3)<br>Manual</th><th>Total</th></tr></thead><tbody><tr><td rowspan="3">F814W</td><td>Top-5</td><td>48</td><td>53</td><td>30</td><td>90</td></tr><tr><td>Top-75%</td><td>51</td><td>61</td><td>29</td><td>94</td></tr><tr><td>PSFEx</td><td>53</td><td>58</td><td>29</td><td>94</td></tr><tr><td rowspan="3">F115W</td><td>Top-5</td><td>1</td><td>20</td><td>1</td><td>21</td></tr><tr><td>Top-75%</td><td>1</td><td>21</td><td>1</td><td>22</td></tr><tr><td>PSFEx</td><td>2</td><td>21</td><td>1</td><td>22</td></tr><tr><td rowspan="3">F150W</td><td>Top-5</td><td>0</td><td>10</td><td>1</td><td>11</td></tr><tr><td>Top-75%</td><td>0</td><td>9</td><td>1</td><td>10</td></tr><tr><td>PSFEx</td><td>0</td><td>8</td><td>1</td><td>9</td></tr><tr><td rowspan="3">F277W</td><td>Top-5</td><td>0</td><td>1</td><td>1</td><td>2</td></tr><tr><td>Top-75%</td><td>0</td><td>0</td><td>1</td><td>1</td></tr><tr><td>PSFEx</td><td>0</td><td>1</td><td>1</td><td>2</td></tr><tr><td rowspan="3">F444W</td><td>Top-5</td><td>0</td><td>3</td><td>5</td><td>8</td></tr><tr><td>Top-75%</td><td>0</td><td>3</td><td>5</td><td>8</td></tr><tr><td>PSFEx</td><td>0</td><td>4</td><td>5</td><td>9</td></tr></tbody></table>
<table><caption>**Table 5**
The Number of Host-galaxy Detections for Each Final PSF</caption><thead><tr><th rowspan="3">PSF</th><th colspan="6">No. of Detected filters</th><th colspan="4">PRIMER</th></tr><tr><th colspan="6">COSMOS-Web</th><th colspan="4">PRIMER</th></tr><tr><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th></tr></thead><tbody><tr><td>Top-5</td><td>2</td><td>3</td><td>3</td><td>14</td><td>66</td><td>14</td><td>0</td><td>0</td><td>2</td><td>0</td><td>3</td></tr><tr><td>Top-75%</td><td>1</td><td>4</td><td>5</td><td>11</td><td>69</td><td>12</td><td>0</td><td>2</td><td>0</td><td>0</td><td>3</td></tr><tr><td>PSFEx</td><td>2</td><td>3</td><td>3</td><td>15</td><td>68</td><td>11</td><td>1</td><td>0</td><td>0</td><td>1</td><td>3</td></tr></tbody></table>

## Appendix B
Detailed Comparison of Different PSF Results

As mentioned in Section 3.1, 2D decomposition is significantly influenced by the differences in PSFs. In this study, we compare three different PSFs and discuss the PSF dependency of the results. Figure 17 (images on the left; (a)–(c)) shows each final PSF image for the four NIRCam filters. We fit each final PSF and each single PSF in the PSF library with the 2D Gaussian model, and measure the FWHM along a semimajor axis and ellipticity $b/a$ (defined as $FWHM_{minor}/FWHM_{major}$). Figures 17(d) and (e) show the distribution of FWHM and $b/a$ in each filter for an example AGN, CID-62, at $z_{spec} \sim 1.92$. As the number of stars used increases in the order of the top-5 stacked, the top-75% stacked, and PSFEx, we can see that the background noise is correspondingly lower. Regarding the FWHM distribution,

the top-5 and the top-75% PSFs are consistent with the FWHM distributions of single PSFs within the PSF library. PSFEx have consistent FWHMs in F115W and F150W, and slightly smaller FWHMs in F277W and F444W, suggesting the possibility of bias between automatic PSF selection by galight. In addition, each final PSF tends to have a higher $b/a$ than individual single PSFs. Regardless of whether $\chi^2$-based stacking or empirical modeling is employed, considering that the final PSF is a more reasonable PSF reconstruction than a single PSF, this suggests a potential bias toward a more elliptical PSF when using a single PSF. Next, Figure 18 compares the estimated $H/T$ using each PSF and filter (the full-filter version of Figure 8(c)). We can find strong positive correlations for all filters, with the distribution along the $y = x$ relation. However, when comparing the distribution with $y = x$, it is evident that PSFEx tends to estimate larger $H/T$ for F277W and F444W compared to top-5 and top-75% PSFs. As mentioned in Section 4.2, this trend can be interpreted by considering the differences in FWHM for each filter and PSF. As shown in Figure 4, PSFEx tends to exhibit smaller PSFs in F277W and F444W than the top-5 and top-75% PSFs. M.-Y. Zhuang & Y. Shen (2023) suggest that different PSFs result in larger $H/T$ in the 2D decomposition analysis. Therefore, it is conceivable that PSFEx shows larger $H/T$ in F277W and F444W with the smaller PSFs than the top-5 and top-75% PSFs.

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Comparison of different methods to evaluate empirical models of the PSF. Columns (a) to (c) show images of the top-5 stacked PSF, the top-75% stacked PSF, and the PSFEx-PSF for an example galaxy, CID-62 ($z = 1.92$), in F115W, F150W, F277W, and F444W, from top to bottom. FWHM along the semimajor axis measured from the 2D Gaussian fits is shown in the top left corner of each image. The white bars in the lower left indicate a scale of 1". Columns (d) and (e) show the distribution of the FWHM (semimajor axis) and ellipticity $b/a$. Gray histograms display the distribution of each single PSF in the library (Section 3.1.1). The values for the top-5 stacked PSF, the top-75% stacked PSF, PSFEx-PSF are marked by the blue, green, and red vertical lines, respectively.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Comparison of the estimated $H/T$ in each filter using different PSFs (columns (a)–(d)). The gray and red colors correspond to whether they were detected in each filter. Spearman's correlation coefficient $\rho$ for each prior is shown in the lower right corner of each panel. High correlation coefficients and the distribution around $y=x$ (black dashed line) suggest consistent results among the fitting with different PSFs.</figcaption></figure>

## Appendix C
Galight Test with Mock Image

We generate mock images of AGN and their host galaxies and perform 2D decomposition on these images to verify whether galight can accurately recover each parameter for the host galaxies. We use lenstronomy (S. Birrer & A. Amara 2018; S. Birrer et al. 2021) to generate 2000 NIRCam/F150W mock galaxy images. These are constructed as a superposition of a PSF and a smooth analytic Sérsic model based on randomly sampled parameters: host-to-total flux ratio ($H/T$), Sérsic index ($n$), effective radius ($r_e$), semimajor axis orientation ($\theta$), and ellipticity ($b/a$). The software tools lenstronomy and galight have different definitions of $r_e$, with the relation $r_{e,galight} = (b/a)^{-0.5} \times r_{e,lenstronomy}$. In this section and consistent with the main part of this paper, all $r_e$ are unified under the definition in galight. The specified range for each parameter is $\log(H/T)$ in [$\sim$ 5, 0], $n$ in [1, 6], $r_e$ in [$\sim$0.05, 1.6], $b/a$ in [0.2, 1.0], and $\theta$ in [0, $\pi$] (in radian). We assume the PSFEx PSF (Section 3.1.2) and that the center of host galaxy components is aligned with the center of PSF components. Additionally, we add Gaussian noise equivalent to the real NIRCam/F150W COSMOS-Web images. Using the generated mock images, we performed 2D decomposition with

the low-$\chi^2$ top-5 PSF reconstruction strategy (Section 3.1.1). Since the PSFEx PSF is assumed when generating the mock images, the results allow us to assess the performance in cases where the PSF may differ from the true PSF, similar to real data analysis. Note that we cannot conclude which PSF reconstruction method is superior from this analysis. In this test, we do not consider the effects of more complex morphological features than single Sérsic components, such as double Sérsic (disk+bulge) components, position offset between AGN and host, substructures (bars, spirals, and clumps), and nearby or merging objects.

Figures 19(a)–(c) compare the estimated major parameters ($H/T$, $r_e$, and $n$, respectively) with the input values. First, focusing on $H/T$, we can find the distinct change in the trend around the input $\log(H/T)$ value $\log(H/T_{\text{ans}}) \sim -2$. For $\log(H/T_{\text{ans}}) \lesssim -2$, the estimated value $\log(H/T_{\text{est}})$ remains relatively constant around -1.5, deviating significantly from the $y=x$ line. This implies that detecting weak galaxy components with $\log(H/T_{\text{ans}}) \lesssim -2$ is practically challenging using the 2D decomposition analysis with galight. On the other hand, for $\log(H/T_{\text{ans}}) \gtrsim -2$, $\log(H/T_{\text{est}})$ exhibits a consistent trend along the $y=x$ relation, while slightly larger than the input values. This trend is consistent with M.-Y. Zhuang &
<figure><img src="image_21.png" /><figcaption>**Figure 19.** Comparison of input (X-axis) and estimated (Y-axis) morphological parameters of mock galaxies for the single Sérsic mock data. Columns (a) to (c) compare log(H/T), $r_e$, and $n$, respectively, and the black dashed line indicate $y = x$. The gray and cyan circles indicate the all mock galaxies and the galaxies with log(H/T$_{est}$) > -1 corresponding to the minimum H/T$_{est}$ for our real galaxies (also see Figures 8, 18). Only mock galaxies with log(H/T$_{est}$) > -1 are used in plotting box plots. As demonstrated, we can accurately reconstruct host galaxy photometry and morphology.</figcaption></figure>

Y. Shen (2023), who reported that different PSFs tend to overestimate H/T. As shown in Figure 7(a), log(H/T$_{est}$) for the real detected galaxies in this study are log(H/T$_{est}$) ≳ -1. If we limit the mock results to log(H/T$_{est}$) ≳ -1, the difference between the estimation and input values (defined as estimation - input values) is Δlog(H/T) ~ 0.06 ± 0.09 dex, suggesting sufficiently consistent estimation of H/T. Next, examining $r_e$ and $n$, we find a consistent distribution along the $y = x$ relation, with an overall tendency to underestimate $r_e$ and overestimate $n$. We can also see that some objects hit the $r_e$ lower limit ($r_e$ ~ 0.05) and $n$ upper limit ($n$ ~ 7). If we limit the analysis to the above realistic cases (log(H/T$_{est}$) ≳ -1), as

shown in Figure 19(b), $r_e$ shows consistent results with Δr = -0.14 ± 0.17 and less hitting to the $r_e$ lower limit. However, for $n$, many objects still exhibit the $n$ upper limit even with log(H/T$_{est}$) ≳ -1. If we additionally cut the upper limit object with $n < 6.8$, the difference between the estimation and the input values is Δn ~ 1.4 ± 1.6, suggesting a consistent estimation with medium scatter. We confirm that larger-$n$ objects tend to have larger Δn, indicating the increasing difficulty of accurately reconstructing $n$ when they have larger $n$. Examining the residual H/T (defined as log(H/T$_{est}$) - log(H/T$_{ans}$)) as a function of $r_{e,ans}$ (Figure 19(d)), we confirm that, with log(H/T$_{est}$) > -1, the values of H/T can be recovered independently over a wide $r_e$
<figure><img src="image_22.png" /><figcaption>**Figure 20.** (a) Comparison of the input and estimated log (H/T) for the double Sérsic mock data. (b) The estimated *n* as a function of the input B/H. Blue and orange dashed lines indicate *n*<sub>disk</sub> and *n*<sub>bulge</sub>. (c) The log (H/T) residual as a function of the input B/H. In every panel, the gray and cyan circles indicate the all mock galaxies and the galaxies with log (H/T<sub>est</sub>) > −1 corresponding to the minimum H/T<sub>est</sub> for our real galaxies (also see Figures 8, 18). Only mock galaxies with log (H/T<sub>est</sub>) > −1 are used for the box plots. As shown, we can accurately reconstruct H/T flux ratio independent of bulge dominance B/H.</figcaption></figure>

range of ~0''1–1''2, effectively covering the typical $r_e$ of the sample (Figure 7(b)).

In this paper, we assume the PSF + single Sérsic model to fit the galaxy images. As discussed in Section 3.2, the single Sérsic model is a first-order approximation and cannot fully describe galaxy morphologies. For example, substructures such as bars and spiral arms are not considered, and both disk and bulge components cannot be described with a single Sérsic component. Even so, we expand the above mock tests using a PSF + double Sérsic model to address the effect of a bulge component on the PSF + single Sérsic model fitting. One of the double Sérsic components corresponds to a disk, and we assume $r_{e,\text{disk}} = 0''35$ and $n_{\text{disk}} = 1.0$. The other Sérsic component corresponds to a bulge, and we assume $n_{\text{bulge}} = 4.0$. $r_{e,\text{bulge}}$, $H/T$, bulge-to-host flux ratio $B/H$, $\theta_{\text{disk}}$, $\theta_{\text{bulge}}$, $b/a_{\text{disk}}$, and $b/a_{\text{bulge}}$ are randomly sampled. The specified range for each parameter is $\log(H/T)$ in [~ − 3, 0], $r_{e,\text{bulge}}$ in [~0''06, 0''2], $b/a_{\text{disk}}$ in [0.6, 1.0], $b/a_{\text{bulge}}$ in [0.8, 1.0], $\theta_{\text{disk}}$ in [0, $\pi$] (in radian), and $\theta_{\text{bulge}}$ in [0, $\pi$] (in radian). Using these models and randomly sampled parameters, we make mock data by following the same method with the above single Sérsic mock data. Then, we performed 2D decomposition with the low-$\chi^2$ top-5 PSF reconstruction strategy (Section 3.1.1), and the results are summarized in Figure 20. Figure 20 compares the estimated $H/T$ with the input values. Due to the inclusion of the bulge model, the deviation from the $y = x$ line around $\log(H/T) \sim -2$ is larger than for the single Sérsic mock data (Figure 19(a)). However, similar to the results for the single Sérsic mock data, $H/T$ is well reconstructed for $\log(H/T_{\text{est}}) > -1$. Figures 20(b)

and (c) indicate the B/H dependence of $n_{\text{est}}$ and $\log(H/T)$ residual. From Figure 20(b), as $B/H_{\text{ans}}$ increases and the system becomes more bulge dominated, the estimated *n* approaches values from 1 (=n<sub>disk</sub>) to 4 (=n<sub>bulge</sub>), indicating that the single Sérsic model fits the bulge for a bulge-dominated system. However, from Figure 20(c), *H/T* is well reconstructed independently on the *B/H* in the range of *B/H* ~ 0–0.6. Therefore, we conclude that we can reconstruct *H/T*, i.e., host photometry, well for systems where the bulge constitutes at most ~50% of the total host galaxy flux.

## Appendix D

### The SED of the Fitted PSF Components

In our decomposition analysis, we model the AGN as a PSF component to subtract its contribution from JWST imaging data. To verify the validity of our analysis, we examine how well the fitted PSF component aligns with the typical SED of AGNs. We use the photometries of the PSF components obtained from fitting with the top-5 PSF + Sérsic models to assess their SED. We fit the photometries of each PSF component with a single power-law model, $F_{\lambda} \propto \lambda^{\beta}$. Figure 21 shows the distribution of the fitted $\beta$, suggesting $\beta = -0.9_{-0.6}^{+0.7}$. While this analysis does not account for dust reddening effects that enlarge $\beta$ (e.g., T. S. Tanaka et al. 2024a), the result is broadly consistent with the typical spectral slope of AGNs ($\beta \sim -1.5$; D. E. Vanden Berk et al. 2001). Therefore, we conclude that our results align well with the typical SEDs of AGNs, supporting the robustness of our decomposition analysis.
<figure><img src="image_23.png" /><figcaption>**Figure 21.** The distribution of the spectral slope $\beta$ ($F_{\lambda} \propto \beta^{\beta}$) for the PSF components. The solid and dashed lines indicate the median and the 1$\sigma$ confidence levels.</figcaption></figure>

## ORCID iDs

Takumi S. Tanaka https://orcid.org/0009-0003-4742-7060

John D. Silverman https://orcid.org/0000-0002-0000-6977

Xuheng Ding https://orcid.org/0000-0001-8917-2148

Knud Jahnke https://orcid.org/0000-0003-3804-2137

Benny Trakhtenbrot https://orcid.org/0000-0002-3683-7297

Erini Lambrides https://orcid.org/0000-0003-3216-7190

Masafusa Onoue https://orcid.org/0000-0003-2984-6803

Irham Taufik Andika https://orcid.org/0000-0001-6102-9526

Angela Bongiorno https://orcid.org/0000-0002-0101-6624

Andreas L. Faisst https://orcid.org/0000-0002-9382-9832

Steven Gillman https://orcid.org/0000-0001-9885-4589

Christopher C. Hayward https://orcid.org/0000-0003-4073-3236

Michaela Hirschmann https://orcid.org/0000-0002-3301-3321

Anton Koekemoer https://orcid.org/0000-0002-6610-2048

Vasily Kokorev https://orcid.org/0000-0002-5588-9156

Zhaoxuan Liu https://orcid.org/0000-0002-9252-114X

Georgios E. Magdis https://orcid.org/0000-0002-4872-2294

Alvio Renzini https://orcid.org/0000-0002-7093-7355

Caitlin Casey https://orcid.org/0000-0002-0930-6466

Nicole E. Drakos https://orcid.org/0000-0003-4761-2197

Maximilien Franco https://orcid.org/0000-0002-3560-8599

Ghassem Gozaliasl https://orcid.org/0000-0002-0236-919X

Jeyhan Kartaltepe https://orcid.org/0000-0001-9187-3605

Daizhong Liu https://orcid.org/0000-0001-9773-7479

Henry Joy McCracken https://orcid.org/0000-0002-9489-7765

Jason Rhodes https://orcid.org/0000-0002-4485-8549

Brant Robertson https://orcid.org/0000-0002-4271-0364

Sune Toft https://orcid.org/0000-0003-3631-7176

## References

Akins, H. B., Casey, C. M., Lambrides, E., et al. 2024, arXiv:2406.10341

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Bailey, P., & Emad, A., 2023 wCorr: Weighted Correlations, Github, https://american-institutes-for-research.github.io/wCorr/

Baker, W. M., Tacchella, S., Johnson, B. D., et al. 2024, NatAs,

Beifiori, A., Courteau, S., Corsini, E. M., & Zhu, Y. 2012, MNRAS, 419, 2497

Bennert, V. N., Auger, M. W., Treu, T., Woo, J.-H., & Malkan, M. A. 2011a, ApJ, 742, 107

Bennert, V. N., Auger, M. W., Treu, T., Woo, J.-H., & Malkan, M. A. 2011b, ApJ, 726, 59

Bennert, V. N., Treu, T., Ding, X., et al. 2021, ApJ, 921, 36

Bertin, E. 2011, in ASP Conf. Ser. 442, Astronomical Data Analysis Software and Systems XX, ed. I. J. I. N. et al. (San Francisco, CA: ASP), 435

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Birrer, S., & Amara, A. 2018, PDU, 22, 189

Birrer, S., Shajib, A., Gilman, D., et al. 2021, JOSS, 6, 3283

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Bouchet, P., García-Marín, M., Lagage, P.-O., et al. 2015, PASP, 127, 612

Bradley, L., Sipcz, B., Robitaille, T., et al. 2023, astropy/photutils: 1.8.0 v1.8.0, Zenodo, doi:10.5281/zenodo.7946442

Brightman, M., Silverman, J. D., Mainieri, V., et al. 2013, MNRAS, 433, 2485

Brusa, M., Civano, F., Comastri, A., et al. 2010, ApJ, 716, 348

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2023, JWST Calibration Pipeline, v1.10.0, Zenodo, doi:10.5281/zenodo.7795697

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Cappelluti, N., Brusa, M., Hasinger, G., et al. 2009, A&A, 497, 635

Casey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2023, ApJ, 954, 31

Cen, R. 2015, ApJL, 805, L9

Chabrier, G. 2003, PASP, 115, 763

Cisternas, M., Jahnke, K., Bongiorno, A., et al. 2011, ApJL, 741, L11

Civano, F., Elvis, M., Brusa, M., et al. 2012, ApJS, 201, 30

Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833

Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486

DeGraf, C., Di Matteo, T., Treu, T., et al. 2015, MNRAS, 454, 913

Di Matteo, T., Colberg, J., Springel, V., Hernquist, L., & Sijacki, D. 2008, ApJ, 676, 33

Ding, X., Onoue, M., Silverman, J. D., et al. 2023, Natur, 621, 51

Ding, X., Silverman, J. D., & Onoue, M. 2022a, ApJL, 939, L28

Ding, X., Silverman, J., Treu, T., et al. 2020b, ApJ, 888, 37

Ding, X., Silverman, J. D., Treu, T., et al. 2022b, ApJ, 933, 132

Ding, X., Treu, T., Silverman, J. D., et al. 2020a, ApJ, 896, 159

Duras, F., Bongiorno, A., Ricci, F., et al. 2020, A&A, 636, A73

Elvis, M., Civano, F., Vignali, C., et al. 2009, ApJS, 184, 158

Fabian, A. C. 2012, ARA&A, 50, 455

Ferrarese, L., & Merritt, D. 2000, ApJL, 539, L9

Garilli, B., Le Fèvre, O., Guzzo, L., et al. 2008, A&A, 486, 683

Graham, A. W., Onken, C. A., Athanassoula, E., & Combes, F. 2011, MNRAS, 412, 2211

Greene, J. E., Labbe, I., Goulding, A. D., et al. 2024, ApJ, 964, 39

Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2009, ApJ, 698, 198

Habouzit, M., Li, Y., Somerville, R. S., et al. 2021, MNRAS, 503, 1940

Harikane, Y., Zhang, Y., Nakajima, K., et al. 2023, ApJ, 959, 39

Häring, N., & Rix, H.-W. 2004, ApJL, 604, L89

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Harrison, C. M. 2017, NatAs, 1, 0165

Hirschmann, M., Khochfar, S., Burkert, A., et al. 2010, MNRAS, 407, 1016

Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008, ApJS, 175, 356

Hunter, J. D. 2007, CSE, 9, 90

Ishino, T., Matsuoka, Y., Koyama, S., et al. 2020, PASJ, 72, 83

Ito, K., Tanaka, T. S., Shimasaku, K., et al. 2024b, arXiv:2408.08492

Ito, K., Valentino, F., Brammer, G., et al. 2024a, ApJ, 964, 192

Jahnke, K., & Macciò, A. V. 2011, ApJ, 734, 92

Jahnke, K., Bongiorno, A., Brusa, M., et al. 2009, ApJL, 706, L215

Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22

Jun, H. D., Im, M., Lee, H. M., et al. 2015, ApJ, 806, 109

Kashino, D., Silverman, J. D., Rodighiero, G., et al. 2013, ApJL, 777, L8

Kass, R. E., & Raftery, A. E. 1995, JASA, 90, 773

Kennedy, J., & Eberhart, R. 1995, in Proceedings of ICNN'95-Int. Conf. on Neural Networks, 4 (Piscataway, NJ: IEEE), 1942

Kocevski, D. D., Barro, G., McGrath, E. J., et al. 2023, ApJL, 946, L14

Kocevski, D. D., Finkelstein, S. L., Barro, G., et al. 2024, arXiv:2404.03576

Koekemoer, A. M., Aussel, H., Calzetti, D., et al. 2007, ApJS, 172, 196

Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511

Koutoulidis, L., Mountrichas, G., Georgantopoulos, I., Pouliasis, E., & Plionis, M. 2022, A&A, 658, A35

Lauer, T. R., Tremaine, S., Richstone, D., & Faber, S. M. 2007, ApJ, 670, 249

Le Fèvre, O., Cassata, P., Cucciati, O., et al. 2013, A&A, 559, A14

Le Fèvre, O., Vettolani, G., Garilli, B., et al. 2005, A&A, 439, 845

Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019, ApJ, 876, 3

Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170

Li, J., Silverman, J. D., Ding, X., et al. 2021, ApJ, 922, 142

Li, J., Silverman, J. D., Izumi, T., et al. 2022, ApJL, 931, L11

Li, J., Silverman, J. D., Merloni, A., et al. 2024b, MNRAS, 527, 4690
Li, J., Zhuang, M.-Y., & Shen, Y. 2024a, ApJ, 961, 19

Li, J. I. H., Shen, Y., Ho, L. C., et al. 2023, ApJ, 954, 173

Lilly, S. J., Le Brun, V., Maier, C., et al. 2009, ApJS, 184, 218

Lilly, S. J., Le Fèvre, O., Renzini, A., et al. 2007, ApJS, 172, 70

Lower, S., Narayanan, D., Leja, J., et al. 2020, ApJ, 904, 33

Lyke, B. W., Higley, A. N., McLane, J. N., et al. 2020, ApJS, 250, 8

Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285

Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2024, A&A, 691, A145

Marconi, A., & Hunt, L. K. 2003, ApJL, 589, L21

Matthee, J., Naidu, R. P., Brammer, G., et al. 2024, ApJ, 963, 129

Mechtley, M., Jahnke, K., Windhorst, R. A., et al. 2016, ApJ, 830, 156

Menci, N., Fiore, F., Bongiorno, A., & Lamastra, A. 2016, A&A, 594, A99

Merloni, A., Bongiorno, A., Bolzonella, M., et al. 2010, ApJ, 708, 137

Nardiello, D., Bedin, L. R., Burgasser, A., et al. 2022, MNRAS, 517, 484

Netzer, H., & Trakhtenbrot, B. 2007, ApJ, 654, 754

Oke, J. B. 1974, ApJS, 27, 21

Ono, Y., Harikane, Y., Ouchi, M., et al. 2023, ApJ, 951, 72

Onoue, M., Inayoshi, K., Ding, X., et al. 2023, ApJL, 942, L17

Pacifici, C., Iyer, K. G., Mobasher, B., et al. 2023, ApJ, 944, 141

Park, D., Woo, J.-H., Bennert, V. N., et al. 2015, ApJ, 799, 164

Peng, C. Y. 2007, ApJ, 671, 1098

Peng, C. Y., Impey, C. D., Ho, L. C., Barton, E. J., & Rix, H.-W. 2006, ApJ, 640, 114

Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, Proc. SPIE, 9143, 91433X

Perrin, M. D., & Soummer, R. 2012, Proc. SPIE, 8442, 84423d

Price-Whelan, A. M., Sipocz, B. M., Gnther, H. M., et al. 2018, AJ, 156, 123

Reines, A. E., & Volonteri, M. 2015, ApJ, 813, 82

Rieke, M. J., Kelly, D. M., Misselt, K., et al. 2023, PASP, 135, 028001

Robitaille, T. P., Tollerud, E. J., Greenfield, P., et al. 2013, A&A, 558, A33

Rodriguez-Gomez, V., Genel, S., Vogelsberger, M., et al. 2015, MNRAS, 449, 49

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525

Schneider, D. P., Richards, G. T., Hall, P. B., et al. 2010, AJ, 139, 2360

Schramm, M., & Silverman, J. D. 2013, ApJ, 767, 13

Schulze, A., Bongiorno, A., Gavignaud, I., et al. 2015, MNRAS, 447, 2085

Schulze, A., Schramm, M., Zuo, W., et al. 2017, ApJ, 848, 104

Schulze, A., Silverman, J. D., Kashino, D., et al. 2018, ApJS, 239, 22

Schulze, A., & Wisotzki, L. 2011, A&A, 535, A87

Schulze, A., & Wisotzki, L. 2014, MNRAS, 438, 3422

Schwarz, G. 1978, AnSta, 6, 461

Sérsic, J.L. 1963, BAAA, 6, 41

Shen, L., Lemaux, B. C., Lubin, L. M., et al. 2020, MNRAS, 494, 5374

Shen, Y. 2013, BASI, 41, 61

Shen, Y., Greene, J. E., Strauss, M. A., Richards, G. T., & Schneider, D. P. 2008, ApJ, 680, 169

Shen, Y., & Kelly, B. C. 2010, ApJ, 713, 41

Shen, Y., & Kelly, B. C. 2012, ApJ, 746, 169

Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, ApJS, 194, 45

Silverman, J. D., Green, P. J., Barkhouse, W. A., et al. 2005, ApJ, 618, 123

Silverman, J. D., Kashino, D., Sanders, D., et al. 2015, ApJS, 220, 12

Silverman, J. D., Mainieri, V., Ding, X., et al. 2023, ApJL, 951, L41

Simmons, B. D., Urry, C. M., Schawinski, K., Cardamone, C., & Glikman, E. 2012, ApJ, 761, 75

Simmons, B. D., Van Duyne, J., Urry, C. M., et al. 2011, ApJ, 734, 121

Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361, 776

Stone, M. A., Lyu, J., Rieke, G. H., & Alberts, S. 2023, ApJ, 953, 180

Stone, M. A., Lyu, J., Rieke, G. H., Alberts, S., & Hainline, K. N. 2024, ApJ, 964, 90

Suess, K. A., Bezanson, R., Nelson, E. J., et al. 2022, ApJL, 937, L33

Suh, H., Civano, F., Trakhtenbrot, B., et al. 2020, ApJ, 889, 32

Sun, M., Trump, J. R., Brandt, W. N., et al. 2015, ApJ, 802, 14

Sersic, J. L. 1968, Atlas de Galaxias Australes (Cordoba: Observatorio Astronomico),

Tanaka, T. S., Silverman, J. D., Nakazato, Y., et al. 2024a, PASJ, 76, 1323

Tanaka, T. S., Silverman, J. D., Shimasaku, K., et al. 2024b, arXiv:2412.14246

Trakhtenbrot, B., & Netzer, H. 2012, MNRAS, 427, 3081

Vanden Berk, D. E., Richards, G. T., Bauer, A., et al. 2001, AJ, 122, 549

Vestergaard, M., & Peterson, B. M. 2006, ApJ, 641, 689

Wang, B., de Graaff, A., Davies, R. L., et al. 2024, arXiv:2403.02304

Weaver, J. R., Davidzon, I., Toft, S., et al. 2023, A&A, 677, A184

Yang, G., Boquien, M., Brandt, W. N., et al. 2022, ApJ, 927, 192

Yue, M., Eilers, A.-C., Simcoe, R. A., et al. 2024, ApJ, 966, 176

Zhang, Y., Ouchi, M., Gebhardt, K., et al. 2023, ApJ, 948, 103

Zhuang, M.-Y., & Ho, L. C. 2023, NatAs, 7, 1376

Zhuang, M.-Y., Li, J., & Shen, Y. 2024a, ApJ, 962, 93

Zhuang, M.-Y., & Shen, Y. 2024b, ApJ, 962, 139

Zou, F., Yang, G., Brandt, W. N., & Xue, Y. 2019, ApJ, 878, 11