Simulated host galaxy analogs of high-z quasars observed with JWST

The hosts of two low-luminosity high-z quasars, J2255+0251 and J2236+0032, were recently detected using JWST's NIRCam instrument. These represent the first high-z quasar host galaxy stellar detections and open a new window into studying high-z quasars. We examine the implications of the measured properties of J2255+0251 and J2236+0032 within the context of the hydrodynamic simulation BlueTides at z = 6.5. We find that these observed quasars fall on the BlueTides stellar to black hole mass relation and have similar luminosities to the brightest simulated quasars. We predict their star formation rates, estimating approximately $10^{2-3}$ $M_{\odot}/ \rm yr$ for both quasar hosts. J2255+0251 and J2236+0032's host galaxy radii also fall within estimates of the radii of the simulated host galaxies of similar luminosity quasars. We generate mock JWST NIRCam images of analogs to the observed quasars within BlueTides and perform a point source removal to illustrate both a qualitative and quantitative comparison of the measured and simulated radii and magnitudes. The quasar subtraction works well for similar luminosity quasars, and the recovered host images are consistent with what was observed for J2255+0251 and J2236+0032, further supporting the success of those observations. We also use our mock imaging pipeline to make predictions for the detection of J2255+0251 and J2236+0032's hosts in upcoming JWST observations. We anticipate that the simulation analogs of future high-z quasar host discoveries will allow us to make accurate predictions of their properties beyond the capabilities of JWST.

In the nearby Universe, there are tight correlations between galactic properties such as stellar mass, black hole mass (i.e., the Magorrian relation, see Magorrian et al. 1998;Croton 2006;Schutte et al. 2019), bulge mass (Kormendy & Ho 2013;Reines & Volonteri 2015), velocity dispersion (Ferrarese & Merritt 2000;Munari et al. 2013;Sharma et al. 2021), size, and luminosity (Kormendy 1977;Danieli & van Dokkum 2019).The relationship between SMBHs and their host galaxy properties has been studied extensively at low redshifts due to their abundance and proximity.In the local Universe, direct measurements of black hole mass using reverberation mapping have allowed measurements of galactic scaling relations such as the Magorrian relation (Bentz & Katz 2015;Bentz & Manne-Nicholas 2018).
Prior to the launch of the James Webb Space Telescope (JWST), observations of high-z quasar host galaxies had only been attainable at sub-mm wavelengths with instruments such as the Atacama Large Millimeter/sub-millimeter Array (ALMA) (Schleicher et al. 2010;Venemans et al. 2012;Wang et al. 2013;Venemans et al. 2016;Decarli et al. 2018;Yue et al. 2021).These ALMA observations were successful in detecting [CII] emission, a strong tracer of the dusty interstellar medium.The rest-frame far-infrared emission observed at sub-mm wavelengths only traces the cold dust of the quasars' host.The rest-frame UV emission observed in the near infrared is a tracer for the stellar light of a quasars' host but requires careful subtraction of the quasar's light.At high-z, black hole mass measurements have been made by combining single epoch emission lines such as MgII and CIV and continuum emission to quantify the velocity and radius of the broad line region of the quasar (McLure & Jarvis 2002).Attempts were made to observe  ∼ 6 quasar hosts with the Hubble Space Telescope (HST), resulting in upper limits on the host galaxy's brightness and stellar mass (Mechtley et al. 2012;McGreer et al. 2014;Marshall et al. 2020).Despite accurate point spread function (PSF) modeling and careful quasar subtraction, the quasar hosts were not detected.Using the BlueTides simulation, Marshall et al. (2022b) predicted that approximately 50% of high-z quasar hosts would be detectable with an equivalent JWST filter and exposure time.
Measurements of detected high-z quasar hosts will allow the precise extension of galactic scaling relations to high-z through direct stellar mass measurements.This will extend our knowledge of the co-evolution of black holes and their host galaxies providing vital information into how the first galaxies formed.It will also shed light on the accelerated growth of SMBHs in early galaxies and their formation histories.
The deconvolving of the quasar from its host galaxy's stellar light requires careful fitting of the point source and galactic light profile.This ensures an accurate estimation of parameters such as the Sérsic magnitude and radius.Ding et al. (2023) used the galight software (Ding et al. 2020) to decouple the stellar light from the quasar by modeling the quasar as a two dimensional PSF and Sérsic profile of the host galaxy.The Sérsic profile is defined as where  1/2 is the half light radius we define as the Sérsic radius, n is the Sérsic index corresponding to increasing steepness in the log brightness and log radius space,   ∼ 2 − 1 3 , and   is the intensity at the half light radius.This is done separately in both the F356W and F150W filters.This was successful except for J2255+0251 which was only detected in F356W.
Soon after, three other high-z quasar hosts were detected in the JWST Emission-line galaxies and Intergalactic Gas in the Epoch of Reionization (EIGER) project.EIGER is a JWST Guaranteed Time Observation (GTO) program combining NIRCam and its accompanying deep wide-field slitless spectroscopy (WFSS) instrument.Part of EIGER's first observations have been published detailing observations of 6 quasars between redshifts 5.9 and 7.1 (Yue et al. 2023).The quasars were imaged in the NIRCam filters F115W, F200W, and F356W with wide-field slitless spectroscopy from F356W.Three out of six of these quasars' host galaxies were detected (Yue et al. 2023).The properties of the EIGER quasars with detected hosts, J0148+0600, J159-02, and J1120+0641, are found similarly except using psfmc (Mechtley 2014).From the PSF subtraction and fit, the UV magnitude, galactic radius, Sérsic index, ellipticity, and stellar mass were also measured (See tables 1 and 2 for an overview of the properties of these quasars).Stone et al. (2023b) also detected the host of J2239+0207 at  ≃ 6.25 by scaling the PSF to match the quasar flux and finding a discrepancy in the difference between the two.We do not include this detection's measurements in our paper to stay consistent with detections that are performed with Bayesian point source removal software, i.e., psfmc (Mechtley 2014) and galight software (Ding et al. 2020).With non-Bayesian methods, Stone et al. (2023a) set upper limits on the galaxy flux from five high-z quasars that were not detected with JWST.These quasar host detections and non-detections have provided a first glimpse into measuring the black hole to stellar mass relationship in high-z quasars.Although Ding et al. (2023)'s black hole mass measurements seem to be in line with the low-z relation, many other observations have measured overmassive black holes (Yue et al. 2023;Stone et al. 2023a).For example, the EIGER quasars are reported to have black hole to stellar mass ratios that are up to 2 dex higher than the local relation (Yue et al. 2023).
As more high-z quasar hosts are detected with JWST, questions about the scatter of key high-z black hole scaling relations have arisen.The errors from the sensitive point source subtraction, calibration uncertainties, and other instrumental effects are difficult to disentangle from the observed properties.This motivates the use of large scale high-z simulations and mock imaging techniques to quantify this bias and explore discrepancies.Other similar simulations do not have large enough volumes to accurately capture high-z quasar formation.For example, the large scale hydrodynamic simulation Illustris (Vogelsberger et al. 2014) is nearly 45 times smaller in volume compared to BlueTides.The subsequent simulation to Blue-Tides, ASTRID, is approximately 4 times smaller in size (Bird et al. 2022).Due to its large box size, BlueTides presents itself as the ideal simulation to study high-z quasars with luminosities high enough (at least at the faint end of the observed quasar luminosity function) to be observable with current telescopes in a statistically robust way.
In this paper, we examine the measured properties of J2255+0251, J2236+0032, and the EIGER quasars with detected hosts within the distribution of  ∼ 6.5 galaxies with black holes simulated in the BlueTides simulation.In section 2, we provide a brief description of the BlueTides simulation.In section 3, we describe the specifics of converting the observed JWST parameters to comparable quantities in BlueTides.We then quantify our comparison between these quasars and the predictions of the BlueTides simulation.In section 4, we describe our methods for generating mock images and spectra of the BlueTides analogs of J2255+0251 and J2236+0032 and also performing a point source removal.We compare the results from the mock spectra, images, and point source removal in section 5. We make predictions for quasar host detectability in upcoming observations of the quasar analogs with the followup filters and exposure times in section 6. Lastly in section 7, we summarize our conclusions and provide future outlooks on upcoming JWST quasar measurements.We assume a ΛCDM cosmology with parameters from WMAP9 (Hinshaw et al. 2013) to maintain consistency with the BlueTides simulation as follows,  0 = 100 (km/sec)/Mpc × h = 67.4(km/sec)/Mpc,Ω  = 0.2814, Ω  = 0.0464, Ω Λ = 0.315,  8 = 0.820, and   = 0.9710.All magnitudes are expressed in AB magnitude.
Table 1.Relevant parameters for J2255+0251 and J2236+0032 as measured or determined from Ding et al. (2023), or in their original discovery papers: Matsuoka et al. (2016) and Matsuoka et al. (2018).All black hole errors are MCMC random errors and do not account for the usual 0.4 dex (Vestergaard & Peterson 2006) systematic errors from the scatter in the empirical relation used in Ding et al. (2023) * Parameters with an * were calculated from observed parameters using methods described in section 3.1.1.66 ± 0.10 * Parameters with an * were calculated from observed parameters using methods described in section 3.1.

