Our knowledge of relations between supermassive black holes and their host galaxies at zgtrsim1 is still limited even though being actively sought out to zsim6. Here we use the high resolution and sensitivity of JWST to measure the host galaxy properties for 107 Xrayselected typeI AGNs at 0.68ltzlt2.5 with restframe opticalnearinfrared imaging from COSMOSWeb and PRIMER. Black hole masses logleftMrm BHModotrightsim6.99.6 are available from previous spectroscopic campaigns. We extract the host galaxy components from four NIRCam broadband images and the HSTACS F814W image by applying a 2D image decomposition technique. We detect the host galaxy for sim90 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 logleftMModotrightsim9.511.6 through SED 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 Mmathrm BHM ratio evolves as left1zright0.480.620.31 thus remains essentially constant or exhibits mild evolution up to zsim2.5. We also see an amount of scatter sigmamu0.300.140.13 similar to the local relation and consistent with lowz studies and a noncausal cosmic assembly history where mergers contribute to the statistical averaging towards the local relation is still feasible. We highlight improvements to come with larger samples from JWST and particularly Euclid which will exceed the statistical power of current wide and deep surveys.
['Astrophysics - Astrophysics of Galaxies']
The Mrm BHM relation up to zsim2 through decomposition of COSMOSWeb NIRCam images
{'No Header': '9', 'The M BH -M ∗ relation up to z ∼ 2 through decomposition of COSMOS-Web NIRCam images': "<!-- image -->  \n- 1 Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan\n- 2 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\n- 3 Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan  \n4  \nCenter for Astrophysical Sciences, Department of Physics & Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA  \n5  \nMax-Planck-Institut fur Astronomie, Konigstuhl 17, D-69117 Heidelberg, Germany  \n6 School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel  \n7 NASA-Goddard Space Flight Center, Code 662, Greenbelt, MD, 20771, USA  \n8 Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, China  \nTechnical University of Munich, TUM School of Natural Sciences, Department of Physics, James-Franck-Str. 1, D-85748 Garching,  \nGermany  \n10 Max-Planck-Institut fur Astrophysik, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany  \n11  \nINAF-Osservatorio Astronomico di Roma, Via Frascati 33, I-00078, Monte Porzio Catone, Italy  \n12 Caltech/IPAC, 1200 E. California Blvd. Pasadena, CA 91125, USA  \n13 Cosmic Dawn Center (DAWN), Denmark  \n14 DTU Space, Technical University of Denmark, Elektrovej, Building 328, 2800, Kgs. Lyngby, Denmark  \n- 15 Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA\n- 16 Institute for Physics, Laboratory for Galaxy Evolution and Spectral Modelling, EPFL, Observatoire de Sauverny, Chemin Pegasi 51, 1290 Versoix, Switzerland  \n17 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  \n18 Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, NL-9700AV Groningen, the Netherlands  \n19 DTU-Space, Technical University of Denmark, Elektrovej 327, DK2800 Kgs. Lyngby, Denmark  \n20  \nNiels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N, Denmark  \n21 Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, Padova, I-35122, Italy  \n22 The University of Texas at Austin, 2515 Speedway Boulevard Stop C1400, Austin, TX 78712, USA  \nDepartment of Physics and Astronomy, University of Hawaii, Hilo, 200 W Kawili St, Hilo, HI 96720, USA  \n24  \nDepartment of Computer Science, Aalto University, PO Box 15400, Espoo, FI-00 076, Finland  \n25  \nDepartment of Physics, Faculty of Science, University of Helsinki, 00014-Helsinki, Finland  \n26 Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA  \n27 Max-Planck-Institut fur extraterrestrische Physik, Gießenbachstraße 1, 85748 Garching b. Munchen, Germany 28 Institut d'Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Universit /acute.ts1 e, 98 bis boulevard Arago, F-75014 Paris, France 29 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91001, USA 30  \n- Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA", 'ABSTRACT': 'Our knowledge of relations between supermassive black holes and their host galaxies at z ≳ 1 is still limited, even though being actively sought out to z ∼ 6. Here, we use the high resolution and sensitivity of JWST to measure the host galaxy properties for 107 X-ray-selected type-I AGNs at  \ntakumi.tanaka@ipmu.jp  \n23  \n0 . 68 < z < 2 . 5 with rest-frame optical/near-infrared imaging from COSMOS-Web and PRIMER. Black hole masses (log ( M BH /M ⊙ ) ∼ 6 . 9 -9 . 6) are available from previous spectroscopic campaigns. We extract the host galaxy components from four NIRCam broadband images and the HST/ACS F814W image by applying a 2D image decomposition technique. We detect the host galaxy for ∼ 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 ⊙ ) ∼ 9 . 5 -11 . 6 through SED 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 . 31 -0 . 62 thus remains essentially constant or exhibits mild evolution up to z ∼ 2 . 5. We also see an amount of scatter ( σ µ = 0 . 30 +0 . 14 -0 . 13 ), similar to the local relation and consistent with lowz studies, and a non-causal cosmic assembly history where mergers contribute to the statistical averaging towards the local relation is still feasible. We highlight improvements to come with larger samples from JWST and, particularly, Euclid, which will exceed the statistical power of current wide and deep surveys.  \nKeywords: AGN host galaxies (2017) - Active galactic nuclei (16) - Active galaxies (17) - Galaxy evolution (594)', '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 σ ∗ and stellar mass M ∗ (Magorrian et al. 1998; Ferrarese & Merritt 2000; Marconi & Hunt 2003; Haring & Rix 2004; Gultekin et al. 2009; Graham et al. 2011; Beifiori et al. 2012; Kormendy & Ho 2013; Reines & Volonteri 2015). How and through what physical processes such a tight relation is formed is still unclear, and the origin of the mass relation can shed light on the evolution of not only SMBHs but also galaxies.  \nA widely considered scenario, in answer to this question, is a co-evolution scheme, where galaxies and black holes mutually increase their mass at a correlated pace. As a potential physical cause for a co-evolution scenario, some studies implement active galactic nuclei (AGN) feedback, where the energy released from AGNs heats the gas and controls star formation or gas accretion through radio jets or an AGN winds (e.g. Springel et al. 2005; Di Matteo et al. 2008; Hopkins et al. 2008; Fabian 2012; DeGraf et al. 2015; Harrison 2017). Additionally, studies support a common gas supply simultane-  \nously fueling both SMBHs and their host galaxies by increasing the BH accretion rate and star formation rate (SFR, Cen 2015; 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 non-casual connection; major mergers have averaged the mass relation statistically (cosmic averaging scenario; Peng 2007; Hirschmann et al. 2010; Jahnke & Macci'o 2011).  \nTo 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., Ding et al. 2020; Habouzit et al. 2021; 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.  \nBefore 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., Peng et al. 2006; Jahnke et al. 2009; Bennert et al. 2011a; Cisternas et al. 2011; Simmons et al. 2011, 2012; Schramm & Silverman 2013; Mechtley et al. 2016; Ding et al. 2020; Bennert et al. 2021; Li et al. 2023b) and ground-based telescopes including Subaru's Hyper Suprime-Cam (HSC) (Ishino et al. 2020; Li et al. 2021), and Pan-STARRS1 (PS1) (Zhuang & Ho 2023). In summary, these studies have concluded that the relation between M BH and M ∗ does not evolve with redshift at z ≲ 2. However, studies using a large statistical and universal sample have yet to be  \nachieved at z > 1. At these redshifts, we can obtain information longer than the 4000 ˚ A 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 ≳ 1 (e.g., Jahnke et al. 2009; Simmons et al. 2012; Mechtley et al. 2016; Ding et al. 2020, utilizing HST/NICMOS or HST/WFC3). Other studies carry out statistical AGN samples using spectral energy distribution (SED) fitting based 1D decomposition method (e.g., Merloni et al. 2010; Sun et al. 2015; Suh et al. 2020).  \nJWST is now revolutionizing the field of AGN - host galaxy relations up to z ∼ 6 and beyond based on its high spatial resolution and unprecedented sensitivity. For instance, 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 ∼ 1 . 6 -3 . 5 from the SDSS DR17Q catalog (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 ∼ 1 . 6. Other studies have also performed 2D decompositions of AGN host galaxies using JWST imaging; for instance, Li et al. (2023a) have analyzed a galaxy that is one of the most promising candidates for having a recoiling SMBH ( z ∼ 0 . 36) while Kocevski et al. (2023) present the host properties of five X-ray-luminous AGNs (3 < z < 5) in CEERS. 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).  \nAt z > 6, 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 ∼ 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 highz AGNs are revealing a higher abundance of lower mass black holes that are actively accreting within very dusty and compact galaxies (e.g., Onoue et al. 2023; Kocevski et al. 2023; Matthee et al. 2023; Harikane et al. 2023; Maiolino et al. 2023; Greene et al. 2023; Kocevski et al. 2024; Wang et al. 2024; Akins et al. 2024). However, in general, the sample sizes of these highz AGNs and the accuracy of M ∗ estimation are still limited. Therefore, the redshift evolution of the mass relation remains highly uncertain.  \nTherefore, 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  \nwhat was statistically analyzed before the JWST era, further bridging the gap between lowz statistical and limited highz studies. We use a sample larger than previous studies using JWST that reaches up to z ∼ 2 -3 when AGN and star formation activities peaked in cosmic history. With NIRCam images of N = 107 AGNs in COSMOS-Web (Casey et al. 2023) and PRIMERCOSMOS, 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 (Lauer et al. 2007; Shen & Kelly 2010; Schulze & Wisotzki 2011, 2014). In addition, we report on the ability to accurately model the JWST PSF in each band and the impact on the derived host galaxy properties.  \nThis 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 constrcution 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 magnitude are AB magnitude (Oke 1974), and we assume a standard cosmology with H 0 = 70 km s -1 Mpc -1 , Ω m = 0 . 30, and Ω Λ = 0 . 70.", '2.1. COSMOS-Web': "COSMOS-Web (PI: Jeyhan Kartaltepe and Caitlin Casey, GO1727, see Casey et al. 2023, for the overview) is a 270-hour treasury survey program in JWST Cycle 1, covering 0.54 deg 2 with NIRCam (Rieke et al. 2023) in four filters (F115W, F150W, F277W, F444W) and 0.19 deg 2 with MIRI (Bouchet et al. 2015) using F770W. Due to the large field, the COSMOS-Web field was split into twenty tiles.  \nThe data are reduced with the JWST Calibration Pipeline 1 (Bushouse et al. 2023) version 1.14.0 and the calibration Reference Data System version 1223. The 5 σ depth in an aperture with a radius of 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 Casey et al. 2023). The mosaic images  \nFigure 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 (Koekemoer et al. 2007) for the COSMOS field. Red stars show the position of each AGN. The region enclosed by the solid red line shows the PRIMER-COSMOS field, and gray shaded region corresponds to the COSMOS-Web field.  \n<!-- image -->  \nhave a resolution of 0 . '' 030/pixel. Details of the image reduction process will be described in Franco et al.,(in preparation). In addition to NIRCam four-band data, we use HST/ACS F814W data (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': "Public Release IMaging for Extragalactic Research (PRIMER, PI: James Dunlop, GO1837) is a 195-hour treasury program of JWST Cycle 1, targeting two equatorial HST CANDELS Legacy Fields: COSMOS and UDS. PRIMER-COSMOS covers 144 arcmin 2 with eight NIRCam (Rieke et al. 2023) filters (F090W, F115W, F150W, F200W, F277W, F356W, F410M, F444W) and 112 arcmin 2 with two MIRI (Bouchet et al. 2015) filters (F770W and F1800W) in the COSMOS field. The processed COSMOS-PRIMER data consists 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 (Bushouse et al. 2023) version 1.8.3 and the calibration Reference Data System version 1017. The 5 σ depth in an aperture with a radius of 0 . '' 15 has a wide range  \nfrom 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 . '' 030/pixel.", '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 Schulze et al. (2015, 2018). Schulze et al. (2015) presents the redshift evolution of AGN population based on spectroscopically observed type-I AGNs from zCOSMOS (Lilly et al. 2007, 2009), VVDS (Le F'evre et al. 2005, 2013; Garilli et al. 2008), and SDSS (Schneider et al. 2010; Shen & Kelly 2012). Schulze et al. (2018) provides the properties of X-ray selected and spectroscopically-confirmed type-I AGNs in the FMOSCOSMOS survey (Kashino et al. 2013; Silverman et al. 2015). Here, we select the targets in Schulze et al. (2015) and Schulze et al. (2018) that are also detected by Chandra (Chandra-COSMOS Survey; Elvis et al. 2009; Civano et al. 2012) or XMM-Newton (XMM-COSMOS; Cappelluti et al. 2009; Brusa et al. 2010). The 2-10 keV flux sensitivity is 7 . 3 × 10 -16 erg cm -2 s -1 for Chandra and 3 × 10 -15 erg cm -2 s -1 for XMM-Newton.  \nWe have M BH estimates from spectra 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 β line is used for calibrating virial mass estimators, we use the M BH from H β (FMOS-COSMOS), H α (FMOS-COSMOS), and Mg ii (zCOSMOS-Bright, Deep) in order of preference; e.g., FMOS H β is used for an object with both FMOS H β and zCOSMOS Mg ii estimation. As shown in Figure 10 of Schulze et al. (2018), there is a very good agreement between the FMOS H α - and H β -based M BH compared to those using MgII. Schulze et al. (2018) also compared FWHM measurements taken at different times for each object (see Fig.7 of 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 (Lyke et al. 2020); we do not use these SDSS Mg ii -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.  \nFor H β , Schulze et al. (2018) used the virial mass estimation relation by Vestergaard & Peterson (2006),  \nM BH (H β ) = 10 6 . 91 ( L 5100 10 44 erg s -1 ) 0 . 5 ( FWHM H β 1000 km s -1 ) 2 M ⊙ , (1)  \nFigure 2 shows the distribution of 2-10 keV X-ray luminosity L [2 -10 keV] (panel a) and M BH (panel b) as a function of redshift, where L [2 -10 keV] is calculated with an X-ray spectral index Γ = 1 . 8 (e.g., Brightman et al. 2013). Since the sample consists of X-ray-selected objects and is flux-limited, there is a tendency for higherz objects to have larger L [2 -10 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 propor-  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n/circledot  \nFigure 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 -10 keV] as a function of z . Black dashed lines indicate the F [2 -10 keV] 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.  \nTable 1. Sample size for each single-epoch M BH estimation  \nwhere L 5100 is continuum luminosity at 5100 ˚ A and FWHM H β is the full-width at half-maximum (FWHM) of H β broad line. Then, the H α -based masses are calculated as given in Schulze et al. (2017); Schulze et al. (2018),  \nM BH (H α ) = 10 6 . 71 ( L H α 10 42 erg s -1 ) 0 . 48 ( FWHM H α 1000 km s -1 ) 2 . 12 M ⊙ , (2)  \nL H α is the H α luminosity and FWHM H α is the FWHM of the broad H α emission line. For Mg ii -based M BH estimation (zCOSMOS-bright Mg ii , and zCOSMOSdeep Mg ii ), the calibration by Shen et al. (2011) is used,  \nM BH (MgII) = 10 6 . 74 ( L 3000 10 44 erg s -1 ) 0 . 62 ( FWHM MgII 1000 km s -1 ) 2 M ⊙ , (3)  \nwhere L 3000 is continuum luminosity at rest-frame 3000 ˚ A and FWHM MgII is FWHM of Mg ii 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 (c.f., 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).  \nFrom 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 broad-band images from HST/ACS (F814W) and JWST/NIRCam (F115W, F150W, F277W, and F444W) for the sample residing in COSMOS-Web. Four more broad-band (F090W, F200W, F356W) and medium-band (F410M) images are available for five galaxies in the PRIMER field.  \nWe also compare the optical color index g -i calculated based on the COSMOS2022 photometry (Weaver et al. 2022) with SDSS quasars and hard-X-ray-detected AGNs at the same redshift range (Figure 6 in 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.  \ntional 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 .  \nFigure 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.  \nFigure 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.  \n<!-- image -->", '3. METHOD': "To extract host galaxy components from the original AGN + host galaxy NIRCam images, we apply a 2D image analysis tool galight (Ding et al. 2020). With galight , we perform forward modeling of each image as a superposition of a PSF component and PSF-convolved S'ersic components corresponding to the light from an AGN and its host galaxy, respectively. We then obtain images of the host galaxy, free 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., Ding et al. 2020, 2022a),  \nthe 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., Suess et al. 2022) or stellar images (e.g., Nardiello et al. 2022; Ding et al. 2022a; Ding et al. 2023; Zhuang & Shen 2023; Baker et al. 2023). The former uses the theoretical PSF model such as WebbPSF (Perrin et al. 2012, 2014), and the latter uses natural stellar images as the PSF directly or modeled PSF with tools such as PSFEx (Bertin 2011). Many previous studies concluded that the synthetic PSF simulated by WebbPSF is narrower than the PSFs reconstructed with the natural stars (Ono et al. 2022; Ding et al. 2022a; Onoue et al. 2023). Note that 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. χ 2 ν -based methods': "First, we follow the strategy of Ding et al. (2020, 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 χ 2 ν and stack the PSFs with the top 3, 5, and 8 χ 2 ν values. Using the single PSFs and the stacked PSFs, they select the PSF with the smallest χ 2 ν as the final PSF. This method is based on the χ 2 ν ; i.e., they assumed that the lower χ 2 ν is indicative of a better (the closer to the more accurate) PSF.  \nFollowing 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-wavelengthchannel 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,  \nwe fit each AGN target with a superposition of a PSFconvolved S'ersic 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 χ 2 ν and top-75% χ 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 χ 2 ν selected for each target.", '3.1.2. Modeling method': "Zhuang & 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, Zhuang & Shen (2023) also suggested that smaller χ 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 Zhuang & Shen (2023), we use PSFEx and compare the results with χ 2 ν -based selected PSFs (Section 3.1.1).  \nPSFEx constructs an empirical PSF model based on the output catalog of SExtractor (Bertin & 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 Zhuang & 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 Zhuang & Shen (2023) tested local PSF reconstruction.", '3.1.3. Comparing the final PSFs': 'Now, we have three final PSFs for comparison; the top-5 χ 2 ν stacked PSF, the top-75% χ 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 semi-major axis. Figure 4 compares the FWHMs of PSFs obtained with each method, target, and filter. The x and y-axis of Figure 4 show the ratio FWHM Top -5 /FWHM PSFEx and FWHM Top -75% /FWHM PSFEx , where FWHM is the value along the semi major axis.  \nFirstly, regardless of the filters, we can see that the distribution extends further 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  \nFigure 4. Comparison of FWHM (semi-major 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 yaxis 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 σ 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 top75% PSFs. Notably, F277W and F444W exhibit a FWHM bias among different PSF reconstruction methods.  \n<!-- image -->  \nthe 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 χ 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 Top -5 /FWHM PSFEx have a larger scatter than FWHM Top -75% /FWHM PSFEx .  \nSecondly, we focus on the FWHM bias between the methods for each filter. As suggested by Zhuang & 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.  \nNote 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 distribu-  \ntion 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'ersic profile (S'ersic 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'ersic components and fitted simultaneously. Our S'ersic model has seven free parameters: amplitude, S'ersic 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'ersic component is ten 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'ersic profile.  \nTo 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'ersic profile, such as bars and spiral arms. Thus, using a S'ersic profile is a first-order approximation to model the global component of AGN host galaxies.  \nIn the fitting process, we cut the image into square regions centered on the target with a radius seven times the standard deviation along the semimajor axis of the 2D Gaussian fitting with photutils (Bradley et al. 2023). Then, the above model is optimized with Particle Swarm Optimizer (PSO; Kennedy & Eberhart 1995). galight also supports Markov Chain Monte Carlo (MCMC) to estimate the posterior parameter distributions. As suggested by 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  \nwith 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 (c.f., Ding et al. 2023).  \nAs described above, using the PSF library constructed with galight , we fit with every single PSF in the library and sort them in order of χ 2 ν . Then we fit with the three final PSFs; top-5 χ 2 ν stacked, top 75% χ 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.  \nFor 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 ˚ A 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 ∼ 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.  \n(1) Bayesian Information Criteria: In addition to the PSF + S'ersic 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; Schwarz 1978) for the two models, PS+SE and PS model, as,  \nBIC = χ 2 + k ln ( n ) , (4)  \nwhere 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 PS+SE is much smaller than BIC PS ,  \nas  \nBIC PS+SE < BIC PS -10 , (5)  \nFigure 5. Example fits for CID-62 ( z ∼ 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 normalized residuals, from left to right. χ 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 disk-like host galaxy which is buried under the PSF before subtraction, regardless of the PSF reconstruction method.  \n<!-- image -->  \nWe decide the threshold value of the BIC difference of 10 based on Kass & Raftery (1995).  \n(2) S/N of the host galaxy: To estimate the significance of the detection of the host galaxy, we calculate the S/N of the host galaxy as done in 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 FoV (Zhuang & Shen 2023; Yue et al. 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 signal-  \nto-noise ratio of host galaxy S/N host within a radius of 2 r e . Then, we define the detection as cases with high S/N host , as  \nS/N host > 5 . (6)  \n(3) Manual inspection and removal: We find that some objects have invalid central values 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 non-detections.  \nCID-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 mismathced 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.  \nWe 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  \nfeatures 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 non-detection cases.  \nWith the above three strategies, we confidently report the detection of an AGN host galaxy for those that fulfill both the conditions (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 top5 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'ersic fits, considering Galactic dust extinction (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., Ding et al. 2022a; Zhuang & Shen 2023; Zhuang et al. 2023) and errors discussed in Section 3.3, considering both observational errors and PSF uncertainty in a radius of 2 r e . We use the 3 σ value for the undetected filters as the upper limit.  \nHere, we also fit with a model containing only a S'ersic component (SE model). For targets with BIC PS+SE > BIC SE -10, the SE model is better or flexible enough to describe the data than the PS+SE model, and we use the S'ersic 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 , Boquien et al. 2019; Yang et al. 2022) SED fitting library. For a stellar population, we use the single stellar population model by Bruzual & Charlot (2003) ( bc03 module) and the Chabrier initial mass function (IMF; Chabrier 2003) with the M ∗ cutoff of 0 . 1 M ⊙ and 100 M ⊙ (Bruzual & Charlot 2003). We assume a delayedτ model for a star-formation history (SFH), where SFR at each look-back time t is modeled as  \nSFR ( t ) ∝    ( t -t age ) exp ( -t -t age τ ) ( t > t age ) , 0 ( t < t age ) , (7)  \nwhere t age 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 Calzetti et al. (2000) law. We set E ( B -V ), M ∗ , t age , and τ as free parameters; their grid values are decided basically following Zhuang et al. (2023) CIGALE run and summarized in Table 2. To avoid unphysical solutions, we set the upper limit of t age to 0 . 95 t H , where t H indicates the cosmic age at each redshift. Otherwise, stellar metallicity and the ionization parameter U are fixed at Solar metallicity and log U = -2.  \nWe also apply the Bayesian-based spectral energy density (SED) fitting code, Prospector (Leja et al. 2017; Johnson et al. 2021), to assess the uncertainty M ∗ with different SED fitting codes. Prospector is based on the Flexible Stellar Population Synthesis ( FSPS , Conroy et al. (2009); Conroy & 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, Calzetti et al. (2000) dust attenuation law, and nebular emission with log U = -2. Prospector can fit with the non-parametric SFH, which separates galaxy formation history into several age bins and assumes a constant SFR in each bin (e.g., Leja et al. 2019). Lower et al. (2020) input cosmological hydrodynamic simulation data into Prospector and concluded that non-parametric SFH tends to reconstruct M ∗ more accurately than parametric SFHs. 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 non-parametric 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 BH or higher Eddington ratio (Lauer et al. 2007; Schulze & 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 Li et al. (2021) and apply the mock  \nTable 2. Major parameter settings in SED fitting  \nFigure 6. Procedure for generating mock galaxies at z ∼ 1 . 5 with an example parameter set of α = α local , β = β local , γ = 0, and σ µ = 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 BH , mock vs. true stellar mass M BH , mock , before adding observational uncertainties, (b) 2 -10 keV flux ( F [2 -10 keV]) , mock ), (c) H α FWHM (FWHM H α, mock ), (d) pseudo-observed virial mass M vir , mock vs. true black hole mass M BH , mock where M vir , mock is the mock black hole mass after considering the observational bias. (e) pseudo-observed stellar mass ˜ M ∗ , mock vs. stellar mass M ∗ , mock where ˜ M ∗ , mock is the stellar mass after considering the observational bias, and (f) pseudoobserved virial mass M vir , mock vs. pseudo-observed stellar mass ˜ M ∗ , mock . The latter is to be compared with the observed M BH -M ∗ plane shown in the left panel of Figure 12. Black dashed lines in panels (a) and (f) indicate M BH -M ∗ local relation obtained from the fitting of the local galaxies (Haring & Rix 2004; Bennert et al. 2011b, , see Section 5), and gray dashed lines in panels (d) and (e) indicate y = x line.  \n<!-- image -->  \n∗  \n∗  \nobservation (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.  \nFirst, we generate the mock redshift z mock and mock true stellar mass M ∗ , mock based on the COSMOS2020 stellar mass function (SMF) by Weaver et al. (2023).  \nNext, we use the M ∗ , mock to generate the mock true BH mass M BH , mock for the mock sample. Here, we assumed the local relation as  \nlog ( M BH M ⊙ ) = α log ( M ∗ M ⊙ ) + β, (8)  \nwhere α and β indicate the slope and the intercept of the local M BH -M ∗ plane. Then, assuming a normal  \ndistribution, we calculate M BH , mock as,  \nlog M BH , mock = N ( α log ( M ∗ M ⊙ ) + β + γ log (1 + z ) , σ 2 µ ) . (9)  \nHere, the parameter γ indicates the strength of redshift evolution of the mass relation, and σ µ is the intrinsic scatter of the mass relation. Figure 6 (a) shows M BH , mock vs. M ∗ , mock distribution with the example parameter set of γ = 0 and σ µ = 0 . 3.  \nThen, based on z mock and M BH , mock , we assign mock Eddington ratio λ Edd , mock by sampling the Eddington ratio distribution function by Schulze et al. (2015). Using the λ Edd , mock and M BH , mock , we calculate the bolometric luminosity L bol , mock .  \nFrom L bol , mock , we obtain L [2 -10 keV] , mock using the bolometric correction by Duras et al. (2020). With the calculated L [2 -10 keV] , mock , mock X-ray flux F [2 -10 keV] , mock is determined with the assumption of Γ = 1 . 8, the same assumption as in Section 2.3. Figure 6 (b) shows F [2 -10 keV] , mock distribution with the example parameter set of γ = 0 and σ µ = 0 . 3.  \nWe also generate mock virial BH mass M vir , mock with the L bol , mock , M BH , mock , and the assumed FWHM distribution. We first calculate mock continuum and line luminosity L 5100 , mock , L 3000 , mock , and L H α, mock . For L 5100 , mock and L 3000 , mock , we used the bolometric correction by Netzer & Trakhtenbrot (2007) and Trakhtenbrot & Netzer (2012), respectively. For L H α, mock , we use the L 5100 , mock and scaling relation between L H α and L 5100 , mock by Jun et al. (2015). These bolometric correction and scaling relations are the same as those in 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 (Shen et al. 2008), and we generate the observed FWHM of H α , H β , 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 β bias = 0 . 6 in calculating FWHM. β bias represents the proportion by which ∆ L , the variation of luminosity from the mean luminosity ¯ L , affects the variation in FWHM. Thus, the FWHM is generated by a lognormal distribution with the scatter of 0.17 dex and the mean value corresponding to the luminosity of ¯ L + β ∆ L (c.f. Shen 2013; Li et al. 2021). Finally, the mock virial BH mass M 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 α, mock histogram and M vir , mock -M BH , mock distribution with the example parameter set of γ = 0 and σ µ = 0 . 3.  \nConsidering possible systematic uncertainties from 2D decomposition analysis, we added an error based on a  \nnormal distribution with 0.2 dex to the M ∗ , mock to consider the uncertainty of M ∗ derived from observations, resulting in the mock observed stellar mass ˜ M ∗ , mock . The resulted ˜ M ∗ , mock -M ∗ , mock distribution is shown in Figure 6 (e).  \nFinally, we get the mock data set of the redshift z mock , observed stellar mass ˜ M ∗ , mock , and virial BH mass M vir , mock . Figure 6 (f) shows the final mock observed mass distribution ( M vir , mock -˜ M ∗ , mock ) with the example parameter set of γ = 0 and σ µ = 0 . 3. In Section 5, we discuss the mass relation, comparing our results with this mock catalog.', '4.1. Morphological parameters': 'Firstly, we examine the results obtained from the 2D decomposition. In Figure 7 (a)-(c), we show the distribution of n and r e (in the unit of pixel and kpc) for each filter using the top-5 stacked PSF. The distribution of n has a peak in the distribution at n ∼ 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 disk-like morphology.  \nIn Figure 7 (d), we compare H/T with S/N host . We can see a strong correlation between the reconstructed H/T and S/N 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 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 Figure 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.  \nIn 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 ≲ 3). However, at larger n ( ≳ 4), the scatter increases. This tendency is likely because larger n implies a compact S'ersic profile that resembles the PSF, making it challenging to distinguish from the PSF. We also find a ten-  \nFigure 7. Summary of the decomposition results in each filter. Columns (a), (b), and (c) respectively show the distribution of the estimated n , log r e in the unit of pixel and kpc. In column (d), the relation between H/T and S/N host is shown. In each panel, red (gray) indicates detection (non-detection) in each filter.  \n<!-- image -->  \nthat 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 .  \nRegarding 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 than the other two PSFs. Zhuang & 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.  \nIn 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.  \nFigure 8. Comparison of the decomposition results ( r e in the unit of arcsec, 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 (Sec. 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 ρ 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.  \n<!-- image -->", '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 RGB) of each host galaxy created by subtracting the PSF and the nearby Sersic component.  \nFirstly, 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 CID307 ( z = 2 . 05), despite their redshifts ∼ 2, we can clearly  \nidentify blue spiral arms with an overall diffuse red broad component.  \nAdditionally, galaxies such as CID-54 ( z = 0 . 97), CID510 ( 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 dustobscured region (e.g., Ito et al. 2024). Furthermore, there are cases with extended red structures (e.g., CID668; z = 0 . 97) which may indicate the presence of dust lanes as seen in the X-ray obscured (type-II) AGNs (Silverman et al. 2023).  \n4.4. Estimated M ∗ and the comparison of different SED fitting methods  \nFigure 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.  \n<!-- image -->  \nAs 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 aver-  \nF277W photometry-toM ∗ ratio in the sample. The assumed H/T of 0 . 2 is because galaxies with H/T ≳ 0 . 2 are almost detected in our methods (Figure 7(d)). The inferred M ∗ for each host galaxy is reported in Table 3.  \nFigure 10 shows four examples of the SED fitting results with the residuals. Generally, for the objects with z ≳ 1, the photometry or upper limit from F814W and F115W fall at a rest-frame wavelength shorter than 4000 ˚ A break and are important in constraining stellar population parameters of host galaxies.  \nWe 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 ∆ log M ∗ = +0 . 13 dex (corresponding to 1 σ ). This offset is not far from a common systematic M ∗ uncertainties among SED fitting methods reported in 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., Zou et al. 2019; Ishino et al. 2020; Shen et al. 2020; Li et al. 2021; Koutoulidis et al. 2022; Zhuang et al. 2023; Li et al. 2024).  \nFinally, 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.  \nWe compare our CIGALE M ∗ with 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 ∆ log M ∗ = +0 . 08 +0 . 20 -0 . 18 . This is very encouraging since their decomposition analysis is independent of our effort.  \nIn Figures 11 (b)-(d), we compare results for the other output SED model parameters: A V , t age , and τ . 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 τ , it is evident that significantly inconsistent values are observed around log ( τ/ Gyr) ∼ 0 . 5. Considering that t age is generally log ( t age / Gyr) ≲ 0 . 6 in our sample, this discrepancy can be attributed to the challenge of accurately determining SFH when τ ≫ t age .  \nAdditionally, we confirmed a strong negative correlation ( ρ = -0 . 60, p ≪ 0 . 05) between ∆log τ and ∆log M ∗ . This relation is because larger τ indicates a longer-lasting SFH, i.e., star formation has been persisting more recently. Consequently, galaxies with larger  \nτ 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 nonBayesian), 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 2 of ρ = 0 . 25. As done in previous local studies, we attempt to fit the observational data with a linear function of  \nlog M BH = α log M ∗ + β. (10)  \nThe red line represents the results with α ≃ 0 . 42 +0 . 12 -0 . 08 and β ≃ 4 . 08 +0 . 86 -1 . 26 , and the orange-shaded region showing the 1 σ 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 α of 0 and compared the BIC values. Although the BIC for the α -free model is lower than the BIC for the α = 0 model, the difference in BIC is less than 10. This means the α -free model is not significantly better than the α = 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 does 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.  \nFor comparison, we also provide the linear fit to the local sample consisting of 30 inactive galaxies (Haring & Rix 2004) and 25 active galaxies (Bennert et al. 2011b). The fit to these 55 galaxies results in α local ≃ 0 . 97 +0 . 10 -0 . 11 and β local ≃ -2 . 48 +1 . 15 -1 . 11 . 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 highz sample have different selection effects. Thus, we cannot directly compare α and β with α local and β local (see Sections 5.1 and 5.2 for the discussion with consideration of selection bias).  \nFor investigating the redshift dependence of the mass relation, we calculate ∆ log ( M BH /M ⊙ ), the relative off-  \nFigure 10. Representative examples of SED fits to PSF-subtracted host galaxies. Source ID and their estimated M ∗ with 1 σ confidence range are shown in the top right corner of each plot. The inverted triangle indicates the 3 σ upper limit due to the non-detection 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.  \n<!-- image -->  \nFigure 11. Comparison of SED fitting results with two different SED fitting codes. Panels (a) to (d) compare the estimated log M ∗ , A V , t age , and log τ . X and Y axis 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.  \n<!-- image -->  \n∗  \n/circledot  \nset of the black hole mass at given M ∗ from the local relation:  \n∆log( M BH /M ⊙ ) = M BH -α local log ( M ∗ /M ⊙ ) -β local , (11)  \nFor α local and β local , we use the above values from the local samples (Haring & Rix 2004; Bennert et al. 2011b), i.e., 0 . 97 and -2 . 48, respectively. We plot ∆log( M BH /M ⊙ ) as a function of z in the right panel of Figure 12. We then parameterize ∆log ( M BH /M ⊙ )  \nTable 3. Data summary of our sample.  \nz  \nFigure 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 σ confidence range. The best-fit relation and Spearman's correlation coefficient are reported at the top of the panel. ( right ) ∆log ( M BH /M ⊙ ) 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 lowz (Haring & Rix 2004; Bennert et al. 2011b), green for intermediatez (Jahnke et al. 2009; Bennert et al. 2011a; Cisternas et al. 2011; Schramm & Silverman 2013) and orange indicating those from Ding et al. (2020) at z ∼ 1 . 5. Red triangles indicate the M ∗ upper limit and ∆ log ( M BH /M ⊙ ) lower limit for five undetected targets.  \n<!-- image -->  \n<!-- image -->  \n∗  \n/circledot  \nto evolve with z as:  \n∆log( M BH /M ⊙ ) = γ log (1 + z ) . (12)  \nHere, we assume there is no redshift change in α and β and ∆log M BH can be described by only the evolution from the local relation (log M BH = α local log M ∗ + β local ). Fitting our data with Equation (12) without considering the selection bias results in γ = 1 . 33 +0 . 13 -0 . 14 , suggesting  \npositive evolution. Although R 2 for this fitting is 0.050, the BIC for this fitting is significantly lower than the γ = 0 model (i.e., ∆ log ( M BH /M ⊙ ) = 0 with a scatter); thus, this model describes the data better than the γ = 0 model. However, because our sample is (X-ray) flux-limited, it's essential to consider the impact of selection bias when determining the mass relation (Ding et al. 2020; Li et al. 2021; Li et al. 2021, e.g.,). Thus, for  \n∗  \nthe rest of this section, we measure the intrinsic slope of the mass relation ( α ) at z ∼ 1 . 5 and the redshift evolution parameter ( γ ) 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 ( α ) of the mass relation at z ∼ 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., Ding et al. 2020; 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.  \nIn particular, we initially assumed that α and β , for constructing the mock samples (Equation (9), have constant values independent of the redshift, and the evolution is expressed simply as γ log (1 + z ). However, the results in the previous Section 5 may indicate a different α from the local value without considering selection biases. In this section, we examine whether the constant α assumption is valid by estimating the intrinsic value of α while considering selection bias and provide evidence for an intrinsic relation between M BH and M ∗ at z > 1 for the first time.  \nSince the parameters α , β , and γ exhibit degeneracies, in this subsection, we constrain the z range to z = 1 -2 and perform fitting for α and β with an assumption of γ = 0 over this redshift range. Consequently, we cannot compare the estimated β with the values in the local relation directly. For the intrinsic dispersion of the mass relation σ µ , we assume 0.3 dex in this subsection, thus matching local studies.  \nWith 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 constrain the evolution parameters. We assume the observation thresholds as;  \nF [2 -10 keV] > F [2 -10 keV] , lim , (13)  \nFWHM line > 1000 km s -1 , (14)  \nM ∗ > 10 10 M ⊙ (15)  \nwhere the first and second conditions correspond to the detection limit of the X-ray observation and broad lines for single epoch M BH estimation. The third condition accounts for the detection limit of the host galaxies, and we confirm that galaxies with masses greater than 10 10 M ⊙ are successfully detected across all redshift ranges in our decomposition analysis. Thus, in comparing mock and real observations, galaxies with  \nmasses below 10 10 M ⊙ 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 [2 -10 keV] above the XMM-Newton F [2 -10 keV] , lim . Note that all of the eight targets with F [2 -10 keV] smaller than XMM-Newton F [2 -10 keV] , lim are targets observed only with Chandra. We changed F [2 -10 keV] , lim depending on in which survey each target was detected. Because M BH 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 H α , FWHM H β , and FWHM MgII .  \nFor 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 BH . Furthermore, we select the mock galaxies with a similar redshift; | z -z mock | < 0 . 1. Then, we calculate the probability that the mock galaxies with the similar z mock as the real galaxy would have the same ˜ M ∗ , mock with | ∆ M ∗ | < 0 . 1 and M vir , mock with | ∆ M BH | < 0 . 1, i.e., calculate the probability p following  \np = N | ∆ z | < 0 . 1 , | ∆ M ∗ | < 0 . 1 , | ∆ M BH | < 0 . 1 N | ∆ z | < 0 . 1 . (16)  \nWe 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 β .  \nThe obtained α -β distribution is shown on the left panel of Figure 13, indicating a strong anti-correlated degeneracy between α and β . Because the 1 σ contour includes values ( α local , β local ) 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 BH -M ∗ distribution for mock observed galaxies (similar to Figure 6f) 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 exhibits 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 BH range in future studies is essential. Even so, we demonstrate that a relation between M BH and M ∗ at highz 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.  \n<!-- image -->  \n∗  \n/circledot  \n∗  \n/circledot  \n∗  \n/circledot  \nFigure 13. ( left ) Posterior distribution in the α -β plane obtained through a comparison between mock observation and our sample. The posterior distributions for α and β are shown in the right and upper 1D histogram. Each contour indicates 1, 2, and 3 σ levels. The median and 1 σ confidence levels of α and β are shown in the top-right corner of each panel. ( right ) Comparison of mock observed masses generated with each α and β 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, 2, and 3 σ levels. Each parameter set, α and β , is shown in the upper left corner of each panel. As shown, we cannot reject the scenario that α does not evolve up to z ∼ 2 . 5.Figure 14. Posterior distribution on the σ µ -γ 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 σ µ and γ are shown in the right and upper 1D histogram. Each contour (red to orange) indicates 1, 2, and 3 σ from inner to outer. The median and 1 σ confidence levels of γ and σ µ are shown in the top-right corner of each panel. Only in panel (b), we plot the contour from Li et al. (2021) in gray, and their γ and σ µ 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.  \n<!-- image -->  \nTherefore, the posterior distribution of α and β that are estimated based on M ∗ also shows no significant PSF dependency.', '5.2. Intrinsic evolution ( γ ) of the M BH /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  \nparameters of γ (evolution rate) and σ µ (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 σ µ = 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 β  \nFigure 14 shows the estimated posterior distributions of γ and σ µ using each prior setting with the top-5 stacked PSF results. As evident in all panels by the orange contours, the intrinsic dispersion σ µ is strongly degenerate with the evolution rate γ where a smaller γ results in a larger σ µ as demonstrated in previous studies (Ding et al. 2020; Li et al. 2021). To reiterate, a smaller value for γ biases the mass relation towards relatively lower M BH values thus a larger σ µ is required to reproduce a certain set of observation data.  \n<!-- image -->  \n∗  \n/circledot  \n∗  \n/circledot  \nFigure 15. ( left ) Comparison of the observed sample and the mock observation of z ∼ 1 . 5 AGN on the M BH -M ∗ plane. The mock data at z ∼ 1 . 5 is constructed with the best-fit γ and σ µ (the result with the top-5 stacked PSFs assuming the 'Gaussian' prior) while assuming α local and β 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, 2, and 3 σ from inner to outer. ( right ) Comparison of the intrinsic M BH -M ∗ relation at z ∼ 1 . 5 (the red line) to the local relation (the gray line).  \nin Equation (9) to the values in the local relation ( α local , β local ) 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.  \nThere is still a possibility that σ µ depends on redshift. Nevertheless, as discussed later, imposing strong constraints on σ µ in our results is challenging due to the sample size and its uncertainties. Therefore, we set σ µ as a constant independent of redshift in the fitting. It means σ µ obtained through this method is considered to be an averaged value over z ∼ 0 . 68 -2 . 5. Even so, this σ µ estimation has the highest statistical significance for such a study at z ≳ 1. We discuss the redshift evolution of σ µ in 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.  \nIn the fitting process, we assume a uniform prior distribution for γ between -1 < γ < 1. Then, we have three different prior settings for σ µ : a uniform distribution between 0 . 01 < σ µ < 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 < σ µ < 1 . 0 with a prohibition of σ µ < 0 . 25 (realistic).  \nConsidering the likelihood distribution for γ , the 'basic' prior setting shows a slight positive-to-no evolution with γ = 0 . 46 +0 . 32 -0 . 63 . Similarly, the 'Gaussian' prior setting results in γ = 0 . 48 +0 . 31 -0 . 62 . If we assume that the intrinsic dispersion should not be significantly smaller than the local dispersion, we limit the allowed range for σ µ to be above 0 . 25 ('Realistic' case). In this case, we find γ = 0 . 22 +0 . 39 -0 . 58 , closer to the case for no evolution ( γ = 0) than the results with the other priors. For the latter, the intrinsic dispersion is slightly higher at 0 . 38 +0 . 12 -0 . 09 . In all cases, our results are consistent with very mild or essentially a lack of evolution with respect to the local relation.  \nThen, the left panel of Figure 15 compares mock observations using the median parameters ( γ = 0 . 48 and σ µ = 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  \ncompares 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 γ and σ µ from Section 5.2 with previous studies. Then, we also discuss the cosmic averaging scenario (Peng 2007; Hirschmann et al. 2010; Jahnke & Macci'o 2011).", '6.1. Comparison to other studies': "First, the conclusion of no or mild evolution with γ = 0 . 48 +0 . 31 -0 . 62 is consistent with studies based on 2D decomposition analysis (e.g., Ding et al. 2020; Li et al. 2021) and studies using a SED-fitting-based decomposition method (e.g., Sun et al. 2015; Suh et al. 2020).  \nIn Figure 14, we compare our results on the estimated γ -σ µ distribution to those of Li et al. (2021). Our sample has higher redshift range than Li et al. (2021), and the change of ∆ log ( M BH /M ⊙ ) is also proportional to log (1 + z ) in our model. Therefore, to reproduce the observational results, when increasing γ , we need to decrease σ µ more than Li et al. (2021). In other words, the slope of the γ -σ µ degenerate relation is steeper in our study. As a result, while the sample size of this study is approximately six times smaller than Li et al. (2021), the uncertainty of γ is only ∼ 2 times larger than Li et al. (2021). On the other hand, due to the steep slope, imposing constraints on σ µ becomes challenging, and the uncertainty becomes ∼ 4 times larger than Li et al. (2021). Even so, our estimated value of σ µ is 0 . 30 +0 . 14 -0 . 13 , which is remarkably similar to Li et al. (2021) with σ µ =0 . 25 +0 . 03 -0 . 04 .  \nIt is worth highlighting that the inference on the value of γ is very close to zero for the 'Realistic' case (Fig. 14c) where we assume that the intrinsic scatter ( σ µ ) cannot be lower than the local dispersion. Interestingly, if σ µ ( ∼ 0 . 4 -0 . 5) is actually higher than the local value, this would push the evolution parameter to negative values ( γ ∼ -0 . 5), thus presenting a scenario where the black holes have to catch up to their host galaxies by a bit.", '6.2. Scatter ( σ µ ) evolution and cosmic averaging': "When assuming a non-casual cosmic averaging scenario (Jahnke & Macci'o 2011), major mergers average and equalize the mass ratio M BH /M ∗ through cosmic history. Thus, σ µ should increase towards high redshift.  \nTo test the cosmic averaging scenario with our data, we generate a mock sample (20,000 parameter sets of M BH and M ∗ ) at z ∼ 1 . 4, the median redshift of our observational sample. We generate M ∗ based on the SMF by Weaver et al. (2023), and calculate the M BH assuming the γ and σ 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 Rodriguez-Gomez et al. (2015) covering z ∼ 0 -1 . 4. We simulate the redshift evolution of M BH and M ∗ by summing them with those of merging partners. Then, we calculate σ µ at each redshift to trace the expected scatter evolution from the assumed conditions. In this simulation, we do not consider the accretion 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 ± 0 . 5 dex.  \nFigure 16 compares the simulated redshift evolution of σ µ with the results from this study and Li et al. (2021). When assuming the evolution only through major mergers, the growth within the 1 σ uncertainty range significantly encompasses the results of Li et al. (2021). Moreover, our median redshift evolution is consistent with the results of Li et al. (2021). Therefore, the σ µ difference between the results from this study and Li et al. (2021) could be interpreted as the major merger-based scatter evolution.  \nHowever, due to the large uncertainty in our results, we cannot draw any definitive conclusions regarding the redshift evolution of σ µ . Furthermore, as evident from Figure 16, our sample has a wide redshift range compared to Li et al. (2021). If σ µ 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 γ -σ µ plane tends to steepen toward higher redshifts, making it relatively challenging to impose constraints on σ µ . Future highz statistical studies will likely require samples of a similar size to Li et al. (2021) with N = 584, or even larger to address the redshift evolution of σ µ . Thus, it will be necessary to conduct comprehensive surveys of highz 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. Ding et al. 2022a; Ding et al. 2023; Stone et al. 2023; Yue et al. 2023; Harikane et al. 2023; Zhuang & Shen 2023; Zhuang et al. 2023; Stone et al. 2023). As discussed above or in the previous studies, the results of the AGN+host galaxy 2D decomposition  \nFigure 16. Evolution of σ µ at z ≲ 2 based on a simple simulation of the cosmic averaging scenario. The vertical red error bars represent the scatter in σ µ within the z range of our sample using the basic prior. Faint grey lines correspond to each sampled parameter set while the red line and the orange filled region represent the median and 1 σ confidence level of σ µ redshift evolution. The green data denotes the constraint from Li et al. (2021) which is consistent within the 1 σ range of our results.  \n<!-- image -->  \ndepend significantly on the PSF reconstruction. Especially, Zhuang & Shen (2023) discussed the effect of different PSF on decomposition results. Zhuang & Shen (2023) compared three PSF modeling methods ( SWarp , photutils , and PSFEx ), but they did not compare them with χ 2 ν based methods directly. Thus, in this paper, we summarize the comparison when using different PSF reconstruction methods.  \nIn this study, we compare three final PSFs: two obtained by χ 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 Zhuang & Shen (2023), could potentially introduce biases in estimating n , r e , 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 e , and overestimation of host fluxes. However, as shown in Figure 8, morphological parameters such as n , r e , H/T generally exhibit consistent relations. Besides, as shown in Figure 8 (d), different PSFs have less impact on M ∗ estimation than on H/T . This could be due to the fact that M ∗ is estimated from SED fitting (Section 3.5) using multi-  \nhotometry that averages the uncertainty in each band. If so, SED fitting with a smaller number of photometric bands (e.g., Ding et al. 2023; 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.  \nModeling 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 χ 2 ν -based methods, as depicted in Figure 17. This method also allows flexible modeling, considering PSF as a function of position or brightness.  \nHowever, 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.  \nχ 2 ν -based selection method: Selecting PSFs based on χ 2 ν 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, Yue et al. (2023) discussed the differences in broad-band 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χ 2 ν single PSFs with a matched shape.  \nHowever, using fewer PSF candidates, like the top5 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 Zhuang & Shen (2023), a lower χ 2 ν 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.  \nGiven 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  \nin 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 e , 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 highz ( z ∼ 0 . 68 -2 . 5) type-I AGNs using the COSMOS-Web (Casey et al. 2023) and PRIMER-COSMOS surveys to measure the black hole - stellar mass relation at highz . Our sample contains 107 targets that are X-rayselected, broad-line AGNs with single-epoch black hole mass estimates (log ( M BH /M ⊙ ) ∼ 6 . 9 -9 . 6) based on H β , H α , and Mg ii from previous spectroscopic surveys (e.g., Schulze et al. 2018).  \nBy utilizing HST/ACS + JWST/NIRCam imaging that covers the rest-frame optical to near-infrared, we obtained multi-band 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 χ 2 ν -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 morphological parameters such as n , r e , 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 investigations of AGN host galaxies, such as double-S'ersic model fitting (decomposition fitting with AGN, bulge, and disk components) or parametric/non-parametric substructure analysis.  \nWith 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 ∼ 0 . 68 -2 . 5 (Figure 12). There is a weak positive correlation between M BH and M ∗ with the correlation coefficient of ρ = 0 . 25 ( p = 0 . 010). We fit the mass relation by a simple (log-)linear relation of log ( M BH /M ⊙ ) = α log ( M ∗ /M ⊙ )+ β with consideration of selection biases and measurement uncertainties. Our results show that the slope of the mass relation at z ∼ 2 ( α = 0 . 89 +0 . 61 -0 . 41 ) is consistent with the local relation ( α local = 0 . 97 +0 . 10 -0 . 11 ) (Figure 13).  \nAssuming 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 σ µ while considering selection biases and uncertainties based comparisons to mock catalogs (Figure 14). Even though the estimated probability distribution shows strong degeneracy between γ and σ µ , we find no or mild evolution with γ = 0 . 48 +0 . 31 -0 . 62 . If we assume that σ µ is not smaller than the local value, we obtain γ = 0 . 22 +0 . 39 -0 . 58 , which is more consistent with the no- or mild-evolution scenario. The estimated γ -σ µ distribution is largely consistent with 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 σ µ is steeper than Li et al. (2021). Therefore, despite the sample being approximately six times smaller, the estimated γ uncertainty is just slightly larger than Li et al. (2021).  \nFurthermore, the estimated value of the intrinsic scatter is σ µ = 0 . 30 +0 . 14 -0 . 13 which is consistent with the local relation and the recent estimate by Li et al. (2021). We show that this value σ µ at high-z may not be in contradiction to a cosmic averaging scenario (Figure 16) as recently put forward by Li et al. (2021) and 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 σ µ are weak. Thus, a larger sample size at z ∼ 1 -3 is needed, especially at highz .  \nFuture 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 AGNhost 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 accu-  \nrate PSF reconstruction, and establishing a framework for evaluating the uncertainties in reconstructed PSFs.  \n- This work is based on observations made with the 1\n- NASA/ESA/CSA James Webb Space Telescope. The 2\n- data were obtained from the Mikulski Archive for Space 3\n- Telescopes at the Space Telescope Science Institute, 4\n- which is operated by the Association of Universities 5\n- for Research in Astronomy, Inc., under NASA contract 6\n- NAS 5-03127 for JWST. These observations are associ7\n- ated with program IDs 1727 and 1837, and the specific 8\n- observations analyzed can be accessed via DOI. Nu9\n- merical computations were in part carried out on the 10\n- iDark cluster, Kavli IPMU. This work was made possi11\n- ble by utilizing the CANDIDE cluster at the Institut 12\n- d'Astrophysique de Paris, which was funded through 13\n- grants from the PNCG, CNES, DIM-ACAV, and the 14\n- Cosmic Dawn Center and maintained by Stephane 15\n- Rouberol. We thank Marko Shuntov, Kei Ito, and 16\n- Mingyang Zhuang for giving us useful advice for recon17\n- structing the PSFs. We thank Mingyang Zhuang, Yue 18  \n19  \nShen, and Junyao Li for giving us their measurements  \n- for comparison. We thank the anonymous referee for 20\n- helpful feedback. Kavli IPMU is supported by World 21\n- Premier International Research Center Initiative (WPI), 22\n- MEXT, Japan. The Cosmic Dawn Center (DAWN) is 23\n- funded by the Danish National Research Foundation 24\n- under grant DNRF140. TT is supported by Forefront 25\n- Physics and Mathematics Program to Drive Transfor26\n- mation (FoPM), a World-leading Innovative Graduate 27\n- Study (WINGS) Program, the University of Tokyo. JS 28\n- is supported by JSPS KAKENHI (JP22H01262) and 29\n- the World Premier International Research Center Ini30\n- tiative (WPI), MEXT, Japan. This work was sup31\n- ported by JSPS Core-to-Core Program (grant number: 32\n- JPJSCCA20210003). BT acknowledges support from 33\n- the European Research Council (ERC) under the Eu34\n- ropean Union's Horizon 2020 research and innovation 35\n- program (grant agreement number 950533) and from the 36\n- Israel Science Foundation (grant number 1849/19). JR 37\n- is supported by JPL, which is run by Caltech for NASA. 38\n- GEM and SG acknowledge the Villum Fonden research 39\n- grant 13160 'Gas to stars, stars to dust: tracing star for40\n- mation across cosmic time,' grant 37440, 'The Hidden 41\n- Cosmos,' and the Cosmic Dawn Center of Excellence 42\n- funded by the Danish National Research Foundation un43\n- der the grant No. 140. 44  \nAPPENDIX", 'Facilities: JWST (NIRCam), HST (ACS)': 'Software: astropy (Robitaille et al. 2013; PriceWhelan et al. 2018; Collaboration et al. 2022), CIGALE (Boquien et al. 2019; Yang et al. 2022), FSPS (Conroy et al. 2009; Conroy & Gunn 2010), galight (Ding et al. 2020), lenstronomy (Birrer & Amara 2018; Birrer et al. 2021), matplotlib (Hunter 2007), wCorr (Bailey & Emad 2023), numpy (Harris et al. 2020), photutils (Bradley 2023), psfr (Birrer et al. in prep), Prospector (Leja et al. 2017; Johnson et al. 2021), jwst (Bushouse et al. 2023)  \nTable 4. The number of the undetected host galaxies for each condition, PSF and filterTable 5. The number of host-galaxy detections for each final PSF  \nTable 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 non-detection 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 .', '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 semi-major 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 ∼ 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 PSFEx and semi-automatic PSF selection by galight . In addition, each final PSF tends to have a higher b/a than individual single PSFs. Regardless of whether χ 2 ν -based stacking or empirical modeling is employed, considering that the final PSF is a more reasonable  \nPSF reconstruction than a single PSF, this suggests a potential bias towards 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. Zhuang & 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.  \nFigure 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 semi-major 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 (semi-major 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.  \n<!-- image -->", '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 (Birrer & Amara 2018; Birrer et al. 2021) to generate 2,000 NIRCam/F150W mock galaxy images. These are constructed as a superposition of a PSF and a smooth analytic S'ersic model based on randomly sampled parameters: host-to-total flux ratio ( H/T ), S'ersic index ( n ), effective radius ( r e ), semi-major axis orientation ( θ ), 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 × r e, lenstronomy .  \nFigure 19 (a-c) compares 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 ans ) ∼ -2. For log ( H/T ans ) ≲ -2, the estimated value log ( H/T 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 ans ) ≲ -2 is practically challenging using the 2D decomposition analysis with galight . On the other hand, for log ( H/T ans ) ≳ -2, log ( H/T est ) exhibits a consistent trend along the y = x relation, while slightly larger than the input values. This trend is consistent  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 18. Comparison of the estimated H/T in each filter using different PSFs (Columns (a) to (d)). The gray and red colors correspond to whether they were detected in each filter. Spearman's correlation coefficient ρ 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.  \n<!-- image -->  \nIn 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 [ ∼ -5 , 0], n in [1 , 6], r e in [ ∼ 0 . '' 05 , 1 . '' 6], b/a in [0 . 2 , 1 . 0], and θ in [0 , π ] (in radian). We assume the PSFEx PSF (Section 3.1.2) and that the center of host galaxy components are 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χ 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'ersic components, such as double S'ersic (disc+bulge) components, position offset between AGN and host, substructures (bars, spirals, and clumps), and nearby or merging objects.  \nwith Zhuang & Shen (2023), who reported that different PSFs tend to overestimate H/T . As shown in Figure 7a, 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 19b, 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 largern 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 19d), we confirm that, with log ( H/T est ) > -1, the values of H/T can be recovered independently over a wide r e range of ∼ 0 . '' 1-1 . '' 2, effectively covering the typical r e of the sample (Figure 7b).  \nIn this paper, we assume the PSF + single S'ersic model to fit the galaxy images. As discussed in section 3.2, single S'ersic 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'ersic component. Even so, we expand the above mock tests using a PSF + double S'ersic model to address the effect of a bulge component on the PSF + single S'ersic model fitting. One of the double S'ersic components corresponds to a disk, and we assume r e, disk = 0 . '' 35 and n disk = 1 . 0. The other S'ersic component corresponds to a bulge, and we assume n bulge = 4 . 0. r e, bulge , H/T , bulge-to-host flux ratio B/H , θ disk , θ bulge , b/a disk , and b/a bulge is randomly sampled. The specified range for each parameter is log( H/T ) in [ ∼ -3 , 0], r e, bulge in [ ∼ 0 . '' 06 , 0 . '' 2], b/a disk in [0 . 6 , 1 . 0], b/a bulge in [0 . 8 , 1 . 0], θ disk in [0 , π ] (in radian), and θ bulge in [0 , π ] (in radian). Using these models and randomly sampled parameters, we make mock data by following the same method with the above single S'ersic mock data. Then, we performed 2D decomposition with the lowχ 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 ) ∼ -2 is larger than for the single S'ersic mock data (Figure 19a). However, similar to the results for the single S'ersic mock data, H/T is well reconstructed for log ( H/T est ) > -1. Figure 20b and c indicate the B/H dependence of n est and log ( H/T ) residual. From figure 20b, as B/H ans increases and the system becomes more bulge-dominated, the estimated n approaches values from 1 (= n disk ) to 4 (= n bulge ), indicating that the single S'ersic model fits the bulge for a bulge-dominated system. However, from figure 20c, 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.", 'REFERENCES': 'Akins, H. B., Casey, C. M., Lambrides, E., et al. 2024,  \narXiv e-prints, arXiv:2406.10341,  \ndoi: 10.48550/arXiv.2406.10341  \nBailey, P., & Emad, A. 2023, wCorr: Weighted  \nCorrelations.  \nhttps:  \n//american-institutes-for-research.github.io/wCorr/  \nBaker, W. M., Tacchella, S., Johnson, B. D., et al. 2023.  \nhttp://arxiv.org/abs/2306.02472  \nBeifiori, A., Courteau, S., Corsini, E. M., & Zhu, Y. 2012, MNRAS, 419, 2497,  \ndoi: 10.1111/j.1365-2966.2011.19903.x  \nBennert, V. N., Auger, M. W., Treu, T., Woo, J.-H., & Malkan, M. A. 2011a, ApJ, 742, 107,  \ndoi: 10.1088/0004-637X/742/2/107  \n-. 2011b, ApJ, 726, 59, doi: 10.1088/0004-637X/726/2/59 Bennert, V. N., Treu, T., Ding, X., et al. 2021, ApJ, 921, 36, doi: 10.3847/1538-4357/ac151a  \nBertin, E. 2011, in Astronomical Society of the Pacific Conference Series, Vol. 442, Astronomical Data Analysis Software and Systems XX, ed. I. N. Evans,  \nA. Accomazzi, D. J. Mink, & A. H. Rots, 435 Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, doi: 10.1051/aas:1996164 Birrer, S., & Amara, A. 2018, Physics of the Dark Universe, 22, 189, doi: 10.1016/j.dark.2018.11.002 Birrer, S., Shajib, A., Gilman, D., et al. 2021, The Journal of Open Source Software, 6, 3283, doi: 10.21105/joss.03283  \nFigure 19. Comparison of input (X-axis) and estimated (Y-axis) morphological parameters of mock galaxies for the single S\'ersic 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.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nQans  \nBoquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103, doi: 10.1051/0004-6361/201834156  \nBouchet, P., Garc´ıa-Mar´ın, M., Lagage, P.-O., et al. 2015, The Mid-Infrared Instrument for the James Webb Space Telescope, III: MIRIM, The MIRI Imager  \nBradley, L. 2023, astropy/photutils: 1.8.0, 1.8.0, Zenodo, Zenodo, doi: 10.5281/zenodo.7946442  \nBradley, L., Sip˝ocz, B., Robitaille, T., et al. 2023, astropy/photutils: 1.8.0, 1.8.0, Zenodo,  \ndoi: 10.5281/zenodo.7946442  \nBrightman, M., Silverman, J. D., Mainieri, V., et al. 2013, MNRAS, 433, 2485, doi: 10.1093/mnras/stt920  \nBrusa, M., Civano, F., Comastri, A., et al. 2010, ApJ, 716, 348, doi: 10.1088/0004-637X/716/1/348  \n<!-- image -->  \nFigure 20. (a) Comparison of the input and estimated log ( H/T ) for the double S\'ersic mock data. (b) The estimated n as a function of the input B/H . Blue and orange dashed lines indicate n disk and n bulge . (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 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 for the box plots. As shown, we can accurately reconstruct H/T flux ratio independent of bulge dominance B/H .  \n<!-- image -->  \nBruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000, doi: 10.1046/j.1365-8711.2003.06897.x  \nBushouse, H., Eisenhamer, J., Dencheva, N., et al. 2023, JWST Calibration Pipeline, 1.10.0, Zenodo,  \ndoi: 10.5281/zenodo.7795697  \nCalzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ,  \n533, 682, doi: 10.1086/308692  \nCappelluti, N., Brusa, M., Hasinger, G., et al. 2009, A&A, 497, 635, doi: 10.1051/0004-6361/200810794  \nCasey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2023,  \nApJ, 954, 31, doi: 10.3847/1538-4357/acc2bc  \nCen, R. 2015, ApJL, 805, L9,  \ndoi: 10.1088/2041-8205/805/1/L9  \nChabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392  \nCisternas, M., Jahnke, K., Bongiorno, A., et al. 2011,  \nApJL, 741, L11, doi: 10.1088/2041-8205/741/1/L11  \nCivano, F., Elvis, M., Brusa, M., et al. 2012, ApJS, 201, 30, doi: 10.1088/0067-0049/201/2/30  \nCollaboration, T. A., Price-Whelan, A. M., Lim, P. L., et al. 2022, The Astrophysical Journal, 935, 167, doi: 10.3847/1538-4357/ac7c74  \nConroy, C., & Gunn, J. E. 2010, ApJ, 712, 833,  \ndoi: 10.1088/0004-637X/712/2/833  \nConroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486, doi: 10.1088/0004-637X/699/1/486  \nDeGraf, C., Di Matteo, T., Treu, T., et al. 2015, MNRAS, 454, 913, doi: 10.1093/mnras/stv2002  \nDi Matteo, T., Colberg, J., Springel, V., Hernquist, L., & Sijacki, D. 2008, ApJ, 676, 33, doi: 10.1086/524921 Ding, X., Silverman, J. D., & Onoue, M. 2022a.  \nhttp://arxiv.org/abs/2209.03359  \nDing, X., Treu, T., Silverman, J. D., et al. 2020, ApJ, 896, 159, doi: 10.3847/1538-4357/ab91be  \nDing, X., Silverman, J., Treu, T., et al. 2020, The  \nAstrophysical Journal, 888, 37,  \ndoi: 10.3847/1538-4357/ab5b90  \nDing, X., Silverman, J. D., Treu, T., et al. 2022b, The Astrophysical Journal, 933, 132,  \ndoi: 10.3847/1538-4357/ac714c  \nDing, X., Onoue, M., Silverman, J. D., et al. 2023, Nature, 621, 51, doi: 10.1038/s41586-023-06345-5  \nDuras, F., Bongiorno, A., Ricci, F., et al. 2020, A&A, 636, A73, doi: 10.1051/0004-6361/201936817  \nElvis, M., Civano, F., Vignali, C., et al. 2009, ApJS, 184, 158, doi: 10.1088/0067-0049/184/1/158  \nFabian, A. C. 2012, ARA&A, 50, 455,  \ndoi: 10.1146/annurev-astro-081811-125521  \nFerrarese, L., & Merritt, D. 2000, ApJL, 539, L9, doi: 10.1086/312838  \nGarilli, B., Le F`evre, O., Guzzo, L., et al. 2008, A&A, 486, 683, doi: 10.1051/0004-6361:20078878  \nGraham, A. W., Onken, C. A., Athanassoula, E., &  \nCombes, F. 2011, MNRAS, 412, 2211,  \ndoi: 10.1111/j.1365-2966.2010.18045.x  \nGreene, J. E., Labbe, I., Goulding, A. D., et al. 2023, arXiv e-prints, arXiv:2309.05714,  \ndoi: 10.48550/arXiv.2309.05714  \nG¨ultekin, K., Richstone, D. O., Gebhardt, K., et al. 2009, ApJ, 698, 198, doi: 10.1088/0004-637X/698/1/198  \nHabouzit, M., Li, Y., Somerville, R. S., et al. 2021, MNRAS, 503, 1940, doi: 10.1093/mnras/stab496  \n- Netzer, H., & Trakhtenbrot, B. 2007, ApJ, 654, 754, doi: 10.1086/509650  \nOke, J. B. 1974, ApJS, 27, 21, doi: 10.1086/190287  \nOno, Y., Harikane, Y., Ouchi, M., et al. 2022.  \nhttp://arxiv.org/abs/2208.13582  \n- Onoue, M., Inayoshi, K., Ding, X., et al. 2023, ApJL, 942, L17, doi: 10.3847/2041-8213/aca9d3\n- Pacifici, C., Iyer, K. G., Mobasher, B., et al. 2023, ApJ, 944, 141, doi: 10.3847/1538-4357/acacff  \nPeng, C. Y. 2007, ApJ, 671, 1098, doi: 10.1086/522774  \n- Peng, C. Y., Impey, C. D., Ho, L. C., Barton, E. J., & Rix, H.-W. 2006, ApJ, 640, 114, doi: 10.1086/499930  \nPerrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, ed. J. Oschmann, Jacobus M., M. Clampin, G. G. Fazio, & H. A. MacEwen, 91433X, doi: 10.1117/12.2056689  \nPerrin, M. D., Soummer, R., Elliott, E. M., Lallo, M. D., & Sivaramakrishnan, A. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8442, Space Telescopes and Instrumentation 2012: Optical, Infrared, and Millimeter Wave, ed. M. C. Clampin, G. G. Fazio, H. A. MacEwen, & J. Oschmann, Jacobus M., 84423D, doi: 10.1117/12.925230  \nPrice-Whelan, A. M., Sip"ocz, B. M., Gunther, H. M., et al.  \n2018, The Astronomical Journal, 156, 123,  \ndoi: 10.3847/1538-3881/aabc4f  \n- Reines, A. E., & Volonteri, M. 2015, ApJ, 813, 82, doi: 10.1088/0004-637X/813/2/82\n- Rieke, M. J., Kelly, D. M., Misselt, K., et al. 2023, Publications of the Astronomical Society of the Pacific, 135, doi: 10.1088/1538-3873/acac53\n- Robitaille, T. P., Tollerud, E. J., Greenfield, P., et al. 2013, Astronomy and Astrophysics, 558, doi: 10.1051/0004-6361/201322068  \nRodriguez-Gomez, V., Genel, S., Vogelsberger, M., et al.  \n2015, MNRAS, 449, 49, doi: 10.1093/mnras/stv264  \n- Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772  \nSchneider, D. P., Richards, G. T., Hall, P. B., et al. 2010,  \nAJ, 139, 2360, doi: 10.1088/0004-6256/139/6/2360  \nSchramm, M., & Silverman, J. D. 2013, ApJ, 767, 13,  \ndoi: 10.1088/0004-637X/767/1/13  \n- Schulze, A., & Wisotzki, L. 2011, A&A, 535, A87, doi: 10.1051/0004-6361/201117564\n- -. 2014, MNRAS, 438, 3422, doi: 10.1093/mnras/stt2457  \nWang, B., de Graaff, A., Davies, R. L., et al. 2024, arXiv e-prints, arXiv:2403.02304,  \ndoi: 10.48550/arXiv.2403.02304  \nWeaver, J. R., Davidzon, I., Toft, S., et al. 2022, arXiv e-prints, arXiv:2212.02512,  \ndoi: 10.48550/arXiv.2212.02512  \n- -. 2023, A&A, 677, A184,  \ndoi: 10.1051/0004-6361/202245581  \nYang, G., Boquien, M., Brandt, W. N., et al. 2022, ApJ, 927, 192, doi: 10.3847/1538-4357/ac4971  \n- Yue, M., Eilers, A.-C., Simcoe, R. A., et al. 2023, arXiv e-prints, arXiv:2309.04614,  \ndoi: 10.48550/arXiv.2309.04614  \n- Zhang, Y., Ouchi, M., Gebhardt, K., et al. 2023, ApJ, 948, 103, doi: 10.3847/1538-4357/acc2c2\n- Zhuang, M.-Y., & Ho, L. C. 2023, Nature Astronomy, 7, 1376, doi: 10.1038/s41550-023-02051-4  \nZhuang, M.-Y., Li, J., & Shen, Y. 2023, arXiv e-prints, arXiv:2309.03266, doi: 10.48550/arXiv.2309.03266  \nZhuang, M.-Y., & Shen, Y. 2023.  \nhttp://arxiv.org/abs/2304.13776  \nZou, F., Yang, G., Brandt, W. N., & Xue, Y. 2019, ApJ, 878, 11, doi: 10.3847/1538-4357/ab1eb1'}