EIGER. V. Characterizing the Host Galaxies of Luminous Quasars at z ≳ 6

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 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 host-to-quasar-flux ratios of ∼1%–5%. Spectral energy distribution fitting suggests that these quasar host galaxies have stellar masses of M * ≳ 1010 M ⊙. 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 BH-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.


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 Corresponding author: Minghao Yue myue@mit.edu* NHFP Hubble Fellow galactic 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).
Meanwhile, 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 et al. 2005;Hopkins et al. 2008), producing a correlated growth of the SMBH and its host galaxy (e.g., Croton 2006;Storchi-Bergmann & Schnorr-Müller 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ò 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.
In 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ñados 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).
Nevertheless, 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.
Although 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).
This 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.
In 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.
This 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 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), andBordoloi et al. (2023) for more information about the EIGER project.
We 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.

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.
There are some noticeable differences between the image reduction of this work and that of previous EIGER papers.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.

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.fitsfiles, 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).

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.

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  use photutils 1 to build effective PSFs, which uses the algorithm described in Anderson & King (2000).The detailed approach is as follows.
For 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, 1 https://photutils.readthedocs.io/en/stable/index.html 2 The specific FWHM limits are (0.′′ 0569, 0. ′′ 0695) for F115W, (0. ′′ 0729, 0. ′′ 0789) for F200W, and (0. ′′ 128, 0. ′′ 152) for the F356W filter.These values correspond to the 3σ limits as measured by Zhuang & Shen (2023) as 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.
It 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.
The 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.
Specifically, 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 (1) where σ 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).
The 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×FWHM are 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.
All PSF models, the PSF star lists, and the code to construct the effective PSF with error maps will be made available online upon publication.

Image Fitting
We use a point source component to describe the quasar and an exponential profile (i.e., a Sérsic 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érsic profiles when there are other bright galaxies close to 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érsic profiles, the half-light major and minor radii and the position angle of the Sérsic profiles, the Sérsic index of the Sérsic profiles (except for the quasar host galaxy which is forced to have n = 1), and the background level.
We fix the Sérsic 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érsic index n ≈ 1 (e.g., Yue et al. 2021).We will discuss more about this point at the end of this Section.
It 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., where ϵ 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.
For 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 Å-break for z ∼ 6 − 7 galaxies.As we will later show in this Section, the quasars with host galaxy detections in our sample have We fit the F356W images from individual visits separately, and determine the best-fit parameters and their

Yue et al.
errors 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 , where χ 2 ν (1p) is the reduced χ 2 for the single-PSF model, and χ 2 ν (1p1e) 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: 1.The error of the quasar host galaxy's magnitude is smaller than 0.3 (i.e., a detection with > 3σ significance); 2. The half-light radius satisfies R e,maj < 0. ′′ 9 and R e,min > 0. ′′ 1; 3. The axis ratio satisfies 0.3 < R e,min /R e,maj < 1; Here, 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.
For 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.
We notice 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, leaving 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.
We 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.
Our 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.
Figure 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%.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.The "(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 on how we validate the host galaxy detections and why we adopt the detection criteria in Appendix A.
The 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 reported 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.
We also try to fit the images without fixing the Sérsic index of the host galaxy.We find that the Sérsic 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érsic 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.
We provide more information about individual quasars in Section 6.1.

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.
We 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 Å, τ 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 normalization of the SED.All other parameters remain largely unconstrained.
We 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.
For 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 ⊙ .
Since 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   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.5for more details.Figure 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.
two-dimensional grism data, which is beyond the scope of this paper.Finally, 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 roughly 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.4). 2 The errors are estimated from prospector posterior distributions.For quasars with non-detections, the upper limits correspond to ∆mF356W = 3.5 (Section 3.3).3 The absolute magnitude at rest-frame 1500 Å from SED fitting.

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: In 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.

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 and 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.
Figure 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 Å (L 5100 = λL λ (5100 Å), 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. The

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.
Figure 10 shows the comparison between the Mg iibased 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 iibased BH masses of high-redshift quasars may be systematically larger by ∼ 0.1dex than their Hβ-based BH masses.

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 Figure 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).
Several previous studies have also used JWST NIR-Cam 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.
Recent 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.
Some 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- Flux [10 −18 erg s 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.
luminosity 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, local 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.
The stellar masses of the high-redshift quasar host galaxies still have large uncertainties.One of the main reasons 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.

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 et 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.
The 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): with 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).Since 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 high- (L 5100 ∼ 10 45 erg s −1 ) and thus have less severe selection biases compared to the luminous quasars in our sample.
Finally, 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 6.1.1.J0100+2802 J0100+2802 (z = 6.31) is a hyperluminous quasar discovered by Wu et al. (2015).Due to its high luminosity, J0100+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 ∼ 850M ⊙ /yr (Wang et al. 2019a), which exhibit clumpy structures (Fujimoto et al. 2020).Eilers 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).Zhuang & Ho (2023).The filled contours represent the sample of mock quasars at z ∼ 6 with L5100 > 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 MBH − M * relation.Nevertheless, the quasars in this work still have larger MBH/M * compared to the median of the mock quasar sample.
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.6.1.4.J159-02 J159-02 (z = 6.381) was initially reported by Bañados 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.Venemans 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.
We 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.Bosman 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.6.1.6.J1148+5251 Yue et al.
J1148+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.

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.
We 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.
In 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.
Another 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-band 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.
It 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.
We now consider the method for image fitting.Most studies of quasar host galaxies (including this work) uses a point source plus a Sérsic profile to describe the image of quasars.The specific setting for the Sérsic profiles varies between the studies.Ding et al. (2023) left all Sérsic parameters free, Zhuang & Shen (2023) fixed the center of the host galaxy to the quasar's position, and this work fix the Sérsic 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).
The 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érsic 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).
Using 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).This 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.