Simulation Overview
BlueTides is a hydrodynamic simulation run between z = 99 and z = 6.5 (Feng et al. 2015;Tenneti et al. 2018;Ni et al. 2020) using the MP-Gadget Smoothed Particle Hydrodynamics (SPH) code.The simulation contains 2 × 7040 3 particles in a box with sidelength = 400ℎ −1 Mpc and includes star formation, feedback processes, and black hole accretion in its computed physics.It contains 200 million star-forming galaxies (Huang et al. 2018).Owing to the finite volume of BlueTides, the most massive black hole is 7.7 × 10 8  ⊙ at z = 7 and 1.4 × 10 9  ⊙ at z = 6.5.Higher black hole masses are unlikely to be created within the 400ℎ −1 Mpc box.Previous studies have shown good agreement between BlueTides and observational constraints on the high-z Universe, including quasar number counts (Tenneti et al. 2017;Ni et al. 2020) and galaxy stellar mass relationships (Wilkins et al. 2017).All these previous studies were at z = 7.We present the first analysis of BlueTides at z = 6.5.
The BlueTides AGN formation and feedback model was developed in Di Matteo et al. (2005).The MassiveBlack I and II simulations (Matteo et al. 2012;Khandai et al. 2015) share the same AGN formation prescription as BlueTides.To align with predictions from the direct collapse scenarios for black holes (Yue et al. 2014), BlueTides uses a black hole seed mass of 10 5.8  ⊙ .They occur within haloes that are larger than approximately 10 10.8  ⊙ .To avoid studying black holes that are too close to the black hole seed mass, we limit our study to 108,000 black holes which have the highest black hole accretion rate at z = 6.5.These black holes are hosted within virial halos with a minimum mass of 10 11.0  ⊙ , maximum mass of 10 13.1  ⊙ , and average mass of 10 11.4  ⊙ .The maximum black hole mass for those not included in the first 108,000 galaxies is approximately 10 6.3  ⊙ .Thus, we limit our sample to include black holes with masses greater than 10 6.3  ⊙ when finding JWST analogs within BlueTides.This BlueTides subset includes black holes with bolometric luminosities (calculated using equation 5) less than 10 14.0  ⊙ hosted within galaxies with stellar masses less than 10 11.3  ⊙ .2020) is overplotted in blue for reference with the observations from the fit.The fit used is referred to as global fit A in Shen et al. (2020) where the lower luminosity portion of the fit is a flexible polynomial.Our smallest bins contain one quasar.

