We report JWSTNIRCam measurements of quasar host galaxy emissions and supermassive black hole SMBH masses for six quasars at 5.9 lt z lt 7.1 in the Emissionline galaxies and Intergalactic Gas in the Epoch of Reionization EIGER project. We obtain deep NIRCam imaging in the F115W F200W and F356W bands as well as F356W grism spectroscopy of the quasars. We use bright unsaturated stars to construct models of the pointspread functions PSFs and estimate the errors of these PSFs. We then measure or constrain the fluxes and morphology of the quasar host galaxies by fitting the quasar images as a point source plus an exponential disk. We successfully detect the host galaxies of three quasars which have hosttoquasarflux ratios of 15. Spectral energy distribution fitting suggests that these quasar host galaxies have stellar masses of M SUBSUB  10SUP10SUP M SUBSUB. For quasars with host galaxy nondetections we estimate the upper limits of their stellar masses. We use the grism spectra to measure the H line profile and the continuum luminosity then estimate the SMBH masses for the quasars. Our results indicate that the positive relation between SMBH masses and host galaxy stellar masses already exists at redshift z  6. The quasars in our sample show a high BHtostellarmass ratio of M SUBBHSUBM SUBSUB  0.15 which is about 2 dex higher than local relations. We find that selection effects only contribute partially to the high M SUBBHSUBM SUBSUB ratios of highredshift quasars. This result hints at a possible redshift evolution of the M SUBBHSUBM SUBSUB relation.
['Quasars', 'Supermassive black holes', '1319', '1663', 'Astrophysics - Astrophysics of Galaxies']
EIGER. V. Characterizing the Host Galaxies of Luminous Quasars at z  6
{'EIGER V. Characterizing the Host Galaxies of Luminous Quasars at z ≳ 6': 'Minghao Yue , 1 Anna-Christina Eilers , 1 Robert A. Simcoe , 1 Ruari Mackenzie , 2 Jorryt Matthee , 2, 3 Daichi Kashino , 4, 5 Rongmon Bordoloi , 6 Simon J. Lilly , 2 and Rohan P. Naidu 1, ∗  \n1 MIT Kavli Institute for Astrophysics and Space Research, 77 Massachusetts Ave., Cambridge, MA 02139, USA 2 Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 27, Zurich, 8093, Switzerland 3 Institute of Science and Technology Austria (IST Austria), Am Campus 1, Klosterneuburg, Austria 4 National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan 5 Institute for Advanced Research, Nagoya University, Nagoya 464-8601, Japan 6 Department of Physics, North Carolina State University, Raleigh, 27695, North Carolina, USA', 'ABSTRACT': 'We report JWST /NIRCam measurements of quasar host galaxy emissions and supermassive black hole (SMBH) masses for six quasars at 5 . 9 < z < 7 . 1 in the Emission-line galaxies and Intergalactic Gas in the Epoch of Reionization (EIGER) project. We obtain deep NIRCam imaging in the F115W, F200W, and F356W bands, as well as F356W grism spectroscopy of the quasars. We use bright unsaturated stars to construct models of the point spread function (PSF) and estimate the errors of these PSFs. We then measure or constrain the fluxes and morphology of the quasar host galaxies by fitting the quasar images as a point source plus an exponential disk. We successfully detect the host galaxy of three quasars, which have host-to-quasar flux ratios of ∼ 1% -5%. Spectral Energy Distribution (SED) fitting suggests that these quasar host galaxies have stellar masses of M ∗ ≳ 10 10 M ⊙ . For quasars with host galaxy non-detections, we estimate the upper limits of their stellar masses. We use the grism spectra to measure the H β line profile and the continuum luminosity, then estimate the SMBH masses for the quasars. Our results indicate that the positive relation between SMBH masses and host galaxy stellar masses already exists at redshift z ≳ 6. The quasars in our sample show a high black hole to stellar mass ratio of M BH /M ∗ ∼ 0 . 15, which is about ∼ 2 dex higher than local relations. We find that selection effects only contribute partially to the high M BH /M ∗ ratios of high-redshift quasars. This result hints at a possible redshift evolution of the M BH -M ∗ relation.  \nKeywords: Quasars', '1. INTRODUCTION': "Supermassive black holes (SMBHs) are ubiquitously found in the centers of galaxies (e.g., Kormendy & Gebhardt 2001; Tremaine et al. 2002; Heckman & Best 2014). Observations of local galaxies have found tight correlations between SMBHs and the properties of their host galaxies, such as stellar masses or velocity dispersions, known as the M BH -M ∗ and the M BH -σ relations (e.g., Kormendy & Ho 2013). These relations suggest a strong co-evolution between SMBHs and their host galaxies, likely through feedback during the active  \nCorresponding author: Minghao Yue  \nmyue@mit.edu  \n- ∗ NHFP Hubble Fellow  \ngalactic nuclei (AGN) phases (e.g., Merloni & Heinz 2008; Ciotti & Ostriker 2001; Ciotti et al. 2010; Fiore et al. 2017). Specifically, AGN activities can produce strong outflows and expel the cold gas content in their host galaxies, quenching the star formation and also exhausting the gas supply to the SMBH (e.g., Fabian 2012; Cicone et al. 2014; King & Pounds 2015). Another important feedback mechanism is that AGNs can inject energy into the galaxy haloes through radio jets, which prevents halo gas from cooling and thereby shutting down both star-formation and BH accretion (e.g., Fabian 2012; Heckman & Best 2014).  \nMeanwhile, other mechanisms might also contribute to these observed relations. For example, SMBH accretion and galaxy star formation can be triggered by the same process (e.g., galaxy mergers; Di Matteo  \net al. 2005; Hopkins et al. 2008), producing a correlated growth of the SMBH and its host galaxy (e.g., Croton 2006; Storchi-Bergmann & Schnorr-Muller 2019). It has also been proposed that the SMBH-host correlation is a statistical effect resulting from the central limit theorem (e.g., Peng 2007; Jahnke & Macci'o 2011). Given the debate of the possible scenarios, the exact origins of the M BH -M ∗ and the M BH -σ are still unclear. One critical unresolved question is whether these correlations are already established in the early universe at very high redshifts or if they gradually take shape throughout cosmic time.  \nIn the past two decades, about 300 quasars at z > 6 have been discovered, which indicates that SMBHs with M BH ≳ 10 9 M ⊙ already exist when the universe is less than 1 Gyr old (e.g., Mortlock et al. 2011; Jiang et al. 2016; Matsuoka et al. 2018; Wang et al. 2019b; Yang et al. 2020; Wang et al. 2021; Ba˜nados et al. 2023). This quasar sample enables studies of SMBH-host coevolution in the early universe. By measuring the properties of these quasar host galaxies, we can characterize the M BH -M ∗ and the M BH -σ relations at z ≳ 6. So far, most of our knowledge about quasar host galaxies at z ≳ 6 comes from sub-millimeter (sub-mm) wavelengths. Recent observations with the Atacama Large Millimeter Array (ALMA) have shown that quasars at z ≳ 6 are hosted by massive starburst galaxies (e.g., Venemans et al. 2016; Decarli et al. 2018; Izumi et al. 2019; Yue et al. 2021). These quasar host galaxies have sizes of ∼ 2 -4 kpc in the sub-mm and show a variety of morphologies and kinematics, ranging from rotation-dominated disks to dispersion-dominated irregular mergers (e.g., Venemans et al. 2020; Neeleman et al. 2021).  \nNevertheless, one important missing piece in our knowledge is the stellar component of high-redshift quasar host galaxies. The M BH -M ∗ and M BH -σ relations of local galaxies describe the connection between the black holes and the stellar components of their host galaxies; however, ALMA observations trace the emission from cold dust and gas, making it hard to make direct comparisons between the above mentioned ALMA observations and the local relations. As luminous quasars are usually several magnitudes brighter than their host galaxies in rest-frame optical, probing the stellar emission of quasar host galaxies at z ≳ 6 is extremely challenging. One viable way to measure the emission from quasar host galaxies is image decomposition, utilizing the fact that quasars appear to be point sources and their host galaxies are extended. This approach requires sharp point spread functions (PSFs) to disentangle the flux from the quasar and its host galaxy.  \nAlthough image decomposition has been successful for quasars at z ≲ 2 (e.g., Mechtley et al. 2016; Ding et al. 2020; Chen et al. 2023), detecting quasar host galaxies at z ≳ 6 is much more challenging, given that the surface brightness of extended objects scales with redshift as (1 + z ) -4 . Even with the sharp PSF of the Hubble Space Telescope (HST) , previous studies have only reached non-detections of quasar host galaxies at z ≳ 6 (e.g., Marshall et al. 2020).  \nThis situation was completely changed, however, by the recent launch of the James Webb Space Telescope (JWST) . With its 6 . 4 -meter aperture and infrared coverage, JWST provides even sharper PSF than the HST and is sensitive to the rest-frame optical emission of high-redshift quasar host galaxies. An early study by Ding et al. (2023) reported the first detection of two quasar host galaxies at z > 6 using NIRCam imaging, indicating that the two quasar hosts are among the most massive galaxies at their redshifts ( M ∗ ≳ 10 10 M ⊙ ). Marshall et al. (2023) characterized the H β and [O iii ] emission lines of two quasar host galaxies at z ∼ 6 . 8 using NIRSpec IFU, showing complicated structures and kinematics of these galaxies. These exciting results motivate us to increase the sample of high-redshift quasars with host galaxy measurements in rest-frame optical, which will set a critical step towards fully understanding the co-evolution between SMBHs and galaxies in the early universe.  \nIn this paper, we report the measurement of the host galaxy emission and SMBH properties for six luminous quasars at z ≳ 6, using deep NIRCam imaging and spectroscopy as part of the Emission-line galaxies and Intergalactic Gas in the Epoch of Reionization (EIGER) project. We use the images to measure the fluxes and morphologies of the quasar host galaxies, and use the grism spectra to measure the black hole masses for the quasars from the H β emission line. Based on these measurements, we discuss the implication of the quasar host galaxies on the M BH -M ∗ relation in the reionization era.  \nThis paper is organized as follows. Section 2 describes the observations and data reduction. Section 3 describes the PSF modeling and the image fitting method we use to detect the rest-frame optical emission of the quasar host galaxies. We present the grism spectra and the SMBH mass measurements of these quasars in Section 4 and discuss the co-evolution between high-redshift SMBHs and their host galaxies in Section 5. We discuss our results in Section 6, and summarize this paper in Section 7. Throughout this paper, we assume a flat ΛCDM cosmology with Ω M = 0 . 3 and  \nH 0 = 70 km s -1 kpc -1 . All magnitudes are AB magnitudes unless further specified.", '2. OBSERVATIONS AND DATA REDUCTION': "The EIGER project (Proposal ID: 1243, PI: Lilly) is a Guaranteed Time Observation (GTO) program targeting six quasars at redshift 5 . 9 < z < 7 . 1, delivering deep NIRCam imaging and wide-field slitless spectroscopy (WFSS) of these quasar fields. The information about these quasars is summarized in Table 1. This Section briefly describes the observations, and we refer the readers to Kashino et al. (2023a), Matthee et al. (2023a), Eilers et al. (2023), and Bordoloi et al. (2023) for more information about the EIGER project.  \nWe obtain NIRCam F115W, F200W, F356W imaging, and F356W grism spectroscopy of the quasars. The observations of each quasar contain four individual visits, forming a mosaic that covers a field of view (FoV) of 3 ' × 6 ' around the quasar. See Kashino et al. (2023a) for more information about the mosaic configuration. The central 40 '' × 40 '' around each quasar is covered by every visits. We adopt the INTRAMODULEX primary dither pattern and the 4-Point subpixel dither pattern to improve the PSF sampling, cover detector gaps, and remove bad pixels. By the time of writing this paper, the observation of four quasars are completed, while we only have two of the four visits yet for J159-02 and three visits for J1030+0524. Our work presented here is based on these observations described above. The exposure time per visit is 4,381s for the F115W imaging, 5,959s for the F200W imaging, 1,578s for the F356W imaging, and 8,760s for the grism spectroscopy.", '2.1. NIRCam Imaging': "The NIRCam images were reduced using the jwst pipeline version 1.8.4. We first run Detector1Pipeline to generate the rate files (*rate.fits), then run Image2Pipeline to obtain calibrated images (*cal.fits). For astrometry, we first align the calibrated images to each other using tweakwcs , then combine all the images and calibrate the absolute astrometry to the Gaia DR2 catalog (Gaia Collaboration et al. 2018). We correct the 1 /f noise, mask snowballs, subtract the wisp patterns, and remove cosmic rays for the images using custom codes (see Kashino et al. 2023a, for more detailed description). We then run Image3Pipeline to stack images with the same visit and the same module. We use a pixel size of 0 . '' 03 for the F356W images and 0 . '' 015 for the F115W and the F200W images.  \nThere are some noticeable differences between the image reduction of this work and that of previous EIGER  \npapers. Previous EIGER papers presented JWST imaging by combining all exposures in each field to form a final co-added image per filter. In this work, we only combine images with the same filter, visit, and module. This approach reduces the systematic uncertainties in PSF modeling and image fitting (Section 3) introduced by imperfect astrometric alignment between the visits. Working on individual visits separately also allows us to estimate the systematic errors of the host galaxy measurements by comparing the results of different visits. In addition, previous EIGER papers used a pixel scale of 0 . '' 03 for the stacked F115W and F200W images, while in this work, we use a pixel scale of 0 . '' 015 to improve the sampling of the PSFs.", '2.2. NIRCam WFSS': "We obtain the NIRCam grism spectroscopy of the quasar fields using the grism 'R' in the F356W filter. This configuration gives an observed wavelength range of 3 . 1 µ m < λ < 4 . 0 µ m, which covers the H β emission line at 5 . 8 < z < 7 . 2. The data was reduced using jwst pipeline version 1.8.5. For each quasar, we first run the Detector1pipeline to obtain the *rate.fits files, then assign the world coordination system (WCS) information to the exposures using the AssignWcsStep and apply flat fielding using the FlatFieldStep . We remove 1 /f noise and sky background variations by subtracting the median value in each column. We then trace and extract the 2D spectra of the quasar from individual exposures using custom scripts utilizing the grismconf module (see Kashino et al. 2023b, for more details). Finally, we extract the spectra from all exposures using optimal extraction (Horne 1986), and combine the extracted 1D spectra using the coadd1d pipeline in the PypeIt package (Prochaska et al. 2020).", '3. MEASURING THE EMISSION OF QUASAR HOST GALAXIES': 'We use the NIRCam images to measure the flux and morphology of the quasar host galaxies. To do this, we fit the images of a quasar as a point source for the AGN plus an exponential disk for the host galaxy, and use the best-fit parameters of the exponential disk to infer properties of the quasar host galaxy.', '3.1. PSF Models and Errors': "We construct PSF models using bright and unsaturated stars in the NIRCam images. This method has been found to provide accurate PSF models for quasar host galaxy detection and outperforms webbpsf (e.g., Ding et al. 2022; Zhuang & Shen 2023). In this work, we  \nTable 1. The quasar sample of the EIGER project  \nFigure 1. The PSFs and the relative error maps, estimated from isolated bright stars in the images. These cutouts have sizes of 3 '' × 3 '' , and the integrated fluxes of the PSFs are normalized. The relative error is larger for brighter pixels, and the central bright pixels have relative errors up to ∼ 30%.  \n<!-- image -->  \nuse photutils 1 to build effective PSFs, which uses the algorithm described in Anderson & King (2000). The detailed approach is as follows.  \nFor each image, we first perform a source detection using the DAOphot algorithm (Stetson 1987), then select objects with magnitudes 18 < m < 21 in each filter and full-width half maxima (FWHMs) consistent with point sources 2 as PSF stars. The magnitude cut is determined to match the fluxes of the quasars and avoid saturation,  \nas the PSF shapes of infrared detectors exhibit flux dependence (e.g., the brighter-fatter effect; Plazas et al. 2018). We visually inspect all the PSF stars and reject those with close companions or bad pixels. We then use the EPSFBuilder class in the photutils package to build the empirical PSFs. Since the PSF of NIRCam depends on the position on the focal plane (e.g., Zhuang & Shen 2023), we construct the PSF models for module A and module B separately. We also notice that the number of suitable PSF stars in a single image is very limited (usually fewer than three), and some images do not have suitable PSF stars. We thus include PSF stars from all quasar fields and visits when fitting the PSF models. The typical number of PSF stars available for one filter and one module is ∼ 10 -20.  \nIt is worth noting that the PSF stars for the three bands are selected independently. Specifically, a PSF star in one band might be too bright or too faint to be selected as a PSF star in the other bands. As a result, the number of PSF stars available for the three bands are different.  \nThe empirical PSF models described above represent the average PSF of all the images. Limited by the number of PSF stars available, we are not able to model the spatial and temporal variations of the PSFs. Instead, we calculate the error maps of the PSF models, which estimate the possible differences between the PSF model and the real PSF of the quasar image.  \nSpecifically, we first compute the differences between the flux-normalized images of the PSF stars (denoted by P i ( x, y )) and the PSF model (denoted by ¯ P ( x, y )), then calculate the PSF error at pixel ( x, y ) as  \nϵ P ( x, y ) 2 = 1 N -1 N ∑ i [ P i ( x, y ) -¯ P ( x, y )] 2 -σ i ( x, y ) 2 (1)  \nwhere σ i ( x, y ) is the random noise (i.e., the ERR extension of the images) of the i -th PSF star. In practice, we perform sigma clipping with σ limit = 3 for each pixel in order to reduce the impact of outliers. The error map gives the standard deviation of the pixels in the PSF models and is added to the noise map in the image fitting (Section 3.2).  \nThe final product of the PSF modeling step contains six PSFs (three filters times two modules). The PSFs have sizes of 3 '' × 3 '' and have the same pixel size as the images. Figure 1 shows the PSF models and their relative error maps. There is a general trend that brighter pixels have larger relative errors; the central pixels of the PSFs have relative errors up to ∼ 30%. The mean relative errors of the PSFs within radii of 3 × FWHMare 0.204, 0.158, and 0.136 for the F115W, the F200W, and the F356W bands. This result is consistent with the finding in Zhuang & Shen (2023), who suggested that the errors of NIRCam PSF decrease toward long wavelengths.  \nAll PSF models, the PSF star lists, and the code to construct the effective PSF with error maps will be made available online upon publication.", '3.2. Image Fitting': "We use a point source component to describe the quasar and an exponential profile (i.e., a S'ersic profile with index n = 1) to describe the quasar host galaxy. We also add a constant background component to model imperfect background subtraction and add additional S'ersic profiles when there are other bright galaxies close  \nto the quasar. We use psfMC (Mechtley et al. 2016) to perform image fitting, which is a Python-based package explicitly designed for quasar host galaxy detections utilizing the Markov chain Monte Carlo (MCMC) method. We assign flat priors to all the free parameters, namely the magnitudes and positions of the point sources and S'ersic profiles, the half-light major and minor radii and the position angle of the S'ersic profiles, the S'ersic index of the S'ersic profiles (except for the quasar host galaxy which is forced to have n = 1), and the background level.  \nWe fix the S'ersic indices of the quasar host galaxies, which can only be poorly constrained due to the errors of the PSF models. The choice of n = 1 (i.e., exponential disks) is based on the following motivations: (1) previous detections of high-redshift quasar host galaxies are consistent with exponential profiles (Ding et al. 2023); (2) high-resolution ALMA observations found that the far-infrared emissions of quasar host galaxies have S'ersic index n ≈ 1 (e.g., Yue et al. 2021). We will discuss more about this point at the end of this Section.  \nIt is worth noting that psfMC incorporates PSF error maps when evaluating the likelihood of an image model. Specifically, the error of a pixel is calculated by combining the PSF error and the random noise of the image, i.e.,  \nϵ 2 all ( x, y ) = F 2 P ϵ 2 P ( x, y ) + σ 2 ( x, y ) (2)  \nwhere ϵ all ( x, y ) stands for the combined error at pixel ( x, y ), F P is the flux of the PSF, ϵ P ( x, y ) is the PSF error at pixel ( x, y ) (described in Section 3.1), and σ ( x, y ) is the random noise of the pixel. By including the PSF error map in the total error, we assign lower weights to pixels with large PSF uncertainties and reduce their impact on the fitting result. This combined error map also helps us to distinguish quasar host galaxy emissions from PSF inaccuracies in the PSF-subtracted images.  \nFor each quasar, we first fit its F356W images to determine whether the host galaxy is detected and to measure the morphology of the host galaxy. This choice is made based on two considerations. First, long-wavelength NIRCam filters exhibit smaller spatial variation, i.e., NIRCam PSFs are more stable at longer wavelengths, as we discuss in Section 3.1. Second, the flux ratio between the host galaxy and the quasar ( F G /F Q ) increases towards longer wavelength (Marshall et al. 2021). The reason is twofold: (1) quasars have blue continua in the rest-frame optical; (2) the F356W filter probes wavelengths longer than the 4000 ˚ A-break for z ∼ 6 -7 galaxies. As we will later show in this Section, the quasars with host galaxy detections in our sample have ( F G /F Q ) F356W ∼ 1 -3 × ( F G /F Q ) F200W .  \nWe fit the F356W images from individual visits separately, and determine the best-fit parameters and their  \nerrors by computing the median and standard deviation of the MCMC samples from all the visits. We also fit the images using a single PSF plus a sky background to evaluate the improvement in reduced χ 2 by including a host galaxy component in the model. Specifically, we use the relative improvement of the reduced χ 2 ,  \n∆ χ 2 ν χ 2 ν = χ 2 ν (1 p ) -χ 2 ν (1 p 1 e ) χ 2 ν (1 p 1 e ) (3)  \nwhere χ 2 ν (1 p ) is the reduced χ 2 for the single-PSF model, and χ 2 ν (1 p 1 e ) represents the reduced χ 2 where the exponential disk component is added to the model. A quasar host galaxy is considered to be detected in the F356W band if the fitting result meets the following criteria:  \n- 1. The error of the quasar host galaxy's magnitude is smaller than 0.3 (i.e., a detection with > 3 σ significance);\n- 2. The half-light radius satisfies R e, maj < 0 . '' 9 and R e, min > 0 . '' 1;\n- 3. The axis ratio satisfies 0 . 3 < R e, min /R e, maj < 1;\n- 4. All visits have ∆ χ 2 ν /χ 2 ν > 0 . 1.  \nHere, criterion (1) ensures that the host galaxy is consistently detected in all visits; criterion (2) ensures that the best-fit host galaxy has a reasonable size and is unlikely to be confused with the PSF component or the background component; criterion (3) further excludes some false detections caused by PSF inaccuracies; criterion (4) ensures that the quasar's image is significantly different from a single PSF. We will discuss why we adopt these criteria in Appendix A with more details.  \nFor quasars with successful host galaxy detections in the F356W band, we fit the F200W and the F115W images by fixing the position and the morphology of the host galaxy to the best-fit values from the F356W images. Again, we fit the images from different visits separately, and estimate the best-fit parameters and their errors using the median and the standard deviation of the MCMC samples from all the visits. We report a host galaxy detection in the F200W or the F115W band if the quasar has ∆ χ 2 ν /χ 2 ν > 0 . 07 and has a host galaxy magnitude error smaller than 0.3 mag. For quasars with non-detections in the F356W band, we fit their F200W and F115W images as a point source plus a constant background.  \nWenotice that the morphology and positions of quasar host galaxies might be wavelength-dependent. Fixing the morphology and positions of quasar host galaxies when fitting the F200W and the F115W band images might introduce additional systematic errors. However,  \nleaving these parameters free will return non-detections of quasar host galaxies in most cases, due to the fainter quasar host galaxies and the stronger PSF errors at short wavelengths. Other observations of high-redshift quasars find that the morphology of quasar host galaxies is consistent in the SW and LW bands (Ding et al. 2023). As such, we determine that fixing the host galaxy morphology and position to the best-fit F356W values is the best approach based on our data.  \nWe demonstrate how we validate the host galaxy detections and estimate the uncertainties of host galaxy properties in Figure 2. The left panel of Figure 2 shows the PSF-subtracted F356W image of J0148+0600 from individual visits. We emphasize that the quasar is located in module A (B) in visits 1 and 2 (3 and 4), and we model the PSF of module A and module B independently. As such, the similar patterns in the PSFsubtracted images from individual visits strongly indicate that the host galaxy detection is reliable. The right panel of Figure 2 shows the host galaxy fluxes measured from the four visits. The histograms represent the posterior distributions from MCMC, and the inter-visit differences reflect the systematic errors of image fitting. We estimate the uncertainties of the host galaxy magnitude by computing the standard deviation of MCMC samples from all the visits, which takes into account both random errors and systematic errors. In this case, the standard deviation of the host galaxy magnitude is ∆ m = 0 . 07, and the host galaxy is successfully detected.  \nOur image fitting procedure takes advantage of the multi-band, multi-visit observations of the EIGER project. The uncertainties of quasar host galaxy measurements are dominated by systematic errors of the PSF model instead of random photon noise. Fitting the four visits individually allows us to validate the result and estimate the systematic uncertainties by comparing the output of all the visits.  \nFigure 3 to 7 present the results of the image fitting. We summarize the best-fit parameters in Table 2. The PSF-subtracted images of these quasars exhibit a variety of features. We successfully detect the host galaxy of J0148+0600 in all the three bands, as well as J159-02 and J1120+0641 in the F356W and the F200W bands. J159-02 and J1120+0641 have ∆ χ 2 ν /χ 2 ν ∼ 0 . 05 in the F115W band, which do not satisfy the reduced χ 2 criterion; nevertheless, we still report tentative detections for their host galaxies in the F115W band because the magnitude errors are smaller than 0.3 mag. The flux ratios between the host galaxies and the quasars range from ∼ 1% -5%.", 'J0148+0600 F356W': "<!-- image -->  \nFigure 2. Fitting the F356W images of J0148+0600 from four visits. This Figure illustrates how we validate the detections of host galaxies. Left : the PSF-subtracted images, which show consistent shapes and brightness across the visits. Right : the MCMCposterior distribution of the host galaxy magnitudes. The small inter-visit differences reflect the systematic uncertainties of the host galaxy flux measurements. We take the median of all the MCMC samples from all four visits as the best-fit host galaxy magnitude, and estimate its uncertainty by computing the standard deviation of all the MCMC samples.  \n<!-- image -->  \nTable 2. The Results of Image Fitting  \nThe '(Image-PSF)/Error' maps clearly show the impact of the PSF error. Pixels that have large PSF errors (mainly those at central regions and on PSF spikes) are assigned lower weights in image fitting, and their influence on the fitting results is suppressed. The regions where the PSF errors are small exhibit significant host galaxy emissions compared to the noise level, confirming that the host galaxy signal is real and is not a result of PSF inaccuracies. We leave a more detailed discussion  \non how we validate the host galaxy detections and why we adopt the detection criteria in Appendix A.  \nThe host galaxies of J1030+0524 and J1148+5251 are not detected according to the detection criteria described above. However, the PSF-subtracted images of the two quasars clearly show extended emissions around the quasar, which might be tidal tails of a recent merger or the diffuse [O iii ] emission from galactic-scale outflows. The analysis of J0100 + 2802 has already been  \nreported in Eilers et al. (2023), and we refer the readers to Eilers et al. (2023) for more details. As a quick summary, J0100+2802 is saturated in all the images due to its extreme brightness, and we do not detect the host galaxy or any extended emission in its PSF-subtracted images.  \nWe also try to fit the images without fixing the S'ersic index of the host galaxy. We find that the S'ersic indices of the host galaxies are poorly constrained. Specifically, the posterior distribution of the indices spans a wide range at 0 . 3 ≲ n ≲ 4, with estimated errors σ n ≳ 1. We notice that the best-fit host galaxy magnitudes of the S'ersic profiles are similar to the exponential disk models (with differences of ≲ 0 . 3 mag). This systematic uncertainty is much smaller than the errors of the stellar mass estimates ( ∼ 0 . 3dex, see Section 3.3) and has no major influence on the main conclusions of this work.  \nWe provide more information about individual quasars in Section 6.1.", '3.3. SED fitting': "One of the main goals of this work is to characterize the position of luminous high-redshift quasars on the M BH -M ∗ plot. To measure the stellar masses of the quasar host galaxies, we perform SED fitting for these galaxies using Prospector (Johnson et al. 2021), with nebular emission treatments based on Cloudy (see Byler et al. 2017, for details). Since we only report tentative detections for the host galaxies of J159-02 and J1120+0641 in the F115W band, we only take the F356W and the F200W magnitudes as the input of the SED fitting for these two quasars.  \nWe assume a delayedτ model for the star formation history (SFH), i.e., SFR( t ) ∝ te -t/τ . We use a Chabrier initial mass function (Chabrier 2003) and assume a dust attenuation following the Calzetti et al. (2000) law. The free parameters and their priors of this SED model include: (1) the stellar mass M ∗ with a log-uniform prior at [10 8 M ⊙ , 10 12 M ⊙ ]; (2) the stellar metallicity log( Z/Z ⊙ ) with a uniform prior at [ -2 , 0 . 2]; (3) the starting time of the star formation t age with a uniform prior at [0 , t ( z )], where t ( z ) is the age of the universe at the quasar's redshift; (4) the exponential decay timescale τ with a uniform prior at [0 . 01Myr , 20Myr]; (5) the dust attenuation (quantified as the optical depth at 5500 ˚ A, τ 5500 ) with a uniform prior at [0 , 2]; (6) the gas-phase metallicity log( Z g /Z ⊙ ) with a uniform prior at [ -2 , 0 . 5]; and (7) the ionization parameter log U with a uniform prior at [ -3 , 1]. With only two or three photometric points for the SED fitting, we are only able to constrain the stellar masses, which roughly give the  \nnormalization of the SED. All other parameters remain largely unconstrained.  \nWe first perform the SED fitting for the detected quasar host galaxies. Figure 8 shows the best-fit SEDs, and Table 3 summarizes the physical parameters of the quasar host galaxies. The quasar host galaxies detected in this work have stellar masses of M ∗ ≳ 10 10 M ⊙ , which are among the most massive galaxies at z ≳ 6. We also plot the tentative F115W detections of the host galaxies of J159-02 and J1120+0641 for reference, but do not include them in the SED fitting. The F115W fluxes of J159-02 is consistent with the best-fit SED model, while the F115W flux of J1120+0641 is ∼ 2 σ lower than the SED model.  \nFor the quasars with only non-detections of their host galaxies, we estimate the upper limits of their host galaxy stellar masses. Specifically, we set conservative lower limits of the host galaxy magnitudes as m host = m QSO + 3 . 5 in the F356W band, then scale the best-fit host galaxy SED model of J0148+0600 to match this magnitude and take the corresponding stellar masses as upper limits. This choice is made because J0148+0600 has ∆ m = m host -m QSO = 3 . 5 in the F356W band and has a significant host galaxy detection. We note that the three quasars with host galaxy non-detections have F356W magnitudes similar to or brighter than J0148+0600; if the host galaxies of these quasars have ∆ m F356W = 3 . 5, we expect that the host galaxies should be detected at significance levels similar to or higher than J0148+0600. The stellar mass upper limits of these quasar host galaxies are also listed in Table 3, which have M ∗ ∼ 10 10 . 5 -10 11 . 5 M ⊙ .  \nSince we only have two or three photometric points available, we are not able to adopt more complicated star formation histories (e.g., a non-parametric star formation history) in the SED fitting. We also notice that the F356W fluxes of the quasar host galaxies contain both stellar continuum and the H β and [O iii ] nebulae lines. Since we do not have the nebulae line fluxes, we leave all emission-line-related parameters free to account for this systematic uncertainty. Future observations with NIRSpec IFU will provide the fluxes of these nebulae lines, which will allow us to use more complicated star formation histories and improve the accuracy of the estimated stellar masses (e.g., Marshall et al. 2023). It is also possible to measure the extended line emissions from quasar host galaxies using the two-dimensional grism spectra of the quasar. However, such measurements require careful analysis of the", 'J0148+0600': "Figure 3. The image fitting results of J0148+0600. From left to right: the original NIRCam image, the host galaxy model, the PSF-subtracted image, the residual image, the PSF-subtracted image divided by the composite error, and the residual image divided by the composite error. We obtain these images by combining images from all the visits. We detect the host galaxy in all three bands. The dashed ellipses mark the apertures corresponding to twice of the half-light radii, and the dashed circles mark the central regions with radii of 2 × FWHM PSF . The 'masked S/N' quoted here are the S/N of the host galaxy signal within the ellipse (excluding the central regions); see Appendix A for more details. J0148+0600 have two close companions, which are modeled as S'ersic profiles when fitting the images. Note that (1) the compisite error is a combination of random noises and PSF errors (see Equation 2); (2) the 'Host Model' images include the two close companions of the quasar.  \n<!-- image -->", 'J1120+0641': "Figure 4. Same as Figure 3, but for J1120+0641. We detect the host galaxy in the F356W and the F200W band, and report a tentative detection in the F115W band. The PSF-subtracted images show consistent shapes in all three bands. We note that we exclude the F356W image from visit 4 due to a large number of bad pixels around the quasar. The host galaxy of J1120+0641 is about 0 . '' 5 away from the quasar, and the PSF-subtracted images exhibit irregular shapes. These features suggest that J1120+0641 might be hosted by an on-going merger. See Section 6.1.5 for more details.  \n<!-- image -->  \nJ159-02  \nFigure 5. Same as Figure 3, but for J159-02. We detect the host galaxy in the F356W and the F200W band, and report a tentative detection in the F115W band. J159-02 has one close companion. Note that we only have two visits for J159-02 when writing this paper. See Section 6.1.4 for more details.  \n<!-- image -->", 'J1030+0524': 'Figure 6. The image fitting results of J1030+0524. From left to right: the original image, the PSF-subtracted residual, the residual normalized by the image error. We do not detect the host galaxy according to the criteria in Section 2.1. The PSFsubtracted images exhibit extended emissions around the quasar, which might be tidal tails of a recent galaxy merger. Note that we only have three visits for J1030+0524 when writing this paper.  \n<!-- image -->', 'J1148+5251': 'Figure 7. Same as Figure 6, but for J1148+5251. We do not detect the host galaxy according to the criteria in Section 2.1. The F356W image exhibits diffused emission extending from the lower left to the upper right corner, which is not seen in the other two bands. This feature might be [O iii ] emission around the quasar. We also see several close companions around the quasar that are detected in all three bands. Note that the central pixels in the F356W band are saturated.  \n<!-- image -->  \n5  \nFigure 8. Fitting the SEDs of the quasar host galaxies. The solid red dots represent the F200W and the F356W magnitudes of the quasar host galaxies. The black line shows the median modeled spectra generated using prospector , and the shaded area gives the 1 σ error. For comparison, we also show the F115W magnitudes of the host galaxies of J159-02 and J1120+0641 in open circles, which are not used when fitting their SEDs. The SED fitting enables constraints on the stellar masses of the quasar host galaxies.  \n<!-- image -->  \ntwo-dimensional grism data, which is beyond the scope of this paper.  \nFinally, we perform a sanity check of the stellar mass estimates by comparing our results with mock galaxies in UniverseMachine Data Release 1 (Behroozi et al. 2020). We select mock galaxies with redshifts 6 < z < 7 and UV magnitudes -23 < M UV < -21 . 5, which  \nroughly match the redshift and luminosity range of the quasar host galaxies. These mock galaxies have stellar masses of 10 9 . 2 M ⊙ -10 10 . 7 M ⊙ (95% confidence interval), with a median of 10 9 . 96 M ⊙ . These numbers are close to the stellar mass estimates for the quasar host galaxies in this work.  \nTable 3. The properties of the quasars and their host galaxies  \nNote - All errors are 1 -σ errors.', '4. BLACK HOLE MASS ESTIMATES FOR THE QUASARS': 'We calculate the black hole masses using the singleepoch virial estimator (e.g., Shen & Liu 2012). This method uses the FWHM of broad emission lines to estimate the velocity of the broad line region (BLR) clouds, and uses the continuum luminosity as a proxy of the distance from the BLR to the SMBH based on the luminosity-radius relation (e.g., Kaspi et al. 2005). Specifically, the SMBH mass can be calculated using the following relation:  \nM BH M ⊙ = 10 a ( FWHM 1 , 000 km s -1 ) b ( λL λ 10 44 erg s -1 ) c (4)  \nIn this work, we use the NIRCam grism spectroscopy to measure the FWHM of their H β emission lines and estimate the black hole masses accordingly. We use the empirical relation suggested by Vestergaard & Peterson (2006), which gives a = 6 . 91 , b = 2 , and c = 0 . 5.', '4.1. Spectral Fitting': 'We run MCMC to fit the spectra of the quasars to measure the profile of the H β lines. The flux model contains the following components: (1) a continuum described by a single power law; (2) the iron emission lines described using the template from Park et al. (2022); (3) the H β and [O iii ] λλ 4959 , 5007 emission lines. All the emission lines are fitted as two Gaussian components to model the complex broad line profiles seen in high-redshift quasars (e.g., Yang et al. 2023). The parameters of the spectral model include the amplitude  \nand power-law index of the continuum, as well as the fluxes, redshifts, and widths of the emission lines. We adopt flat priors for all these parameters. We further fix the flux ratio between the [O iii ]4959 and the [O iii ]5007 lines to be 1 : 3 (Storey & Zeippen 2000) and require that the two lines have the same redshifts and widths for both Gaussian components. The wavelength ranges to be fitted depend on the redshift of the quasars; specifically, we fit the window 3 . 20 µ m < λ < 3 . 85 µ m for J0148+0600, 3 . 40 µ m < λ < 4 . 05 µ m for J1120+0641, and 3 . 30 µ m < λ < 3 . 95 µ m for all the other quasars. For J1120+0641, the H γ line is redshifted to the wavelength window, and the [O iii ] lines fall on the edge of the transmission curve. We thus add the H γ line to the flux model and fix the [O iii ] redshift to ensure a successful fit.  \nFigure 9 shows the spectra and the best-fit flux models of the quasars. We measure the FWHM of the H β line and the continuum luminosity at rest-frame 5100 ˚ A ( L 5100 = λL λ (5100 ˚ A), and calculate the SMBH masses using Equation 4. We also estimate the bolometric luminosities of the quasars using L 5100 assuming a bolometric correction of 9.26 (Runnoe et al. 2012), and compute the Eddington ratios of the quasars. These results are listed in Table 3, and all the data will be available online after the paper is accepted for publication.  \nThe [O iii ] line profiles have been used to indicate possible outflows driven by the quasars (e.g., Yang et al. 2023). The quasars in the EIGER sample exhibit a variety of [O iii ] line profiles. J0100+2802, J0148+0600, J1030+0524, and J159-02 exhibit two broad [O iii ] com-  \nponents, while J1148+5251 only has one broad [O iii ] component. The [O iii ] lines of J1120+0641 are redshifted to the edge of the F356W filter and cannot be well-measured. This work focuses on the H β -based black hole masses, and we leave more detailed analysis of the [O iii ] emission lines to future studies.', '4.2. Comparison of Different Black Hole Mass Indicators': 'Previous studies have suggested that H β is a more reliable SMBH mass indicator compared to other broad emission lines (e.g., Shen & Liu 2012). However, the H β line of z ≳ 6 objects is not observable with ground-based facilities due to atmospheric absorption, and previous studies of z ≳ 6 quasars have been using the Mg ii line to measure the SMBH masses. The infrared coverage of JWST makes it possible to probe the rest-frame optical emission lines from high-redshift quasars. By analyzing the NIRCam grism spectroscopy of eight quasars at z > 6, Yang et al. (2023) showed that Mg ii -based SMBH masses are systematically higher than the H β -based SMBH masses. It is thus important to investigate the differences between the SMBH mass estimates of high-redshift quasars indicated by different emission lines.  \nFigure 10 shows the comparison between the Mg ii -based and the H β -based SMBH mass estimates for the six quasars in this work. We also include z > 6 quasars from Yang et al. (2023) and SDSS quasars from Wu & Shen (2022) for comparison. Given the systematic uncertainties of the BH mass estimators ( ∼ 0 . 3dex), the Mg ii -based and the H β -based SMBH masses agree with each other. The EIGER quasars have a mean value of log M BH, H β -log M BH, Mg ii = -0 . 093; this result is in line with Yang et al. (2023), who reported a mean value of log M BH, H β -log M BH, Mg ii = -0 . 13 for their z ≳ 6 . 5 quasar sample. These results suggest that the Mg ii -based BH masses of high-redshift quasars may be systematically larger by ∼ 0 . 1dex than their H β -based BH masses.', '5. THE SMBH-HOST GALAXY COEVOLUTION IN THE REIONIZATION ERA': "With the host galaxy stellar masses measured in Section 3 and the black hole masses measured in Section 4, we now characterize the position of z ≳ 6 luminous quasars on the M BH -M ∗ plot. Figure 11 illustrates the positions of the EIGER quasars on the M BH -M ∗ plot. We also include other high-redshift quasar host galaxies and local galaxies from the literature in this  \nFigure 3 . For the high-redshift quasar sample, there is a clear trend that more massive galaxies host larger black holes, indicating that the correlation between SMBHs and their host galaxies already exists at z ≳ 6. Meanwhile, the luminous quasars in the EIGER sample have M BH /M ∗ ∼ 0 . 15, which is ∼ 2 dex higher than the local M BH -M ∗ relation. The quasars with only nondetections of their host galaxies also lie above the local relation according to the upper limits of their stellar masses. This comparison suggests that SMBHs in luminous quasars might have experienced early growth compared to their host galaxies' star formation (e.g., Volonteri 2012).  \nSeveral previous studies have also used JWST NIRCam imaging to measure the stellar masses of highredshift quasars and AGNs. Stone et al. (2023a) reported NIRCam observations of a sub-Eddington quasar at z = 6 . 25, which also has an overmassive black hole compared to local galaxies. Ding et al. (2023) measured the host galaxy fluxes of two low-luminosity quasars at ∼ 6 . 3, which have small black holes ( M BH ∼ 10 8 -10 9 M ⊙ ) and massive host galaxies ( M ∗ ≳ 10 10 . 5 M ⊙ ). Ding et al. (2023) argued that the two quasar host galaxies are consistent with the local M BH -M ∗ relation after correcting for the selection effects.  \nRecent JWST observations have revealed a population of faint broad line AGNs at z ≳ 4 (also known as the 'little red dots', e.g., Kocevski et al. 2023; Matthee et al. 2023b; Maiolino et al. 2023). Harikane et al. (2023) analyzed a sample of low-luminosity AGNs at 4 < z < 7, which have relatively lower M BH /M ∗ ∼ 0 . 01. The broad-line AGN at z = 8 . 679 reported by Larson et al. (2023) have M BH /M ∗ ≈ 0 . 3%. However, Furtak et al. (2023) reported a low-luminosity AGN at z = 7 . 0451 with M BH /M ∗ ≳ 3%. These results show the large diversity of high-redshift SMBHs and their host galaxies, which might have experienced a variety of growth histories.  \nSome previous studies have used phenomenological approach to understand the redshift evolution of the M BH -M ∗ relation. For example, Caplar et al. (2018) developed an analytical method to derive the M BH -M ∗ relation that matches the observed star formation rate density and quasar luminosity function, who suggested that the BH-to-host mass ratio should be larger at higher redshift, following an evolution of M BH /M ∗ ∝ (1 + z ) 2 . 5 . Pacucci et al. (2023) showed that low-  \nFigure 9. The NIRCam grism spectra of the six quasars in our sample and the best-fit models. We fit the quasar spectra as a continuum power law (gray) plus iron emission lines (orange), the H β lines (blue), the [O iii ] lines (green), and the H γ lines (magenta) for J1120+0641. We use the widths of the H β lines to measure the mass of the SMBHs.  \n<!-- image -->  \nty AGNs at 4 < z < 7 recently discovered by JWST violate the local M BH -M ∗ relation by 3 σ . In contrast, Zhang et al. (2023a) showed that the M BH -M ∗ only evolves mildly at z < 10 using the empirical model Trinity , which was designed to match observables including quasar luminosity functions, quasar probability distribution functions, active black hole mass functions,  \nlocal SMBH mass-bulge mass relations, and the SMBH mass distributions of high-redshift bright quasars. This comparison again shows that it is not a trivial task to correctly characterize the high-redshift M BH -M ∗ relation.  \nThe stellar masses of the high-redshift quasar host galaxies still have large uncertainties. One of the main  \nFigure 10. Comparison between the Mg ii -based and the H β -based SMBH masses. We also include the quasars at z > 6 . 5 from Yang et al. (2023) and the SDSS DR16 quasars (Wu & Shen 2022) for comparison. Both the EIGER quasars and the quasars from Yang et al. (2023) have Mg ii -based BH masses that are larger than their H β -based BH masses. Note that the errorbars represent the random errors from spectral fitting; the systematic uncertainties of BH masses are ∼ 0 . 3dex.  \n<!-- image -->  \n/circledot  \nreasons is that we only have two or three photometric points for the SED fitting. By improving the PSF models for image fitting, it should be possible to improve the accuracy of the host galaxy flux measurements and include more photometric points for the SED fitting. Future observations with NIRSpec IFU will provide emission line fluxes (like the [O iii ] line) of the quasar host galaxies (e.g., Marshall et al. 2023), which will put stronger constraints on the stellar masses and other properties of the quasar hosts.", '5.1. Estimating the Selection Bias': 'It is worth noticing that selection effects can lead to apparent overmassive black holes in flux-limited samples of quasars (Lauer et al. 2007). Specifically, quasars with larger SMBHs are generally more luminous and are more likely to be detected in flux-limited surveys. Zhang et al. (2023a) showed that, for luminous AGNs with L bol > 10 46 erg s -1 and M ∗ ∼ 10 10 M ⊙ , this selection effect leads to a SMBH sample that is ∼ 1dex more massive than the true M BH -M ∗ relation. However, it is still unclear how to quantify this selection effect and unveil the intrinsic redshift evolution of the M BH -M ∗ relation. Following the method in Li et al. (2022), Ding  \net al. (2023) found that the two z ∼ 6 . 4 quasars in their sample are consistent with no redshift evolution of the M BH -M ∗ relation. In contrast, Stone et al. (2023a) argued that the large M BH /M ∗ values of high-redshift quasars can only be partially explained by selection effects.  \nThe EIGER quasars are among the most luminous quasars at z ≳ 6 (e.g., Fan et al. 2022) and are thus subject to significant selection bias. We follow the method described in Li et al. (2022) to estimate the impact of the selection bias on our results. To do this, we first generate a sample of mock galaxies by randomly drawing M ∗ values from the z ∼ 6 stellar mass function suggested by Song et al. (2016). We then assign a M BH to each mock galaxy following the low-redshift M BH -M ∗ relation in Zhuang & Ho (2023):  \nlog M BH = 1 . 49 log( M ∗ / 10 11 M ⊙ ) + 7 . 81 (5)  \nwith a scatter of 0.51 dex. We further assume that the Eddington ratios of these SMBHs follow the distribution of z ∼ 4 . 75 AGNs measured by Kelly & Shen (2013). We convert the bolometric luminosities of these SMBHs to L 5100 using the bolometric correction in Runnoe et al. (2012). In addition, we add random errors of 0 . 3 dex to M BH and 0 . 5 dex to M ∗ , which mimic the measurement errors of these properties (i.e., the scatters of the single-epoch Virial black hole mass estimator and the systematic error in SED fitting).  \nSince the quasars in our sample have L 5100 > 10 46 erg s -1 , we focus on mock galaxies in this luminosity range. The result is shown in Figure 12. The solid and filled contours represent low-redshift AGNs and mock high-redshift quasars, respectively. Due to selection bias, high-redshift mock quasars have M BH /M ∗ ratios that are ∼ 1dex higher than the local relation. Individual EIGER quasars locate outside the 68% confidency contour but fall inside the 95% confidency contour. Nevertheless, the M BH /M ∗ ratios of the entire EIGER sample are still ∼ 1dex higher than the median of the mock quasars. This comparison disfavors the hypothesis that the overmassive black holes seen in highredshift quasars are purely due to selection effects, and hints at a possible redshift evolution of the M BH /M ∗ ratio. This result is consistent with Pacucci et al. (2023), who previously inferred a significant redshift evolution of the M BH -M ∗ relation based on the high-redshift faint AGNs discovered by JWST .  \nThe high-redshift quasars in Stone et al. (2023b) have similar M BH /M ∗ ratios to the EIGER quasars, while the two z ∼ 6 . 3 quasars reported by Ding et al. (2023) are consistent with the mock quasars. However, we notice that the quasars in Ding et al. (2023) are fainter  \nFigure 11. The M BH -M ∗ relation. The filled and open red circles represent EIGER quasars with host galaxy detections and non-detections, respectively. For the EIGER quasars, we use a typical SMBH mass error of 0 . 3 dex and use stellar mass errors from the prospector models. We include high-redshift quasars from Ding et al. (2023) and Stone et al. (2023b), low-luminosity AGNs from Harikane et al. (2023) and Maiolino et al. (2023), and low-redshift galaxies from Kormendy & Ho (2013) and Reines & Volonteri (2015) for comparison. EIGER quasars have M BH /M ∗ ∼ 0 . 15, which is ∼ 2dex larger than the low-redshift relations. Note that the Kormendy & Ho (2013) relation uses bulge masses instead of total stellar masses and is located on the left side of the Reines & Volonteri (2015) relation.  \n<!-- image -->  \n∗  \n/circledot  \n( L 5100 ∼ 10 45 erg s -1 ) and thus have less severe selection biases compared to the luminous quasars in our sample.  \nFinally, we note that the EIGER quasars are not selected to be a flux-limited sample. The mock catalog method in this subsection can only provide a rough estimate of the selection effect. Accurate characterization of the selection bias requires a more uniformly selected sample and a more careful analysis (see, e.g., Pacucci et al. 2023; Zhang et al. 2023a,b, for relevant discussions), which is beyond the scope of this work. Figure 12 demonstrates the need for a larger sample of highredshift quasars with more accurate M BH and M ∗ measurements.', '6. DISCUSSION': '6.1. Notes on Individual Objects  \n6.1.1. J0100+2802  \nJ0100+2802 ( z = 6 . 31) is a hyperluminous quasar discovered by Wu et al. (2015). Due to its high luminosity,  \nJ0100+2802 is saturated in all three bands in EIGER imaging, and we are not able to detect its host galaxy. ALMA observations have suggested that the host galaxy of J0100+2802 is a starburst galaxy with a star formation rate of ∼ 850 M ⊙ / yr (Wang et al. 2019a), which exhibit clumpy structures (Fujimoto et al. 2020).  \nEilers et al. (2023) reported the NIRCam observations and ground-based spectroscopy of J0100+2802, who showed that the images of this quasar are welldescribed by a point source, and the C iv , Mg ii , and H β emission lines give consistent estimated BH masses. We note that the data used in this work was produced by a later version of the jwst pipeline compared to Eilers et al. (2023). The spectral fitting of this work is also slightly different from that of Eilers et al. (2023); we fix the line ratios and profiles of the [O iii ] doublet in the spectral fitting, while Eilers et al. (2023) left all emission line parameters free. Consequently, the spectral properties reported in this work are slightly different from that in Eilers et al. (2023).  \nFigure 12. Evaluating the selection bias of the luminous quasar sample. The black contours represent the distribution of low-redshift AGNs in Zhuang & Ho (2023). The filled contours represent the sample of mock quasars at z ∼ 6 with L 5100 > 10 46 erg s -1 , following the method described in Li et al. (2022). From inside to outside, the contours contain 68% and 95% of objects (corresponding to 1 σ and 2 σ levels). Due to the selection bias, the mock quasar sample lies above the low-redshift M BH -M ∗ relation. Nevertheless, the quasars in this work still have larger M BH /M ∗ compared to the median of the mock quasar sample.  \n<!-- image -->  \n∗  \n/circledot', '6.1.2. J0148+0600': 'J0148+0600 ( z = 5 . 977) was initially discovered by Jiang et al. (2015). This quasar is known for the large Gunn-Peterson trough in its spectruum (e.g., Becker et al. 2015) The PSF-subtracted images of J0148+0600 (Figure 3) indicate that its host galaxy has a regular elliptical shape. J0148+0600 has two projected companions in the north and the south directions, and it is unclear whether these two objects are associated with the quasar.', '6.1.3. J1030+0524': 'J1030+0524 ( z = 6 . 304) was initially reported by Fan et al. (2001). The image fitting process returns a nondetection of the quasar host galaxy according to the criteria in Section 3.2. Nevertheless, the PSF-subtracted images clearly show extended emissions around the quasar, which have consistent shapes in all three bands. These features look like tidal tails and suggest a recent galaxy merging event. Note that we only got the data from three visits for J1030+0524 by the time of writing this paper.  \nJ159-02 ( z = 6 . 381) was initially reported by Ba˜nados et al. (2016). The best-fit image model suggests a regular-shaped host galaxy with large radii. The PSFsubtracted image shows a projected companion galaxy. By the time of writing this paper, J159-02 has only been observed by two visits, where the quasar is located on module A.', '6.1.5. J1120+0641': "J1120+0641 ( z = 7 . 08) was initially reported in Mortlock et al. (2011). The PSF-subtracted images of J1120+0641 show irregular features, suggesting that the quasar might have experienced a recent galaxy merger. We also notice that the quasar is ∼ 0 . '' 5 away (about 3.9 kpc at z = 7 . 08) from the center of the host galaxy.  \nVenemans et al. (2017) reported the high-resolution ALMA observation of J1120+0641. The sub-mm continuum of J1120+0641 has a deconvolved FWHM of 1 . 24 kpc × 0 . 83 kpc, which is much smaller than the near-IR emission as revealed by NIRCam imaging. This result indicates that the sub-mm emission and the restframe optical emission come from different regions in the quasar host galaxy. Venemans et al. (2017) estimated the dynamical mass ( M dyn < 4 . 3 × 10 10 M ⊙ ), dust mass ( M dust ∼ (0 . 8 -4) × 10 8 M ⊙ ), and gas mass ( M gas ≲ 2 × 10 10 M ⊙ ) of the quasar host. These numbers are in line with the SMBH mass ( M BH = 1 . 2 × 10 9 M ⊙ ) and the stellar mass ( M ∗ = 5 . 6 × 10 9 M ⊙ ) estimated in this work.  \nWe notice that the Venemans et al. (2017) also reported an offset of ∼ 0 . '' 5 between the quasar (measured by ground-based imaging) and the host galaxy (measured by ALMA). Specifically, the sub-mm emission locates 0 . '' 22 west and 0 . '' 49 south from the optical quasar. This offset is roughly consistent with our result (see Figure 4). We note that, due to the limited spatial resolution and astrometric accuracy, Venemans et al. (2017) were not able to confirm this offset with sufficient significance.  \nBosman et al. (2023) reported the JWST Mid-Infrared Instrument (MIRI) spectrum of J1120+0641, which covers observed wavelengths 4 . 9 < λ obs < 27 . 9 µ m. Bosman et al. (2023) reported an H α -based BH mass of 1 . 5 × 10 9 M ⊙ , consistent with the H β -based results within scatters. Bosman et al. (2023) also found a redshifted H α core component (∆ v = -315 ± 37 km s -1 ), which was interpreted as a possible sign of a recoiling black hole by the authors. This scenario is in accordance with the offset between the quasar and the host galaxy emission.  \nJ1148+5251 ( z = 6 . 42) was initially reported by Fan et al. (2003). The image fitting procedure returns a nondetection of the host galaxy. Nevertheless, the PSFsubtracted image in the F356W band exhibit diffused emission extending from the northeast to the southwest of the quasar. We notice that this emission is absent in the F200W and the F115W images, which suggests that the emission might be dominated by H β or [O iii ] nebular lines.", '6.2. Systematic Errors and Possible Improvements of the Image Fitting Method': "Measuring the host galaxies of luminous quasars is a challenging task. Given the strong fluxes of the quasars, optimal PSF modeling and image fitting are needed to reveal the emission from the quasar host galaxies. In this Section, we discuss the PSF modeling and image fitting methods in previous studies and this work, as well as possible improvements to these methods.  \nWe first consider the method to build PSF models. Although recent JWST observations of quasar host galaxies have been using bright stars to construct PSF models, the exact methods adopted by these studies have some noticeable differences. For example, Ding et al. (2023) identified bright stars in the quasar's image, then chose the best five stars that give the smallest χ 2 in the image fitting as PSFs. Stone et al. (2023a) explicitly obtained images of a bright star and use the images as PSF, instead of using stars in the quasar's image. Zhuang & Shen (2023) tested three tools for PSF modeling, including SWarp (Bertin et al. 2002), PSFEx (Bertin 2011), and photutils . Zhuang & Shen (2023) suggested that PSFEx had the best performance, and that host galaxies with F G /F Q ∼ 10% can be securely detected with sufficient imaging depth.  \nIn this work, we use photutils to construct PSF models, where we gather stars from all quasar fields and all visits as PSF stars. The key feature of this work is that we estimate error maps of PSF models and add these errors into image fitting. We suggest that this step is critical and should be included in similar studies in the future. Specifically, if we only consider random noises when fitting the quasar images, the inaccuracies of the PSF models will have a substantial effect on the fitting result and will bias the estimated host galaxy fluxes. By including the PSF error maps in image fitting, we give lower weights to pixels that have larger PSF uncertainties and make the resulting host galaxy measurements less biased.  \nAnother systematic error is the SED mismatch between the quasars and the PSF stars. Since quasars and stars have different SEDs in the near-IR, the broad-  \nband PSF of these objects should have different shapes. To investigate this effect, we use webbpsf to generate two PSFs corresponding to two types of SEDs: (1) an M dwarf from the stellar spectrum library by Castelli & Kurucz (2003); (2) a quasar at z = 6 with the Vanden Berk et al. (2001) spectral template. For both the F356W and the F200W band, the differences between the two PSFs have similar levels to the PSF error maps described in Section 3.1. We thus expect that the impact of SED differences should be comparable to that of PSF error maps. The uncertainties introduced by the PSF error maps can be estimated by the MCMC samples from a single visit, which introduces an error of σ m, host ∼ 0 . 02 for the host galaxy magnitude measurements. This value is much smaller than the magnitude errors we report in Table 2, which are dominated by the inter-visit systematic uncertainties.  \nIt is still unclear what is the optimal way to model the PSF of NIRCam. In particular, the PSF of NIRCam depends on the position on the detector, the flux of the source, and the time of the observation. These effects are ignored in this work due to the limited number of PSF stars available, and we leave the detailed analysis of these effects to future studies.  \nWe now consider the method for image fitting. Most studies of quasar host galaxies (including this work) uses a point source plus a S'ersic profile to describe the image of quasars. The specific setting for the S'ersic profiles varies between the studies. Ding et al. (2023) left all S'ersic parameters free, Zhuang & Shen (2023) fixed the center of the host galaxy to the quasar's position, and this work fix the S'ersic index of the host galaxy to be n = 1 (i.e., an exponential disk). Studies based on ground-based imaging have used one-dimensional profile fitting, which is a powerful tool for low-redshift quasar host galaxies (e.g., Matsuoka et al. 2014; Yue et al. 2018).  \nThe result of this work demonstrates the limitation of the existing image fitting method. Several quasars in the EIGER sample show irregular emissions in their PSF-subtracted images. For these quasars, a regular S'ersic profile might not give the correct description of the host galaxy. It is thus desirable to develop image fitting methods that can describe irregular galaxy shapes. Additionally, we suggest that the position of the host galaxy should be left free instead of fixed to the quasar's position since some quasars show large offsets from their host galaxies (like J1120+0641).  \nUsing the bluetides simulation (Feng et al. 2016), Marshall et al. (2021) estimated that the detection limit of high-redshift quasar host galaxies is F G /F Q ∼ 1% for NIRCam (assuming an exposure time of 10 ks).  \nThis work successfully detects quasar host galaxies that have F G /F Q ∼ 1% with exposure times of ∼ 1 . 6 ks for the F356W imaging. With improved PSF models in the future, it is possible to achieve detection limits of F G /F Q ≲ 1%, allowing us to measure quasars with fainter host galaxies.", '7. CONCLUSION': 'We present NIRCam observations of six quasars at z ≳ 6 observed by the EIGER project. We use NIRCam imaging to measure the host galaxy emissions of the quasars, where we fit the quasar images as a point source plus an exponential disk. We construct PSF models and their error maps using bright stars in the images, and run MCMC to perform image fitting to estimate the fluxes of the quasar host galaxies. We use NIRCam grism spectra to measure the profile of the broad H β emission lines and calculate the SMBH masses of the quasars. The main results of this work are as follows:  \n- 1. We detect the host galaxy of J0148+0600 in all the three bands, as well as the host galaxies of J15902 and J1120+0641 in the F200W and the F356W bands. We report tentative detections for the host galaxies of J159-02 and J1120+0641 in the F115W band. These quasars have host-to-quasar flux ratios of ∼ 1% -5%. SED fitting shows that these quasar host galaxies have M ∗ ≳ 10 10 M ⊙ .\n- 2. We report non-detections for the host galaxies of J0100+2802, J1030+0524, and J1148+5251. We also estimate the upper limits of the fluxes and stellar masses of their host galaxies. The PSFsubtracted image of J1030+0524 and J1148+5251 show diffused emissions around the quasar, which might come from the tidal tails of the ongoing galaxy merger or extended [O iii ] emissions around the quasar.\n- 3. We compute the black hole masses of the six quasars using their H β emission lines. The H β -based black hole masses of these quasars are slightly smaller than the Mg ii -based ones by ∼ 0 . 1dex, which is consistent with the results by Yang et al. (2023).\n- 4. The quasars in the EIGER sample have M BH /M ∗ ∼ 0 . 15, which is ∼ 2 dex larger than low-  \nft galaxies with similar stellar masses. This comparison suggests that high-redshift quasars might have experienced early SMBH growth compared to the star formation of their host galaxies. Selection bias also contributes to the high M BH /M ∗ ratios of luminous quasars. However, the M BH /M ∗ ratios of EIGER quasars are larger than the mock luminous quasar sample even after the selection bias is considered, hinting at a possible redshift evolution of the M BH -M ∗ relation.  \nThe EIGER quasars are among the most luminous quasars at z ≳ 6. The detection of their host galaxies illustrates the promising potential for JWST to build a large sample of high-redshift quasars with host galaxy detections in the near-IR. The PSF models of NIRCam will be improved by future observations and analyses, which will enable more accurate characterization of the high-redshift M BH -M ∗ relation. Meanwhile, JWST observations have discovered several AGNs at even higher redshifts (e.g., Furtak et al. 2023; Goulding et al. 2023; Larson et al. 2023; Kokorev et al. 2023). Follow-up analysis of these AGNs will characterize the co-evolution of SMBHs and their host galaxies at even earlier cosmic times, approaching the origin of the M BH -M ∗ relation that already has its shape at z ∼ 6.  \nWe thank the referee for the valuable comments on this paper. We thank John Silverman, Madeline Marshall, Ming-Yang Zhuang, Weizhe Liu and Jinyi Yang for inspiring discussions and suggestions. DK thanks the support from JSPS KAKENHI Grant Number JP21K13956. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observa1 2 3 4 5 6 7 8 9 10 11 12  \n- tions are associated with program ID #1243. 13  \nFacilities: JWST(NIRCam)  \nSoftware: astropy (Astropy Collaboration et al. 2013, 2018), psfMC(Mechtley 2014), webbpsf(Perrin et al. 2014), jwst', 'APPENDIX': 'Figure 13. Testing the reduced χ 2 of the quasars and stars. The histograms present the distribution of ∆ χ 2 ν /χ 2 ν (Equation 3), and the black dashed line marks the criterion for successful host galaxy detection. The ∆ χ 2 ν /χ 2 ν of the quasars with host galaxy detection, averaged over all visits, are listed in red text. In the F356W and the F200W bands, all three quasars have ∆ χ 2 ν that are much higher than the PSF stars. In the F115W band, only J0148+0600 passes the detection criteria in the F115W band, and J159-02 and J1120+0641 have ∆ χ 2 ν that are similar to some stars.  \n<!-- image -->  \n-  \n-  \n-', 'A. VALIDATING THE QUASAR HOST GALAXY DETECTIONS': "Given the small flux ratios between the detected host galaxies and the quasars, we need to ensure that the detected quasar host galaxy emissions are real and are not results of PSF model inaccuracies. We perform this test by running the image fitting procedure on the PSF stars selected in Section 3.1. For each PSF star, we first build a new empirical PSF model by excluding the star from the PSF star list, then fit the star as a point source plus an exponential disk and a sky background. If the PSF models are accurate, we will get non-detections for the exponential disk components. In reality, we may get positive fluxes for the exponential disk components due to the inaccuracies of the PSF models, according to which we can examine whether the quasar host galaxy detections we report are reliable.  \nWe first examine the reduced χ 2 values for stars and quasars. The result is shown in Figure 13. Quasars have significantly higher ∆ χ 2 ν /χ 2 ν values than stars. Specifically, all stars have ∆ χ 2 ν /χ 2 ν < 0 . 1 in the F356W band and ∆ χ 2 ν /χ 2 ν < 0 . 07 in the F200W and the F115W band; in contrast, the three quasars with host galaxy detections have χ 2 ν /χ 2 ν > 0 . 1 in the F356W and F200W bands. The differences in ∆ χ 2 ν /χ 2 ν between stars and quasars are smaller in the F115W band. Figure 13 suggests that the host galaxy detections for the three quasars are real at least in the F356W and the F200W bands. These results also explain why we adopt the ∆ χ 2 ν /χ 2 ν criteria in Section 3.2.  \nAnother useful indicator of the quasar host galaxy detections is the signal-to-noise ratio ( S/N ). For each quasar, we evaluate the S/N of the host galaxy signal within an ellipse aperture that corresponds to two times the half-light radii. The dashed ellipses in Figure 3 to 5 mark these apertures. We notice that, due to the contribution of the PSF errors, the noises at central pixels are much larger than the noises at larger radii and can even be larger than the total flux of the host galaxy. Meanwhile, the '(Image-PSF)/Error' columns in Figure 3 to 5 show that the host galaxy signal is significant at larger radii where the PSF errors are smaller. Based on this consideration, when evaluating the S/N of quasar host galaxies, we exclude the central regions with r < 2 × FWHM PSF where the PSF errors are high. The dashed circles in Figure 3 to 5 represent the masked central regions. We then compute the signal ( S ) and noise ( N ) in the aperture:  \nS = ∑ ( x,y ) I ( x, y ) -P ( x, y ) , N = √ ∑ ( x,y ) ϵ all ( x, y ) 2 (A1)  \nwhere I ( x, y ) is the original image, P ( x, y ) is the best-fit PSF component, and ϵ all ( x, y ) is the composite noise map (Equation 2). The S/N of these quasar host galaxies are presented in Figure 3 to Figure 5; these quasar host galaxy detections have S/N > 10, which again confirms that the host galaxy signals cannot be explained by PSF errors.  \nWe further examine the best-fit parameters of the exponential disk component for the stars. The result of the experiment is shown in Figure 14. Again, we note that the three bands have distinct PSF star samples due to the magnitude requirements (Section 3.1); in particular, there is no overlap between the PSF stars for the F356W band and the F200W band. The open black circles mark stars where the best-fit exponential profile does not satisfy the radiirelated criteria (criteria (2) and (3) in Section 3.2), for which our image fitting method would report non-detections of the host galaxies. The solid circles mark stars that satisfy the radii-related criteria. Note that we do not apply  \nFigure 14. Testing the host galaxy detection reliability using field stars. We fit the stars as a PSF plus an exponential disk and a background component, following the same method we used for quasar image fitting. Open circles mark stars that do not satisfy the radii-related criteria, filled black circles mark stars that pass the radii-related detection criteria, and the red circles represent the quasars with host galaxy detections. Most stars in the F356W and the F200W band fail to pass the radii-related criteria, while some stars in the F115W band pass the radii-related detection criteria and have m host -m PSF similar to the real quasars. This Figure demonstrates the reason why we adopt the radii-related detection criteria for the quasar host galaxies.  \n<!-- image -->  \nm  \nhost  \n-  \nm  \nPSF  \ncriterion (1) to the stars because most stars were only observed by one visit, meaning that we cannot estimate the inter-visit systematic uncertainties for the stars.  \nFor the majority of stars, the best-fit exponential profiles have either a very large radii (mimicking a sky background), a very small radii (mimicking a point source), or a small axis ratio (mimicking the PSF residual, especially on the spikes). As a result, most ( > 90%) PSF stars in the F356W band and the F200W band fail to pass criteria (2) and (3) and thus return non-detections of the host galaxy component. This result explains why we adopt these criteria in Section 3.2. We also notice that stars in the F356W and F200W bands that pass criteria (2) and (3) have faint best-fit exponential profiles (∆ m = m exp -m PSF ≳ 5, where m exp stands for the best-fit magnitude of the exponential disk component). Meanwhile, some stars in the F115W band return an exponential disk component that is ∼ 3 . 5 -4 . 5 magnitudes fainter than the point source component. By visually inspecting their PSF-subtracted images, we find that these stars have strong PSF wings and spikes, which mimic the extended exponential disk. This result is also consistent with the finding in Section 3.1 and in Zhuang & Shen (2023), i.e., the PSF spatial variations are stronger in bluer bands.  \nThe red circles in Figure 14 represent the quasar host galaxies detected in Section 3.2. These quasars pass criteria (2) and (3) in the F356W band, indicating that the host galaxy detections are not results of PSF inaccuracies. Another piece of evidence for the host galaxy detections is the consistent results from individual visits. Figure 2 presents the images of J0148+0600 as an example. For the other two quasars, we also find that the PSF-subtracted images from individual visits show similar patterns and that the best-fit host galaxy magnitudes from different visits are consistent within ≲ 0 . 3 magnitudes.  \nWe notice that a few stars in the F115W band return ∆ m and R e values similar to the quasar host galaxies, meaning that it is possible that PSF inaccuracies contribute significantly to the F115W fluxes of the quasar host galaxies. This result is consistent with our finding in Section 3.2 that the PSF errors are larger in the F115W band. Furthermore, J1120+0641 and J159-02 fail to pass the ∆ χ 2 ν criteria in the F115W band. We thus exclude the F115W magnitudes in the SED fitting for these two quasar host galaxies.  \nThe next question is whether the reported quasar host galaxy properties (i.e., magnitudes and radii) and their errors are reliable. We answer this question using mock observations. For each quasar, we generate a set of mock images by adding an exponential profile to the PSF star images. We scale the images and the error maps of the PSF stars to match the total flux of the quasar, and use the best-fit parameters of the quasar host galaxy to generate the exponential profile. We notice that the errors of image fitting can be large for stars that are much fainter than the quasar. To mitigate this effect, we only use stars that have | m star -m quasar | < 1 to generate mock observations. The", 'Yue et al.': 'Figure 15. Testing the reliability of the reported quasar host galaxies and radii. See text for details about mock observation generation. Top panel: the F356W band test, where the morphology and positions of the exponential profile are left free. Each black dot represents one realization, and the red dot represents the real observations on the quasar (which is also the input value to generate mock observations). Middle and Bottom : the F200W and the F115W band test, where the morphology and positions of the exponential profile are fixed to the input values. The histograms represent the realizations of mock observations, and the red line and the shaded area represent the best-fit host galaxy magnitude and its error from real observations on quasars. The mock observations give exponential profile magnitudes that are consistent with the input value (except for J159-02 in the F115W band, which we report as a tentative detection). This test validates the quasar host galaxy magnitudes and their errors we report in this work.  \n<!-- image -->  \nonly exception is J0148+0600 in the F200W and F115W bands, where we adopt -1 < m star -m quasar < 2 to include sufficient numbers of stars.  \nWe then run image fitting for the mock images to examine whether we can recover the input host galaxy parameters. Following the method used for the quasars, we leave the exponential disk morphology free when fitting the F356W images, and keep the morphology fixed to the true value when fitting the F200W and the F115W images. We mimic the multiple visits to the quasars by randomly drawing mock observations. Specifically, for one realization, we randomly select N mock observations (where N is the number of visits), then compute the median and the standard deviation of the host galaxy parameters among these observations as the best-fit values and their errors. We generate 100 realizations for each quasar and each band.  \nFigure 15 presents the results of this test. In the top panel, we examine the best-fit exponential profile magnitudes and radii from mock F356W observations. Each gray point represents one realization, and the red point represents the real observation of the quasar (which is also the input parameter to generate mock observations). The best-fit exponential profile magnitudes of the mock observations are consistent with the input values within the errors, and the radii of the exponential profiles can be recovered with an accuracy of ≲ 50%. In the bottom two panels, we show the best-fit exponential profile magnitudes for the mock observations in the F200W and the F115W bands. Note that the morphology of the exponential profile is fixed when fitting these two bands. In most cases, the best-fit exponential profile magnitudes from mock observations are consistent with the input value within errors, except for J159-02 in the F115W band, which we report as a tentative detection. These results indicate that the host galaxy magnitudes and their errors reported in this work are reliable.  \nTo summarize, the host galaxy detections of J0148+0600, J159-02, and J1120+0641 are reliable and cannot be explained by PSF inaccuracies. We report host galaxy detections for J0148+0600 in all three bands. For J159-02 and J1120+0641, we report host galaxy detections in the F356W and the F200W bands, as well as tentative detections in the F115W band. Mock observations show that the reported host galaxy magnitudes and their errors are accurate, validating our image fitting method.', 'REFERENCES': '```\nAnderson, J., & King, I. R. 2000, PASP, 112, 1360, doi: 10.1086/316632 Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068 Astropy Collaboration, Price-Whelan, A. M., Sip"ocz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f Ba˜nados, E., Venemans, B. P., Decarli, R., et al. 2016, ApJS, 227, 11, doi: 10.3847/0067-0049/227/1/11 Ba˜nados, E., Schindler, J.-T., Venemans, B. P., et al. 2023, ApJS, 265, 29, doi: 10.3847/1538-4365/acb3c7 Becker, G. D., Bolton, J. S., Madau, P., et al. 2015, MNRAS, 447, 3402, doi: 10.1093/mnras/stu2646 Behroozi, P., Conroy, C., Wechsler, R. H., et al. 2020, MNRAS, 499, 5702, doi: 10.1093/mnras/staa3164 Bertin, E. 2011, in Astronomical Society of the Pacific Conference Series, Vol. 442, Astronomical Data Analysis Software and Systems XX, ed. I. N. Evans, A. Accomazzi, D. J. Mink, & A. H. Rots, 435 Bertin, E., Mellier, Y., Radovich, M., et al. 2002, in Astronomical Society of the Pacific Conference Series, Vol. 281, Astronomical Data Analysis Software and Systems XI, ed. D. A. Bohlender, D. Durand, & T. H. Handley, 228 Bordoloi, R., Simcoe, R. A., Matthee, J., et al. 2023, arXiv e-prints, arXiv:2307.01273, doi: 10.48550/arXiv.2307.01273 Bosman, S. E. I., \' Alvarez-M\'arquez, J., Colina, L., et al. 2023, arXiv e-prints, arXiv:2307.14414, doi: 10.48550/arXiv.2307.14414 Byler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, ApJ, 840, 44, doi: 10.3847/1538-4357/aa6c66\n```'}