CONCLUSION
We present NIRCam observations of six quasars at z ≳ 6 observed by the EIGER project.We use NIR-Cam 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: 1. We detect the host galaxy of J0148+0600 in all the three bands, as well as the host galaxies of J159-02 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 ⊙ .
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.
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).
4. The quasars in the EIGER sample have M BH /M * ∼ 0.15, which is ∼ 2 dex larger than low-redshift 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.
The 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.
We 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 observations are associated with program ID #1243.2013,2018), psfMC (Mechtley 2014), webbpsf (Perrin et al. 2014), jwst APPENDIX Figure 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.
To 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.

Figure 1 .
Figure 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%.

Figure 2 .
Figure2.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 MCMC posterior 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.

F356WFigure 3 .
Figure3.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 × FWHMPSF.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érsic profiles when fitting the images.Note that (1) the compisite error is a combination of random noises and PSF errors (see Equation2); (2) the "Host Model" images include the two close companions of the quasar.

F356WFigure 4 .
Figure 4. Same as Figure3, 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.5for more details.

F356WFigure 5 .
Figure5.Same as Figure3, 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.4for more details.

F356WFigure 6 .
Figure6.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.

F356WFigure 7 .
Figure7.Same as Figure6, 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.
[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-ponents, 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.

Figure 10 .
Figure 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.

Figure 11 .
Figure 11.The MBH − 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 fromKormendy & Ho (2013) andReines & Volonteri (2015) for comparison.EIGER quasars have MBH/M * ∼ 0.15, which is ∼ 2dex larger than the low-redshift relations.Note that theKormendy & Ho (2013) relation uses bulge masses instead of total stellar masses and is located on the left side of theReines & Volonteri (2015) relation.

Figure 12 .
Figure12.Evaluating the selection bias of the luminous quasar sample.The black contours represent the distribution of low-redshift AGNs inZhuang & Ho (2023).The filled contours represent the sample of mock quasars at z ∼ 6 with L5100 > 10 46 erg s −1 , following the method described inLi 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 MBH − M * relation.Nevertheless, the quasars in this work still have larger MBH/M * compared to the median of the mock quasar sample.

Table 1 .
The quasar sample of the EIGER projectThe SMBH masses calculated using the Mg ii broad emission line.The errors in the table only contain random errors; the systematic errors of the black hole mass is about 0.4 dex.References for the quasar properties.

Table 2 .
The Results of Image Fitting

Table 3 .
The properties of the quasars and their host galaxies The errors listed in the Table are random errors from MCMC.The systematic errors of log MBH is ∼ 0.3dex is dominated by the scatter of the empirical relation (Equation