Quasar Luminosity Function (QLF)
In Figure 1, we show the z = 7 BlueTides QLF plotted alongside the z = 6.5 QLF.We compare this to observed quasars, including the parameterized QLF from Shen et al. (2020), which is overplotted for reference (global fit A in Shen et al. 2020).It is described by the following single power law: where  = 1.509 ± 0.058, log  * = −5.452,and log  * = 11.978± 0.055.Since the  ∼ 6 quasar luminosity function does not include many low luminosity quasars (with UV luminosities greater than about -21 mag), the typical double power law used to describe the quasar luminosity function is replaced with a single power law to avoid degeneracies in the parameter fit.Very few quasars have been found above z = 7, and at z = 6.5, BlueTides provides a sample of quasars that is more representative of those that have been discovered so far.We see that BlueTides produces a realistic quasar population, with quasars matching the lowest known luminosity quasars and spanning the lower end of the range of SHELLQs quasars.Due to the finite volume of BlueTides, the BlueTides QLF does not extend far into the brightest SDSS quasars or include many quasars matching the luminosity of the EIGER quasars.We also find that the lower redshift QLF peaks an order of magnitude higher in quasar bolometric luminosity as expected for older black holes.

Black hole bolometric and Eddington luminosity
To identify mock quasars within BlueTides with the properties of observed high-z quasars, we use the mass accretion rate of the black hole to calculate its corresponding bolometric luminosity.We can then calculate their corresponding UV magnitudes, which is accomplished using a fit from measured quasar spectral energy distributions (SEDs).
The UV magnitudes of J2255+0251 and J2236+0032 were previously measured in the SHELLQs survey.The absolute UV quasar magnitudes are −23.87±0.04 and −23.66±0.10 for J2255+0251 and J2236+0032, respectively.For the EIGER quasars, the UV magnitudes for the quasars are -27.08,-26.47, and -26.44 for J0148+0600, J159-02, and J1120+0641, respectively.As described in Qin et al. (2017), we can convert a quasar bolometric luminosity to a magnitude in UV at 1450 Å using Hopkins et al. (2007) extrapolated from the B-band 4344 Å magnitude using Blanton & Roweis (2007), where There is no analytic form for the inverse of equation 4. We solve numerically for the bolometric luminosity, which was calculated for all quasars explored in this work and is tabulated in 1.
From the black hole accretion rate, we approximate  bol of the black hole with where  BH is the quasar mass accretion rate and  is the radiative efficiency.We take  = 0.1 assuming an equipartition magnetic field (Gruzinov 1998).We compare the bolometric luminosity to the Eddington luminosity of a black hole, which is where  is the opacity of the black hole.If we assume the highly accreting material around a black hole to be mostly ionized hydrogen so that Thomson scattering describes the opacity reasonably well, we can set  =  T / p , where  T is the Thomson cross section and  p is the mass of a proton.We note that Bondi (or spherical) accretion is assumed in BlueTides.

Black hole mass
Figure 1 shows that with luminosities near 3 × 10 12  ⊙ , J2255+0251 and J2236+0032 fall within the most luminous quasars in BlueTides as described in section 2.1.This allows us to compare other properties of the observed quasars with analogs in BlueTides.The EIGER quasars are in the BlueTides luminosity range where there are only a handful of quasars.With so few quasars, the binned luminosity function does not extend to such a high luminosity as seen in Figure 1.
Figure 2a shows the simulated black hole masses for quasars in BlueTides with the measured quasar masses overplotted.Both Ding et al. (2023) quasars' measured black hole mass and luminosity fall within the expected BlueTides distribution.The quasars have tight errorbars and fall within the expected black hole mass range.J2236+0032 has a black hole mass above that predicted by BlueTides for quasars with similar bolometric luminosities.The EIGER quasars J159-02 and J1120+0641 also fall within the expected black hole mass range, but the most massive quasar, J0148+0600, has a mass of 10 9.892  ⊙ which is above the mass limit of BlueTides.

Stellar mass and SFR
Figures 2b through 3a show the joint distribution of black hole masses and luminosity with stellar mass.The stellar mass was measured as 10 10.53  ⊙ and 10 11.12  ⊙ for J2255+0251 and J2236+0032, respectively (Ding et al. 2023).The black hole masses and stellar masses of J2255+0251 and J2236+0032 follow the black hole mass to stellar mass ratio.However, the EIGER quasars lie 1-2 dex above the Kormendy & Ho (2013) relation and do no lie within the BlueTides distribution.The EIGER quasars have bolometric luminosities that are brighter than the BlueTides quasars which may bias the comparison.The uncertainty and difficulty of removing the quasar from its host (as described in section 4.3) could lead to biased stellar mass measurements.Understanding this observational bias will be the goal of future work.
In Figure 3b, we show where the quasars would lie on the star formation rate versus bolometric luminosity distribution.This allows us to estimate their SFR using the BlueTides distribution of quasars.We anticipate a star formation rate between approximate 10 2−3  ⊙ /yr for both quasars.For the more luminous EIGER quasars, we can place an approximate lower limit of 10 2.75  ⊙ /yr.Followup measurements in the sub-mm with instruments such as ALMA of the H emission line with JWST will directly measure the SFR.This shows the potential of using high-z simulations to both constrain and support future quasar observations.

Galaxy radii
The half light radius (  1 and Table 2 and range between approximately 0.5 to 3 kpc. In Figure 4, we compare the half light radii estimates of BlueTides galaxies from Marshall et al. (2022a) against their bolometric luminosity.These radii were calculated by generating rest-frame images in the standard FUV top hat filter (1500 Å).They do not account for a PSF or any noise.The orientation and size of the field of view were found to make a difference less than 0.04 dex and are thus neglected in our comparisons.To calculate these half light radii, Marshall et al. (2022a) performed aperture photometry with 20 circular apertures logarithmically spaced between 0.05 to 2 kpc and interpolated between them to find the half light radius.We find that the half light radii measured in BlueTides and those fit in Ding et al. (2023) are consistent as seen in Figure 4.For the EIGER quasars, we find that J1120+0641 is consistent with BlueTides.J159-02 and J0148+0600 have radii measurements of 2.64 ± 0.1kpc and 2.23 ± 0.11kpc, respectively, above our sample in radius.
In recent results with the JWST ERO SMACS 0723 Field for galaxies between redshifts 7 and 12 as seen in Adams et al. (2022), the maximum half light radius measured is near 0.8 kpc.This result is in line with similar HST measurements found by Holwerda et al. (2015) which measures sizes in the CANDELS survey.Holwerda et al. (2015) measured an average radius size of 0.5 ± 0.1 kpc.Similarly in Yang et al. (2022) which explore massive galaxies from the JWST GLASS survey with redshifts above 7, the half light radii found in the NIRCam filter F200W does not exceed 1.11 ± 0.11kpc (Galaxy ID 4397).These radii values fall near the average of those found for BlueTides galaxies that host high bolometric luminosity quasars and recently detected quasar hosts from Ding et al. (2023).
In summary, the quasars hosts detected in Ding et al. ( 2023) are described well by the results of the BlueTides simulation and are in line with other high-z galaxy radii measurements.

BlueTides analog selection
Finding the analogs of detected quasar hosts in BlueTides allows us to make predictions beyond what is possible with existing observations, and to assess the success of the simulation.We select the BlueTides analogs within the 4D property plane of black hole bolometric luminosity, black hole mass, stellar mass, and host galaxy radius.We find one analog per quasar within our joint distribution constraints.Each of these analogs is seen within the distributions shown in figures 2a through 4 in the same color as the respective observed JWST quasar with a black circle around the BlueTides galaxy.
We adjust the range to explore for each property to ensure that we find at least one analog per observed quasar.Our analogs are within 0.4 dex of the measured bolometric luminosities for both quasars.They are within 0.4 and 0.01 dex of the measured black hole masses and 0.2 and 0.01 dex of the measured stellar masses for J2236+0032 and J2255+0251, respectively.The J2255+0251 Blue-Tides analog is within the error bars of the measured galaxy radius, but the J2236+0032 analog is about 0.7 kpc larger in galaxy radius than the observed host.In Figure 2a, the J2255+0251 analog is seen with a black hole mass and bolometric luminosity nearly matching the observed quasar host.For J2255+0251, we choose the BlueTides galaxy with the closest properties as our analog, but for J2236+0032, we select the brightest quasar from approximately 5 analogs to remain closer to the magnitude of the detected host in Ding et al. (2023).This allows for comparison of the point source subtraction between two analogs with varying magnitudes for J2236+0032.With a BlueTides galaxy analog for both J2236+0032 and J2255+0251, we are able to make comparisons and predictions between BlueTides and these high-z quasars.

Dust Calibration of BlueTides Star Particles
A key component in comparing BlueTides to observations is ensuring accurate dust modeling during the mock imaging process described in section 4.3.Each star has a simple exponential dust model for the stellar natal cloud if stars are younger than 10 Myr.The optical depth of the natal cloud is calculated as follows where  is the metallicity of the birth cloud,  = −1, and  is the wavelength at the particular optical depth.The interstellar medium (ISM) dust contribution is calculated for each star particle regardless of age as follows where  = 10 4.6 ,  = −1, and  metal is metal density in the line of sight as a function of  which is integrated out to , the distance to the star. and  are calibrated from the observed UV galaxy luminosity function at z = 7 as shown in Marshall et al. (2022a).

Mock images and spectroscopy
To compare model quasars directly with existing observations, we created mock JWST NIRCam images and processed these to recover host properties in analogy with observations.We generate mock images using the Python package, SynthObs (Wilkins 2019) including both a mock background and simulated dust.The effects of adding both background and dust to a raw simulated galaxy image increase the galactic magnitude by ∼ 0.4.Within SynthObs, we assign a flux to each star particle based on each particle's metallicity, mass, and age using Binary Population and Spectral Synthesis (BPASS) (Eldridge et al. 2017;Stanway & Eldridge 2018) processed through Cloudy (Ferland et al. 2017).We select the BPASS model which uses the standard Salpeter stellar initial mass function (IMF) (Salpeter 1955) up to 300 ⊙ .We assume 90% of Lyman-continuum photons escape.The dust model for each stellar particle is described in section 4.2.We then use the star particle flux and redshift to generate an SED for each star.We also add a point source to emulate a quasar from Ding et al. (2023) with the same observed apparent magnitude for each filter as shown in Table 1.Although we could model the SED of the BlueTides black hole in the analog galaxy, we choose to match the quasar magnitudes in the mock images to ensure a more similar point source removal process.The point source removal is highly dependent on the ratio of galaxy to quasar flux where a relatively brighter host typically allows for an easier host detection.
The filter transmission curves are convolved with each SED to get the appropriate filter specific flux for each particle.These SEDs get added together to form the overall spectrum of the galaxy.Finally, we convolve the spatial grid of filter corrected fluxes with mock JWST NIRCam PSFs from WebbPSF (Greenbaum et al. 2015).For the higher resolution F150W filter, we do not incorporate any supersampling and generate a 25 x 25 kpc field with 0.031 arcsec/pixel with 144 pixels.The F356W filter's images are also generated with a 25 x 25 kpc field with 0.063 arcsec/pixel supersampled to 0.031 arcsec/pixel with 140 pixels.For each image at z = 6.5, we note that 1 pkpc = 0.18 arcsec.
We also add background noise to the image to allow a more realistic fitting process.The background is derived from the aperture flux limit for each JWST filter at a 10 depth and the size of the aperture (Rieke et al. 2023).The aperture flux limit for a particular filter is the minimum observable flux for a 10 observation of a point source for a specific exposure time.More explicitly, the background at each pixel for a given filter is background(pixel) = √︃ (F limit /10) 2 /A aperture × N(0, 1), ( where  limit is the aperture flux limit,  aperture is the area of the aperture, and N is a random value drawn from a normal distribution at variance 1.We use an exposure time of 3.1ks to maintain consistency with Ding et al. (2023) for our mock point source removal, unless specified otherwise.

Mock Spectra of BlueTides Analogs
We make redshifted mock spectra for the JWST analogs of J2255+0251 and J2236+0032 and 15 other close analogs with similar bolometric luminosities (some of their mock images can be seen in Figure 9).We infer the SED in the measured two-band photometry as described in section 4.3.The magnitudes and flux ratios for both the observed and simulated quasars are summarized in Table 3. Notably, the BlueTides host analogs have magnitudes broadly similar to the magnitudes of the observed quasar hosts.The upper limit for the host of J2255+0251 in F150W falls significantly below the magnitude of the BlueTides analog.For the filters where hosts were detected, this provides evidence that the primary analogs chosen are representative of the observed quasar hosts.For both observed quasar hosts, the spectra fit must be lower in F150W and higher in F356W to match the observed magnitudes.This may indicate that the observed host brightness was underestimated, due to the difficulty of the quasar subtraction.
In Figure 5, we plot the galaxy spectra for both the primary analogs and additional analogs to compare to Extended Data Figure 5 in Ding et al. (2023).The mock spectra also look similar to high-z galaxy spectra observed with JWST in Maiolino et al. (2023) and Boyett et al. (2024).This indicates reasonably accurate spectral modeling with SynthObs of the analogs.
Each analog has a similar spectral shape.The increasing amplitude of the spectra is roughly correlated with stellar mass but also depends on differences in average age, metallicity, and metal surface density for stars in the galaxies.All analog mock spectra also show the diversity of emission lines potentially present in the high-z quasar hosts which are much more difficult to detect in high-z quasar systems compared to the quasar spectra.
Although we do not perform a quantitative comparison, the spectra of the galaxy analogs look significantly different than the median template from SED inference in Ding et al. (2023).Blueward of the Balmer break, the analogs show a much steeper spectra, instead being more consistent with their template corresponding to the 1 upper limits.This results in the significantly higher flux in both primary BlueTides analogs in the F150W filter than the observed hosts.We predict this difference in spectral shape between all analogs and the observed hosts to be primarily due to a difference in stellar age in modeled spectra.The median inferred spectra in Ding et al. (2023) uses a stellar population with an age of 370 Myr and 500 Myr for the galaxies of J2236+0032 and J2255+0255, respectively.For the primary BlueTides analogs, the average age is approximately 100 Myr resulting in a much steeper UV slope.However, the BlueTides galaxy spectra appear broadly consistent with the inferred spectra from the lower age range modeled in Ding et al. (2023).

Mock Rest Frame Images of BlueTides Analogs
With the intrinsic properties simulated in BlueTides, we also create rest frame images of galaxies.These provide a glimpse into the evolution of detected hosts of the Ding et al. (2023) and other SHELLQs quasars that is unachievable with observations.We use SynthObs to generate rest frame images of the BlueTides analogs.These images do not have the observational effects of the other mock images, such as noise and PSF effects, and have a much higher resolution.We show the stellar mass distribution and recent star formation in Figure 6.We also show images in idealized top hat filters in the far UV, visible (V), and infrared (H) as well as the stacked RBG version.These images provide a higher resolution snapshot of the star particles in the BlueTides analog than the blurred redshifted images shown in Figures 7a through 8b.The galaxies are notably compact in the left three panels in stellar mass, SFR, and intrinsic FUV.The intrinsic FUV panel traces the star formation and stellar mass.This is in line with results from Marshall et al. (2022a) which predicted that galaxies hosting quasars are more compact intrinsically relative to the entire BlueTides sample.When we include dust attenuation in the four right panels of Figure 6, we see a flatter emission profile as predicted in Marshall et al. (2022a).This results in the appearance of an extended galaxy.The ring around the center of these galaxies shows a region where the dust is particularly strong.

Point source removal
After the creation of mock images described in section 4.3, we compare the mock imaged BlueTides galaxies before and after the point source removal process to further increase the confidence that the Ding et al. (2023) hosts appear to be plausible detections.To do this, we perform quasar subtraction using the software psfmc on our model images, using a mock star image as our input PSF model.For our mock star image, we use the same model PSF from WebbPSF used to create our quasar image but with a different realization of the Figure 5. Mock spectra from SynthObs of the two primary analog galaxies and 15 other galaxies hosting black holes with similar bolometric luminosities as those in Ding et al. (2023).We also show the magnitudes of the primary analog and observed host in both F356W and F150W for both quasars for comparison.All mock spectra have similar shape but differ in amplitude due to stellar mass, stellar age, metallicity, and metal surface density.background and shot noise (see Marshall et al. (2021)).This is an idealized case, where the PSF shape is known perfectly up to noise variation.While the most significant challenge in high-z quasar host detection is often the PSF mismatch error, our approach shows a more ideal case, which is useful for understanding the optimal outcome from these subtractions.Even with a perfect PSF shape, we will inevitably end up with some PSF subtraction error, which we begin to understand qualitatively here, and our aim is to see how this might appear in the Ding et al. (2023) host images.
We recover the properties of the host galaxy and quasar from the imaged galaxies using the software psfmc to do point source extraction.psfmc was written specifically for quasar subtraction unlike other algorithms that perform similar 2D surface brightness fitting such as galfit in Peng et al. (2002).psfmc uses the Affine Invariant Markov chain Monte Carlo (MCMC) ensemble sampler implemented in emcee (Foreman-Mackey et al. 2013).The posterior sampled includes priors specified by the user and the following Gaussian joint likelihood for all pixels Simulated host galaxy analogs of high-z quasars observed with JWST 9 (10) where (pixel) includes the variance at that pixel in both the PSF and raw image and I CM (pixel) is the intensity of the convolved model at that particular pixel given some set of model parameters .For each run of psfmc, we pass in a PSF generated by an imaged point source with a total flux of 6 × 10 5 nJy and convolved with the observed filter's mock PSF from WebbPSF using SynthObs.
The output of psfmc contains the posteriors of all parameters and the final galaxy image with the quasar subtracted.The intermediate images derived are the raw model and the convolved model.The raw model contains the highest probability galaxy and quasar surface brightness distribution.We convolve the raw model with the PSF to form the convolved model which is subtracted from the input observed image to determine the residuals.The third column of the panels labeled "model" in figures 7a through 7d show the convolved model.
We employ the sampling techniques from psfmc in the same filters as those observed in JWST.We model each quasar and its host as a point source and Sérsic profile defined galaxy system.We specify uniform priors for their positions, the quasar magnitude, half light semi-minor and major axis (the latter is taken as the half light radius), magnitude, and proper angle.We set the Sérsic index to 1 which is commonly done (i.e., in Ding et al. 2023 andYue et al. 2023) and has been shown to be consistent with simulated high-z galaxies (Tacchella et al. 2015).We sample the following free parameters with uniform priors:  filter,quasar ,  filter,Sérsic ,  1/2 ,  1/2,b , and Sérsic angle.The uniform prior ranges are as follows for all Sérsic profile fitting and point source removal in this work: We initiate 64 walkers and discard the initial 200-400 samples to allow proper convergence of our sampler.The values of our results are drawn from the last 1k sample ensuring accurate parameter error estimation.The mock images are shown and compared qualitatively in the subsequent sections 5.4 and 5.5.

Mock Images of BlueTides Analogs
We also generate realistic JWST NIRCam mock observations for each analog.This serves as a useful comparison as we are able to compare the exact BlueTides galaxy with and without the source extracted.This is evidently impossible with actual observations but allows us to test how accurately the point source subtraction reveals the underlying stellar light.
In Figures 7a through 8b, we show six panels of our mock observed images, depicting I) the BlueTides galaxy without a quasar imaged, II) the same BlueTides galaxy with a quasar imaged, III) the convolved model from psfmc described in section 4.3, IV) the host galaxy image with the quasar subtracted via psfmc, V) a normalized residual between the data and point source subtracted image from psfmc (i.e., (data -model) /  data−model ), and VI) the detected host JWST image from Ding et al. (2023).With psfmc, we are successful in extracting the 2D galactic brightness distribution and quasar parameters for all of our imaged BlueTides galaxies.This is promising as the analogs are particularly faint in the F150W filter.The actual host (galaxy only image) versus quasar subtracted (data -quasar image) all recover the underlying host with moderate visible error.In Figures 7b and 7c, the center of the host is missing flux.Each quasar subtracted host also features subtle to significant morphological changes from the actual host image.Through the normalized residuals in the fourth panel, we can quantitatively compare the data + quasar versus galaxy only images which ideally would be the same.There is a slight underextraction and overextraction around the core area containing the point source.This is an artifact of the point source removal and should be considered with point source subtraction on real images.
In Figures 7a through 7c, the J2236+0032 and J2255+0251 host analogs are similar to the detected hosts in Ding et al. (2023) in all filters.In Figure 7d, the J2255+0251 host analog detection in F150W is shown alongside the non-detection in Ding et al. (2023).Notably, this is also the dimmest host analog detection in our sample.
In Figure 8a, we show the mock images from a dim analog matching properties of J2236+0032.This analog is only slightly different in magnitude and black hole bolometric luminosity compared to the analog chosen in Figures 7a through Figure 7b.The dim analog hosts a black hole within 20% versus 50% of J2236+0032's bolometric luminosity and is approximately 1 mag fainter than the standard analog we chose.This analog is detected but with a less visibly accurate galaxy reconstruction compared to the original analog shown in Fig- ure 7a.This is seen qualitatively by comparing the galaxy and dataquasar panels where the dimmer galaxy has a larger corrupted quasar core.
In Figure 8b, we show the brightest galaxy in our BlueTides sample that hosts one of the 100 most luminous black holes to show the point source removal on a different galaxy morphology.This galaxy looks more elongated than any of the other analogs and may be undergoing a merger.The elongated morphology of this galaxy hosting one of the most luminous black holes in BlueTides supports further host galaxy comparison to those detected in Ding et al. (2023) and Yue et al. (2023).The size and luminosity of the galaxy seems to play a significant role in point source subtraction and will be explored in a future work.
In Figure 9, we show a wide array of mock images of BlueTides galaxies with similar bolometric luminosities to J2236+0032 and J2255+0251 in the F356W filter for comparison.At z = 6.5, these galaxies look substantially more diffuse than those detected in Ding et al. (2023) with some showing evidence of mergers.The EIGER galaxy hosts which contain quasars with higher luminosity resemble these larger black hole hosts.The J0148+0600 and J159-02 detected hosts show evidence of companion galaxies, and J1120+0641 is notably more diffuse.This supports further statistical studies on physical properties of host galaxies of high-z quasars to constrain evolutionary models.

Measured Properties in Mock Images
Qualitatively, the imaged galaxies look similar to the detected galaxies in Ding et al. (2023), with similar shapes and compactness.However for most of these mock images, our point source subtraction algorithm described in section 4.3 creates moderate to large errors in predicted radii and magnitude when extracting a point source.We compare the Sérsic radii and magnitudes quantitatively for the imaged galaxy with and without the point source.We extract both Sérsic radii and magnitudes as defined in equation 1 using the MCMC framework in psfmc to maintain consistency.We tabulate the mean and one sigma value from our sampled posteriors for each Sérsic radius in Table 4 for both the original galaxy before adding the point source and the radius estimated after point source subtraction.We   2023) is due to supersampling in both filters such that the pixel size in each filter varies by a factor of two (0.0315 arcsec/pixel for F356W and 0.0153 arcsec/pixel for F150W).Our pixels size is 0.031 arcsec/pixel in both filters.The corresponding pixel numbers between our images and the Ding et al. (2023) are not equivalent.
find that the results do not always agree to one sigma for both the original and point source subtracted galaxies indicating that point source subtraction is introducing bias.
The radii have a large discrepancy between simulated and observed values.We note that the largest bias in radius is found for the J2236+0032 analog in F150W shown in Figure 7b with the radii differing by about 0.4 pkpc and a percent difference of ∼ 30%.The radii are always underestimated from the actual value in our analog sample shown in Table 4.However, the radius of the brightest galaxy is recovered relatively accurately.This indicates that there is a larger error for smaller hosts that are more hidden behind the corrupted quasar core.This exemplifies the caution needed when using measured radii values from high-z hosts.We also compare these radii to the dust attenuated half light radii from Marshall et al. (2022a) as described in section 3.4.These radii are on the same order as the radii found through the other methods but do not match even those radii measured without an added quasar.However, this is expected as the radii from Marshall et al. (2022a) are not from a Sérsic profile fit.
We also measure the galaxy magnitude for the BlueTides analog with and without the point source removal process.We find that the magnitudes are predicted more accurately with each measured magnitude falling within approximately one sigma of the pre-point source removal value.Similarly to the radii measurements, we find that the post point source removal magnitude measurement is underestimated approximately 50% of the time compared to the actual Sérsic magnitude of the host.
Lastly, we compare the magnitude and radius pre-and post-point source removal for the brightest BlueTides galaxy within the first 100 most massive black holes and imaged in Figure 8b.As seen in Table 4, the radius is predicted the most accurately for this galaxy.Similar to the other quasar analogs, the observed magnitude shows a slight underestimation post point source removal.
We note that the errors we find for our analogs are smaller than those seen for the observed quasars and their hosts in Ding et al. (2023).Our errors include errors from our more idealized point source removal, as well as the noise in the image.In the observations in Ding et al. (2023), the uncertainty also includes larger systematic errors due to PSF mismatching and other effects caused by the real telescope imaging process.
Although we have not conducted a robust statistical study in this work, we see a trend of more accurately predicted galactic properties for brighter galaxies.With our mock imaging pipeline and point source removal algorithm, we can use BlueTides to quantify bias and make estimates for the over or underestimation of quasar host properties more generally.This will be the subject of a future paper using the large sample of BlueTides galaxies.

DETECTION PREDICTIONS FOR J2255+0251 AND J2236+0032
Follow up observations with JWST of J2255+0251 and J2236+0032 (detected in both F150W and F356W by JWST) will allow further constraints on the properties of their hosts.In JWST Cycle 2 GO # 3859 (PI: Onoue), these quasars will be followed up in five NIRCam filters.In Table 5, the filters and exposure times for this proposal are shown.In Figure 10, we show the mock images of galaxies without an added quasar for these filters and exposure times.We also extend the exposure time of the original filters (F356W and F150W) to 10ks in Figure 11 for comparison to the 3.1ks exposure times shown in figures 7a through 7d.We supersample in all filters with wavelengths larger than 2m (filters above F200W) to maintain consistency with our mock imaging in section 4. The analogs are visible in most filters but are notably fainter in the medium band filters (F250M and F480M) and at filters including wavelengths below 2 microns.We expect that a host detection for the more compact and faint host of J2255+0251 in GO #3859 will be very difficult in F115W, F200W, F250M, and F480M.The brighter and slightly elongated host J2236+0032 may be relatively easier to detect in all filters.In Figure 11 for F150W, the analog of J2255+0251 is still extremely faint which may indicate host detection in this filter is impossible.4. Mean and one sigma error of sampled values for each Sérsic radius and magnitude for BlueTides galaxies.Values are shown with and without the added quasars for each mock imaged analog.We show both quasar analogs, the dimmer J2236 analog, and the brightest galaxy hosting one of the first 100 most luminous black holes.Each Sérsic radius and magnitude was generated using uniform priors and psfmc.All errors shown on the properties are one sigma.The BlueTides galaxy indicates which of the two galaxies was used in each image.The right column shows the dust half light radius described in section 3.4 from Marshall et al. (2022a).Although a simulated version of the true hosts of J2255+0251 and J2236+0032, our detection predictions will be useful in future JWST observations in filter and exposure time selection of high-z quasars and as a test of the predictive power of BlueTides.

CONCLUSIONS
The launch of JWST has opened a new era in high-z quasar observation with J2255+0251 and J2236+0032 recently becoming the first high-z quasars with host galaxy detections (Ding et al. 2023).Soon after, three more quasar hosts were detected (Yue et al. 2023).To understand the properties of high-z quasars, we compare the recently detected quasar hosts to systems hosting black holes with similar bolometric luminosities in the BlueTides simulation.To summarize, we find the following results: (i) We find the measured host stellar mass, radii, and black hole mass of J2255+0251 and J2236+0032 are consistent with similar BlueTides quasars.The expected star formation rates of J2255+0251 and J2236+0032 from comparisons with quasars of similar luminosities is between log 10 SFR ∼ 2 − 3M ⊙ /yr.We also set a lower limit of 10 2.75  ⊙ /yr for the EIGER quasars.
(ii) We find that the newly detected EIGER quasar hosts Yue et al. (2023) do not match the BlueTides simulation host predicted properties and thus do not align with the low-z Magorrian relation.However, the bolometric luminosities of these quasars are at the bright end of the BlueTides sample (See the quasar luminosity functions in Figure 1).It is not yet clear whether the disagreement is due to an observational bias or a bias owing to the finite simulation volume where we do not have a statistical comparative sample.
(iii) With the criteria described in section 4.1, we select two analogs similar to J2255+0251 or J2236+0032 in stellar mass, black hole mass, and galactic radius.We compare the BlueTides analogs to the actual JWST detected quasar hosts and find the photometric galaxy magnitudes match closely between simulation and observation.We use mock imaging and a point source subtraction algorithm to make qualitative comparisons and begin to understand the effects of observational biases with these two analogs.We find that the mock images of BlueTides analogs are qualitatively similar to the detected hosts in Ding et al. (2023).However, the host detection remains difficult.This supports the lack of quasar host detections in Stone et al. (2023a).We do detect all hosts, but note that the success is highly dependent on the luminosity of the host galaxy.When we determine the half light radius and magnitudes for mock images of galaxies, we find that the recovered magnitudes agree mostly to within one sigma.The radii have a larger discrepancy between actual and post point source removal values with a maximum percent difference of  ∼ 30% or about 0.4 pkpc.The point source removal process consistently underestimates the radii in the galaxies tested in Table 4.This motivates future work to quantify these biases.
(iv) We make predictions for the detections of J2236+0032 and J2255+0251 in upcoming JWST observations.We predict that longer exposure times and longer wavelength filters will provide the highest probability for successful host detections but the J2255+0251 host may remain elusive even with a longer exposure time in the F150W filter.
Our matched comparisons to BlueTides show the utility of highz simulations to study high-z quasar hosts.Despite being tuned to low-z observations, BlueTides is still able to accurately predict high-z galaxy and black hole properties.In addition to confirming agreement of observed and simulated properties, BlueTides will serve as a tool to infer the properties of high-z quasar hosts from observations.

Figure 1 .
Figure 1.The BlueTides bolometric quasar luminosity function (QLF) is plotted for all black holes at  = 6.5 in black with Poisson errors in shaded grey.Similarly in red with errors in pink, the  = 7 BlueTides QLF is plotted for comparison.The Ding et al. (2023) quasars have bolometric luminosities similar to the brightest quasars in BlueTides.We also overplot the range of bolometric luminosity depth for two quasar surveys: the SHELLQs and SDSS quasars with  = 5.7 − 6.9 from Matsuoka et al. (2018) and Jiang et al. (2016) in tan and light blue, respectively.The z = 6 global quasar luminosity function from Shen et al. (2020) is overplotted in blue for reference with the observations from the fit.The fit used is referred to as global fit A inShen et al. (2020) where the lower luminosity portion of the fit is a flexible polynomial.Our smallest bins contain one quasar.

Figure 2 .Figure 3 .
Figure 2.Each plot displays a combination of stellar mass, black hole mass, and bolometric luminosity for the BlueTides black holes and the J2255+0251, J2236+0032, and EIGER (Yue et al. 2023) quasar properties overplotted.The J2255+0251 and J2236+0032 BlueTides analogs (see section 4.1) are plotted in green or purple for J2255+0251 and J2236+0032, respectively with a black circle.Left: Bolometric luminosity and black hole mass for simulated and observed quasars.The observed quasar mass is the measured value from NIRSpec observations.The black line indicates the luminosity limit for BlueTides quasars which is approximately 3 times the Eddington luminosity.Right: The same as above but for stellar mass versus bolometric luminosity.

Figure 4 .
Figure 4. Same as Figure 2 but for the half light radius versus bolometric luminosity for BlueTides galaxies.The measured properties of recently detected high-z hosts overplotted.

Figure 6 .
Figure 6.Rest frame images of BlueTide analogs.From left to right in each panel, we show the following i) stellar mass distribution in each pixel, ii) star formation rate of the z = 6.5 galaxy, iii) intrinsic far UV image without dust, iv-vi) far UV (FUV) images, visible (V), and near infrared (H) images with dust, and v) the stacked FUV, V, and H images.
(a) The BlueTides analog to J2236+0032 imaged in the F356W filter.(b) The BlueTides analog to J2236+0032 imaged in the F150W filter.(c) The BlueTides analog to J2255+0251 imaged in the F356W filter.(d) The BlueTides analog to J2255+0251 imaged in the F150W filter.

Figure 7 .
Figure 7. Left to right of each panel: The host galaxy without the quasar imaged (very faint in F150W), the mock imaged quasar, the PSF convolved model which includes both the point source and galaxy, the data with the quasar extracted, the residual (data -model) normalized by the standard deviation of the residuals, and a copy of a portion of Figure 2 from Ding et al. (2023) showing the corresponding detected quasar host.Color bar units are MJy/str.Each image is a cut out from the full image and is equivalent in size to the rightmost plot in each panel with a size of 1.5" and 3" for F150W and F356W, respectively.The smaller size of F150W in Ding et al. (2023) is due to supersampling in both filters such that the pixel size in each filter varies by a factor of two (0.0315 arcsec/pixel for F356W and 0.0153 arcsec/pixel for F150W).Our pixels size is 0.031 arcsec/pixel in both filters.The corresponding pixel numbers between our images and theDing et al. (2023) are not equivalent.

Figure 8 .
Figure8.Same panels as Figure7.Top: The closest analog to J2236+0032 with a moderately successful point source subtraction in the F356W filter.This analog is within 20% of the bolometric luminosity of J2236+0032 rather than within 50% like the analog presented in this work.Bottom: Brightest galaxy among the first 100 BlueTides galaxies with the most highly accreting black holes.All color bar units are MJy/str.Note that each mock image is 3" in width and roughly the same size as theDing et al. (2023) .

Figure 9 .
Figure 9. Mock images in the F356W filter of 15 galaxies in BlueTides with similar bolometric luminosities to J2236+0032 or J2255+0251.Half of the images show compact morphology similar to the observed high-z quasar hosts detected in Ding et al. (2023).These images are smoothed prior to imaging.Color bar units are MJy/str.

Figure 10 .
Figure 10.Mock imaged BlueTides analogs in the filters and exposure times described in Table5.These images allow us to predict positive detection of the host in all filters and exposure times but are notably fainter for J2255+0251.These images are smoothed prior to imaging.

Figure 11 .
Figure 11.Same as Figure 10 but for a mock 10ks exposure of the F356W and F150W filters which are imaged with 3.1ks exposure times in Ding et al.(2023).This filter and exposure time is not included in followup proposals but shows a prediction for a longer exposure time with Figure10.Color bar units are MJy/str.These images are smoothed prior to imaging. .

Table 2 .
Relevant parameters for the EIGER quasars with host galaxies detected as measured or determined from Yue et al. (2023), or in their original discovery papers: Bañados et al. (2016), Yang et al. (2021), Farina et al. (2022), D'Odorico et al. (2023), and Mazzucchelli et al. (2023).All black hole errors are MCMC random errors and do not account for the usual 0.4 dex (Vestergaard & Peterson 2006) systematic errors from the scatter in the empirical relation used in Yue et al. (2023).

Table 3 .
The Ding et al. (2023)g galaxy magnitudes and ratios between galaxy and total (quasar + host) fluxes measured in SynthObs.The observed magnitudes and flux ratios are taken from Extended Data Table1inDing et al. (2023)using errors on the found using the dispersion on the results from five different PSF models.