Characterization of JWST NIRCam PSFs and Implications for AGN+Host Image Decomposition

We present a detailed analysis of the point spread function (PSF) of JWST NIRCam imaging in eight filters: F070W, F115W, F150W, F200W, F277W, F356W, F444W, and F480M, using publicly available data. Spatial variations in the PSF FWHM generally decrease with wavelength: the maximum and RMS fractional variations are $\sim20\%$ and $5\%$ in F070W, reduced to $\sim3\%$ and $0.6\%$ in F444W. We compare three commonly-used methods (SWarp, photutils, and PSFEx) to construct model PSFs and conclude that PSFEx delivers the best performance. Using simulated images of broad-line AGNs, we evaluate the impact of PSF mismatches on the recoverability of host galaxy properties. Host fluxes are generally overestimated when adopting mismatched PSF models, with larger overestimation for more AGN-dominated systems. Broader PSFs tend to produce less concentrated hosts while narrower PSFs tend to produce more concentrated and compact hosts. Systematic uncertainties in host measurements from PSF and model mismatches are generally larger than the formal fitting uncertainties for high signal-to-noise ratio data. Image decomposition can also lead to an artificial offset between the AGN and host centroids, which is common (e.g., $>1\sigma$ [$3\sigma$] detection in $\sim 80%$ [$\sim 20-30\%$] of systems), and scales with the mean host surface brightness. Near the surface brightness limit, this artificial offset can reach as large as $\sim80\%$, $26\%$, and $7\%$ of $R_e$ in systems with $R_e=$0.12", 0.48", and 1.92", respectively. We demonstrate our PSF construction and image decomposition methods with an example broad-line quasar at $z=1.646$ in the CEERS field.


INTRODUCTION
Tight correlations between the mass of supermassive black holes (BHs) and properties of their host galaxies, such as the stellar velocity dispersion (e.g., Gebhardt et al. 2000;Ferrarese & Merritt 2000) and the luminosity/mass of the bulge (e.g., Kormendy & Richstone 1995;Magorrian et al. 1998), suggest that they may coevolve with each other (e.g., Kormendy & Ho 2013).The energy feedback from rapidly accreting active galactic nuclei (AGNs) is widely proposed to regulate the growth and evolution of itself and its host galaxy (e.g., Somerville et al. 2008;Weinberger et al. 2018;Davé et al. 2019).
Due to the difficulties of measuring velocity dispersion from stellar absorption features, the studies of BH mass mingyang@illinois.edu(M BH ) and host property relations at high redshift mainly focus on the stellar mass (M * ) from broadband images or the dynamical mass as a surrogate from gas kinematics (e.g., Peng et al. 2006;Neeleman et al. 2021).In unobscured broad-line AGNs, M BH can be estimated using emission from the broad-line region using the reverberation mapping technique or single-epoch methods (e.g., Peterson et al. 2004;Greene & Ho 2005;Wang et al. 2009;Shen et al. 2015).However, the emission from the accretion disk can easily dominate the total observed emission at rest-frame ultraviolet and optical wavelengths, hampering the measurements of host properties.Fortunately, image decomposition, which makes use of the different spatial extent of the AGN (unresolved) and its host galaxy (can be resolved depending on its size, redshift, and telescope), offers a feasible way to extract the emission from the host galaxy.
At z 2, even the reddest filter F160W of the Wide Field Camera 3 (WFC3) onboard the Hubble Space Telescope (HST) can only probe host emission at rest-frame 5000 Å.At shorter wavelength, the AGN-to-host contrast is higher, rendering it more difficult to decompose host emission from the AGN.In addition, at rest-frame 5000 Å, host emission is dominated by the young stellar population and becomes insensitive to the total stellar mass.Meanwhile, emission from the host galaxy could be significantly blended with that of the AGN with ∼kpc (∼0.2) resolution.The unprecedented resolution (0. 03-0.16), high sensitivity, and near-infrared (NIR) coverage up to ∼ 5 µm of the Near Infrared Camera (NIRCam) onboard the James Webb Space Telescope (JWST; Gardner et al. 2006) make it possible to detect rest-frame 0.5 − 1.6µm emission of AGN host galaxies at z > 2 for the first time.For example, with a full-width-at-half-maximum (FWHM) of the point-spreadfunction (PSF) of ∼0.09 in the F277W filter, we can probe rest-frame ∼ 0.93 µm emission from z = 2 galaxies at a spatial resolution of 0.75 kpc.The significantly improved resolution and rest-frame NIR coverage are essential to studying the properties of host galaxies of high-z broad-line AGNs.
No instrument has the perfect PSF that does not vary with time and position on the focal plane, even for the NIRCam of JWST.Recent work has shown that the NIRCam PSF has strong spatial (up to ∼ 15 − 20%) and temporal variations (∼ 3 − 4%) based on images of individual exposures (Nardiello et al. 2022).However, as AGN-host decomposition and most of the measurements are performed on images after combining individual exposures, it is important to analyze and quantify the properties of PSF for the final combined images.In this work, we compare the performance of commonly-used methods for constructing model PSFs and utilize the best PSF models to characterize the properties of the NIRCam PSF and their spatial variation.We then quantify the effect of PSF mismatch on the derived properties of AGN host galaxies from image decomposition using mock images of AGNs.As a demonstration of our overall methodology of PSF modeling and AGN-host decomposition, we present host galaxy measurements for a z = 1.646 broadline AGN.
For this purpose, we use NIRCam imaging observations pointing in the south continuous viewing zone (S-CVZ; R.A.= 90.0000 • and Decl= −66.5607 • , i.e., the south ecliptic pole).This particular field has the advantage of having a much higher star density compared with other high-galacticlatitude fields, but much less source blending compared to observations of star clusters.These observations provide a large sample of stars from which we construct model PSFs and quantify their characteristics.
This paper is structured as follows.Section 2 describes the data, data reduction, and PSF construction methods.Sec-tion 3 compares the performance of different PSF construction methods and quantify the properties of PSF and their spatial variations.Section 4 studies the effects of PSF mismatch on the recovered host properties.Section 5 applies our methodology to a real AGN.Our main conclusions are summarized in Section 6.We adopt a cosmology with H 0 = 70 km s −1 Mpc −1 , Ω m = 0.3, and Ω Λ = 0.7 and a Chabrier (2003) initial mass function (IMF).

Data
We use public NIRCam imaging observations from the Cycle 0 calibration program "Stray Light Pointed Model Correlation" (PI: Erin Smith, ID=1448) pointing at the S-CVZ field on 6th May 2022.The observations used Module B of NIRCam with the FULL subarray setting, covering a field-of-view (FoV) of ∼2. 2 × 2. 2. Module B consists of 4 short wavelength (SW; 0.6-2.3µm)detectors and 1 long wavelength (LW; 2.4-5.0µm)detector.SW detectors are separated by ∼ 4 − 5 , leaving gaps of cross shape at the center.Observations were taken in four filter pairs (SW+LW simultaneously) of three integrations each, with a total exposure time of 450.944 s for each filter pair: F070W+F277W, F115W+F356W, F150W+F444W, and F200W+F480M.A three-point dither pattern was used to optimize the subpixel sampling and to mitigate bad detector pixels.In particular, Mid-Infrared Instrument (MIRI) imaging was taken simultaneously with NIRCam imaging as coordinated parallel observation for this program.Therefore, a 3-POINT-WITH-MIRIF770W dither pattern is adopted due to the requirement of a larger dither step by MIRI imaging.

Data Reduction
We use version 1.8.5 of the jwst pipeline with the Calibration Reference Data System (CRDS) version of 1039, with custom steps for data reduction.The detailed procedures are as follows: (1) We reduce the uncalibrated ( uncal.fitsfiles) raw data using the Stage 1 pipeline calwebb detector1.Detector1Pipeline, which applies basic detection-level corrections to each exposure, by adopting the default configuration with the exception of turning on the snowball correction.(2) We then apply 1/f noise (horizontal and vertical striping patterns) subtraction to the products of stage 1 pipeline (i.e., count rate data; rate.fitsfiles) using scripts developed by the Cosmic Evolution Early Release Science Survey (CEERS; ERS 1345, PI: Steven Finkelstein) team (Bagley et al. 2022).(3) We run Stage 2 pipeline calwebb image2.Image2Pipeline with default parameters, which involves steps including wcs assignment, flat-fielding, and photometric calibration and returns fully calibrated individual exposures ( cal.fits files).( 4 We perform 2-dimensional background subtraction on calibrated data after masking sources and bad pixels with a box size of 50 pixels using SExtractorBackground implemented in photutils (Bradley et al. 2022).Astrometry is corrected using coordinates of sources from the GAIA data release 3 (DR3; Gaia Collaboration et al. 2022) after accounting for their proper motions.(5) We combine individual exposures into mosaics using resample step in Stage 3 pipeline calwebb image3.Image3Pipeline.The output mosaics are drizzled to a common pixel scale of 0. 03 (30 mas) pixel −1 without shrinking the input pixels (pixfrac= 1) in all the filters following the choice by the CEERS team (Bagley et al. 2022).All other pipeline parameters are set with their default values.The final round of background subtraction is performed to remove the remaining background (if any) in the mosaics.All the JWST data used in this paper can be found in MAST: 10.17909/6btv-br09.

PSF Construction
We use SExtractor (Bertin & Arnouts 1996) to perform source detection and to select candidate PSF sources that are of high signal-to-noise ratio (SNR; SNR WIN > 100), non-blended (FLAGS < 2), non-irregular (ELONGATION < 1.5), point-like (CLASS STAR > 0.8), and without bad pixels (IMAFLAGS ISO = 0).We then construct a master point source catalog by selecting objects that are detected in at least two adjacent filters.The final point source catalog consists of 1338 objects with CLASS STAR 0.96 for 90% of them.For each object in each filter, we generate cutout image with a size of 135 pixels (∼ 4 ) for SW images and 203 pixels (∼ 6 ) for LW images, which encircle at least 98% and 97% of the light (Rigby et al. 2023), respectively.Given the extremely high stellar density in this field (i.e., proximity to the Large Magellanic Cloud), the vast majority of these point sources are stars.
For high-resolution imaging analyses with HST or JWST data, it is customary to build empirical PSF models using a library of PSF sources.The major benefit of these PSF models is that they can mitigate the noise in the low surface brightness wings of the PSF, compared with PSF templates directly based on observed stars.We construct PSF models using two different methods: one by building a pixel-based model and the other by stacking PSF star images.Following Zhuang & Ho (2022), we use the Python package photutils (Bradley et al. 2022), which is based on the algorithm of Anderson & King (2000), and the software PSFEx (Bertin 2011) for the first method, and SWarp (Bertin et al. 2002) for the second method 1 .We adopt default configurations for these packages, in particular a quadratic smoothing kernel for photutils, a PIXEL basis type for PSFEx, and a lanczos3 resampling method for SWarp.
We construct two types of PSF models for each method, one using all of the point sources across the entire FoV (global), and the other using point sources inside a particular region (region) for photutils and SWarp to probe variations of the PSF across the FoV.As illustrated in Figure 1, we divide the final mosaic into 36 regions.One of the advantages of PSFEx is that it supports modeling of spatial variations of the PSF model across the FoV as a function of pixel coordinates (x and y) in the form of a polynomial.In this work, we adopt a third-order polynomial with ten vector images (constant, x, x 2 , x 3 , y, xy, x 2 y, y 2 , xy 2 , and y 3 ).For PSFEx, we construct the PSF model for each region using the central position of that region.
The PSF models are constructed with an oversampling=2, i.e., the pixel size of the PSF model is 1/2 of that of the input image (30 mas).The final image size of PSF models is 265×265 pixels for SW filters and 401×401 pixels for LW filters.We verify that PSF models with oversampling=1 produce the same trends shown in Section 3. The only exception is for photutils, oversampling=1 produces systematically larger FWHMs due to kernel smoothing, which occasionally fails, and in some cases produces unsatisfactory results.
1 we ignore the reproject package as its reproject exact function is currently known to have precision issues for images with resolutions < 0. 05.   6) and ( 8): Minimum, median, and maximum values of PSF FWHM and q among 36 regions.Col. ( 7) and ( 9): 3σ-clipped standard deviations of PSF FWHM and q among 36 regions, with fractional RMS in the parentheses (measured with respect to the median value across all 36 regions).

PSF Characteristics
We measure the shape of each PSF model by fitting a 2dimensional Gaussian model to its core, with a fitting size of 2×N +1, where N marks the index of pixel that is just larger than the half of the PSF FWHM.For region PSF models, we use the fitting size determined from the global ones.The geometric properties of the PSF models are presented in Table 1.Three methods generally produce consistent results on the geometry of the PSF models.PSF models constructed using photutils tend to have larger FWHMs in SW filters, likely due to the more significant kernel smoothing effect in less oversampled PSF models with small FWHMs (∼ 2 pixels).
Using the global PSF model as the typical representation for each filter, we find that PSFs in all SW and LW filters are critically-or well-sampled (FWHM 2 pixels), given the drizzle parameters adopted in this paper (pixel scale=0.03 pixel −1 and pixfrac=1).Note that the FWHMs reported here are larger compared with those reported in the JWST user documentation2 , especially toward shorter wavelength filters.As shown in Bagley et al. (2022), the recovered FWHM of PSF in the final mosaic is affected by the output pixel scale and pixfrac, with larger FWHMs from larger output pixel scale and larger pixfrac, due to increased pixel-to-pixel correlations.The PSF cores exhibit a high degree of axial-symmetry, with axis ratio (major-axis/minoraxis) q 0.97 in all filters, except for the F070W filter (q ≈ 0.9).
We find clear spatial variations of the PSF across the FoV, as indicated by the spread of FWHMs of region PSF models from all three methods (Table 1).The degree of spatial variation depends on wavelength, with maximum fractional difference (maximum/minimum−1) and 3σ-clipped standard deviation (RMS) of the PSF FWHM decreasing from ∼ 20% and 5% in the F070W filter to ∼ 3% and 0.6% in the F444W filter.This level of variations is similar to that reported by Nardiello et al. (2022).Using F070W as an example, Figure 2 and Figure 3 show the region PSF model images across the FoV and their normalized radial surface brightness profiles, respectively.There are significant radial variations up to 100% at 60 − 90 mas (2 − 3× FWHM) from the center.The variations persist toward larger radii and are dominated by random noise at r 1 due to background fluctuations.Thanks to the simul-taneous modeling of point sources across the FoV, PSF models constructed using PSFEx have much higher SNR at large radii compared with the other two methods, with less contamination from background noise and low-level, small-scale gradients that are not accounted for by background subtraction.Together with better performance of PSF models constructed from PSFEx in terms of modeling the light profiles of point sources, as well as measuring their magnitudes (Section 3.2), we use PSF models constructed using the PSFEx method as the reference in the following analyses.
Figure 4 shows the spatial variations of FWHM, axis ratio, and position angle of the PSF models in eight filters across the FoV of NIRCam.We find that the PSF models with the broadest and narrowest profiles are preferentially located toward the edge of the FoV.Different filters show distinct patterns: a horizontal pattern in F070W, a vertical pattern in F115W, F277W, F444W, and F480M, and a diagonal pattern in F150W, F200W, and F356W.These patterns could be affected by stochastic aliasing of the PSF during the process of drizzling, which is more significant for SW filters where the output pixel scale is close to the input value with experience from HST (Rhodes et al. 2007).Axis ratio and position angle exhibit somewhat less spatial variation, with extreme values more commonly found toward the corners.The variation in position angle is partially due to the large uncertainty of the position angle measurement at large axis ratios.
It is possible that the spatial variation of the PSF shown in this work is somewhat underestimated due to the smoothing from the adoption of the third order polynomial function for the PSFEx (Section 2.3).While these spatial variations of the NIRCam PSF are modest in most filters (e.g., fractional RMS of a few percent), they may impact the detection and measurements of host galaxies in unobscured AGNs where the light is dominated by the centrarl AGN.Quantifying the effects of PSF mismatches on AGN+host image decomposition is one of the main motivations of this work (Section 4).

The Effect of PSF Spatial Variations on the Measurements of Point Sources
To investigate the effect of PSF spatial variations on the measurements of point sources, we fit their images with gloabl and the corresponding region PSF models constructed from the three methods using galfit (Peng et al. 2002(Peng et al. , 2010)).In particular, as PSFEx can provide a PSF model at any location based on interpolation, we also include a local PSF model for each source using its position on the mosaic.

Light Distribution of Point Sources
Figure 5 shows the normalized probability density distribution of reduced χ 2 , which is an indicator of goodness of the fit in terms of the light profile.We divide the sources into 4 categories: (a) 10 brightest sources, (b) sources with magnitude between 0th and 5th percentiles, (c) sources between 5th and 16th percentiles, and (d) sources between 16th and 100th percentiles.We also limit our tests to regions that contain at least 5 point sources used for PSF construction to avoid the problem of "fitting one object with itself".All PSF models provide overall good fits (χ 2 ν ≈ 1 − 2) to the majority (84%) of the sources (Figure 5d).However, the performance of different PSF models tends to vary as source becomes brighter.For three global PSF models, PSFEx performs slightly better than SWarp, and both are noticeably better than photutils.For region PSF models, photutils provides the least satisfying results among 7 types of PSF models, likely due to the relatively low quality PSF models built based on a small number of objects (see Figure 2 for examples of PSF models using 11-26 objects in the F070W filter).region PSF models from SWarp perform similarly well as global PSF models in F070W, F115W, and F150W, but slightly worse in F200W and the three LW filters.re-10 -2 10 -1 10 0 Semi-major axis (arcsec)   gion and local PSF models from PSFEx have slightly improved performances compared to global PSFEx in F070W and F115W, likely due to the relatively larger PSF spatial variations in these bands that are properly accounted for with the region and local models.Based on these results, we conclude that the three PSF models (global, region, and local) from PSFEx provide the best performance among 7 PSF models in describing the light distribution of point sources.However, we note that although PSFEx constructs PSF models by modeling multiple point sources simultaneously, its high order vector images may still have relatively poor SNR due to limited number of objects located close to the edges (Figure 2).Therefore, we check the performance of region PSF models from PSFEx using the median χ 2 ν for all point sources in each region.We find that region PSF models for regions that are close to the edges perform equally well as those for inner regions, albeit with lower SNR.This demonstrates the overall robustness and reliability of PSFEx in modeling the spatial variations of the PSF.

PSF Magnitude of Point Sources
Figure 6 compares the magnitude of point sources measured using different PSF models.We use PSF magnitude (m PSF ) measured using the global PSFEx PSF model as the reference, given its high SNR in the PSF wings and overall good performance in describing the light distribution of point sources (Section 3.2.1).We find that m PSF measurements from three PSFEx PSF models are very consistent with each other, with |∆m PSF | 0.02 mag in systematic offset and ≤ 0.05 mag in random scatter across all the filters.Both PSF spatial variation and larger noise toward the fainter end of sources could contribute to the observed total scatter.The largest scatter observed in the F070W filter is most likely driven by its significant PSF spatial variations (Table 1).The systematic differences, although rather mild, could result from the uneven distribution of sources in the FoV.In summary, our results suggest that the spatial variations of the PSF have a mild ( 0.05 mag or 5%) effect on measuring the magnitude of point sources.
For global PSF models from the other two methods, we find that they are in general good agreements, albeit with consistently (scatter 0.02 mag) larger systematic differences with respective to the reference magnitudes.Noticeably, the offsets are higher for photutils PSF models (|∆m PSF | 0.19 mag) than SWarp PSF models (|∆m PSF | 0.07 mag).The consistence is much worse for region PSF models from SWarp (|∆m PSF | up to 0.09 mag and scatter up to 0.14 mag) and photutils (|∆m PSF | up to 0.15 mag and scatter up to 0.28 mag).We find that the systematic magnitude offset is mainly driven by inadequate background subtraction.As we are using a large cutout image size (135×135 pixels for SW filters and 203×203 pixels for LW filters) for PSF construction, the outskirt of the majority of the sources is dominated by background and thus could be affected by imperfect background subtraction.To test this hypothesis, we reduce the size of the image of PSF models by half and re-normalize it by subtracting the median value of pixels close to the edges (5 pixel wide).As shown in Figure 7, the substantial magnitude offsets (|∆m PSF | 0.19) for global PSF models are significantly reduced to |∆m PSF | 0.03 mag for both SWarp and photutils in all filters except F480M.The residual systematic offsets could arise from subtle differences in FWHM (as a proxy of the PSF model profile) and minor over-(under-)estimation of background noise.Similarly, region PSF models, in particular those of photutils, are more likely to suffer from the factors mentioned above and exhibit a larger  degree of inconsistency, mainly due to fewer sources available for PSF construction in individual regions.The small offsets (≤ 0.01 mag) between PSFEx PSF models of two sizes in all the filters suggest that this method is insensitive to background variations and subtraction.

Conclusions on PSF Characterization and Recommendations
Our results show that PSF models constructed using point sources in the drizzled mosaics with a pixel scale of 0. 03 pixel −1 and pixfrac=1 in eight NIRCam filters (F070W, F115W, F150W, F200W, F277W, F356W, F444W, and F480M) have FWHM 2 pixels.The NIRCam PSF has significant spatial variations across the FoV, with maximum and RMS variations of FWHM decreasing from ∼ 20% and 5% in F070W to ∼ 5% and 0.6% in F444W.Different filters have diverse variation patterns with extreme FWHMs being predominately found close to the edges.Neglecting PSF spatial variations would have minor impact on the magnitude of point sources across the FoV (≤ 0.02 mag systematic offset and ≤ 0.05 mag scatter), depending on the required photometric accuracy.
We recommend the use of PSFEx to construct the PSF model, as well as modeling its spatial variations.It can deliver PSF models with higher SNR toward larger radii, with better performance in recovering the light profile and magnitude of point sources compared to the other two methods.If a sufficient number of point sources are available, it is recommended to construct region or local PSF models.Otherwise, adopting the global PSF model still provides satisfactory results under most circumstances, with an average 0.02 mag systematic offset and 0.05 mag random scatter.For the use of the other two methods (SWarp and photutils) to construct PSF models, adopting a smaller PSF model size would significantly mitigate the effects of background fluctuations and imperfect background subtraction.

Future Work
As our analysis is based on combined mosaics, the characteristics of the PSF and its spatial variations presented here are subject to the specific configuration of observations.Observations with various dither patterns and steps, number of exposures, and drizzling parameters (e.g., output pixel size, pixfrac, and flux distribution kernel) are required to fully explore their effects on the geometry and spatial variation of the PSF.

IMPLICATIONS ON THE AGN-HOST IMAGE DECOMPOSITION
For our work, the main motivation of PSF characterization of NIRCam imaging is to assess the reliability of host measurements in unobscured AGNs, or their high-luminosity counterparts, quasars, using imaging decomposition techniques.Below we apply our methodology of PSF model construction discussed above to simulated AGN+host images, and investigate the recoverability of host properties.
To achieve this goal, we generate mock AGN+galaxy images with a central point source (AGN) and a Sérsic profile (host galaxy), convolved with the measured global PSF model from PSFEx.The centers of the AGN and the host are tied together in the decomposition with zero offset.These simulated systems, although oversimplified, provide critical information on the true uncertainties and biases of recovered host properties that are often underestimated from fitting (e.g.Häussler et al. 2007;van der Wel et al. 2014;Zhuang & Ho 2022).For simplicity, we do not consider more complex substructures of the host galaxy, which could lead to additional systematics (Zhuang & Ho submitted).
We generate mock systems across a broad range of parameters.We consider AGN magnitude (m AGN ) from 16.5 to 30.0 mag in steps of 0.5 mag; AGN-to-host ratio from 0.1 to 10 in steps of 0.2 dex, which corresponds to a step of 0.5 mag; effective radius (R e ) from 4 to 64 pixels (0. 12-1.92) in steps of 0.3 dex, which covers the range of ∼ 1 − 10 kpc at z = 1 − 7.5; Sérsic index n = 0.5, 1, 1.5, 2, 3, 4, and 6; axis ratio between major-and minor-axis q = 0.3, 0.6, and 0.9; a fixed position angle of 45 • .For the convenience of applying our results to other fields with different exposure times (hence different depths), we measure the 5σ point source depth and 3σ surface brightness (SB) depth, where σ is the standard deviation of 50 flux measurements randomly distributed in the background image within an aperture of radius = PSF FWHM (for point sources) and 0. 5 (for extended sources), respectively.The measurements for point sources are then corrected for flux loss outside the aperture to deter- mine the final point-source depth.The results are presented in Table 2.
In total, we generated 32,340 mock AGN+galaxy images in each of the 7 wide filters (except F480M) with a cutout size of 6 × R e to enclose most of the emission from the host galaxy.Two types of noise are added to the mock images: (1) Poisson noise for the source and (2) background noise.For the background noise, we use real images from which the PSF models are constructed.Our final background cutout contains no bright sources with only ∼ 0.3% of pixels having fluxes larger than 5× background standard deviation.The original cutout has a size of ∼ 180 pixel (5.4) and is padded in a 3 × 3 array to construct the final background image.
We fit images of mock AGNs with four sets of fitting scenarios using galfit: (1) the fiducial model with the same PSF model as the one used to generate the mock images (global PSFEx); during the decomposition, the centers of the AGN and Sérsic components are tied together without offset; (2) same as (1) but with a broader PSF model; (3) same as (1) but with a narrower PSF model; (4) same as (1) but allow the centers of AGN and Sérsic to offset from each other.For all four fitting scenarios, we set limits of fitting parameters: magnitude to 10-35 mag, R e to 0.5-100 pixels, n to 0.3-8, q to 0.1-0.99,and center position constrained to within ±1 pixel of its input value for the first three sets and ±3 pixels for the last set.The fitting results for all mock images are made available to the public in the github.
Fitting scenarios (1) and (4) will evaluate the recoverability of host properties as functions of input parameters assuming a perfect PSF model, but with different fitting methodologies (whether or not tie together the centers of the two components).Fitting scenarios (2) and (3) will evaluate additional uncertainties from PSF mismatches.When using broader or narrower PSF models, we adopt two region PSF models with the largest and smallest FWHMs to illustrate the maximum variations in each filter -they differ from the fiducial PSF FWHM by ∼ −7% and +12% in the F070W and ∼ −2% and +1% in the F444W filter.For illustration purposes, we show the results for host galaxies with 19 ≤ m Sérsic ≤ 27.5 mag in the F150W filter in the following part of Section 4. This filter displays a moderate level of spatial variations among all filters (Table 1), so that the results presented in Section 4.2 can be considered a typical case for a given filter.

Reliability of Host Measurements
Figure 8 compares the recovered parameters of mock AGNs with the input values, as a function of m Sérsic using the global PSF model from PSFEx in the F150W filter.This test represents the situation where the PSF model is perfect (i.e., the mock images are generated using the same PSF model).We find overall good consistency between the recovered parameters and the input values for bright host galaxies.As host galaxies become fainter, the recovered host properties tend to show larger measurement uncertainties, larger offset (scatter) compared with input values, and a larger fraction of objects have failed measurements (i.e., measurement/uncertainty< 3; open symbols).In addition, the consistency between re- covered host parameters (m Sérsic , R e , and n) and the input values degrades toward higher AGN-to-host ratio, while the recovery of m AGN becomes worse toward lower AGN fractions.No obvious differences are found among systems with different q values.As long as m Sérsic is successfully measured (measurement/uncertainty≥ 3), q can be estimated correctly as well, albeit with larger scatter toward higher AGNto-host ratio (e.g., fractional scatter increases by a factor of ∼ 3 at q = 0.6, from 3% for AGN-to-host ratio=0.1 to 10% for AGN-to-host ratio=10).
At fixed R e , we find systematic overestimation of the host flux (smaller m Sérsic ) due to a combination of overestimated R e and underestimated n. m AGN tends to be underestimated (flux overestimated) when host is more concentrated (larger n).At fixed n, both small and large R e can lead to slight underestimation of host flux and larger scatter in all four parameters presented here.
Among the three host parameters, m Sérsic has the highest fraction of successful measurements, followed by R e and n.The exact success rate depends on the AGN-to-host ratio and host structure.Galaxies with higher n and both small and large R e are more likely to have unsuccessful measurements.m Sérsic and R e in objects with low AGN-to-host ratios, small n and moderate R e can be successfully measured down to the 5σ point source depth.However, no successful measurement can be obtained at m Sérsic > 25 mag (∼ 2 mag above the 5σ point source depth) for n, with little dependence on the AGN-to-host ratio.
The deviations from input values (∆) are larger than the formal measurement uncertainties (error bars in Figure 8) for to-host ratio=10, with extreme values as high as ∼ 100.This urges the need of generating object-tailored mock AGNs with input parameters from the best fit to evaluate systematic biases and to obtain more realistic measurement uncertainties.

The Effect of PSF Mismatch on Host Measurements
We now investigate the impact of PSF mismatch on the recovery of host parameters.Figure 9 and 10 show the results from a broader (3.6% broader in the FWHM selected from region 6) and a narrower (4.2% narrower in the FWHM selected from region 31) PSF model in the F150W filter, respectively.
Recovered parameters from the broader and narrower PSF models generally follow the same trends as those using the perfect PSF model for mock AGNs with AGN-to-host ra-tio=0.1 and 1, albeit with slightly larger scatter and offsets.Significantly larger offsets are observed in AGN-dominated systems (AGN-to-host ratio=10).
We find that although the total fluxes of host galaxies are more likely to be overestimated (smaller m Sérsic ) for both mismatched PSF models, the reasons are different.The  broader PSF model tends to underestimate n, predicting more extended host galaxies.This behavior is expected since the emission from the central part of the host is mistaken as part of the unresolved emission from the AGN with the broader PSF model, causing the fitted host to be more extended.Since Sérsic models with larger n contain a larger fraction of emission outside R e , the more extended emission at large radii coupled with the true n from the host galaxy will upscale the overall flux from the Sérsic model, when the fitted n is underestimated using a broader PSF model.On the other hand, the narrower PSF model tends to underestimate R e and overestimate n, predicting more concentrated and smaller host galaxies.Therefore, with the narrower PSF model, the emission from the AGN is split into contributions from a narrower point component and an extended component.Consequently, AGN flux is underestimated and host flux is overestimated for the narrower PSF model.For systems with AGN-to-host=10 (i.e., AGNdominated systems), the underestimation of AGN flux can reach as much as ∼ 0.6 mag, and the overestimation of host flux can reach as much as 1.8 mag!Moreover, in systems with AGN-to-host=10 and successful m Sérsic measurements from the narrower PSF model, the measured q values are generally around ∼ 0.6 regardless of their input values, making it difficult to estimate the inclination angle of the host.
The structure of the host galaxy plays a more important role on the host parameter recovery when the PSF model is mismatched.Measurements are more likely to fail (measurement/uncertainty< 3) in small host galaxies.The broader PSF model tends to predict smaller m Sérsic , smaller R e , and higher n toward larger host galaxies (larger R e ) and lower n toward more concentrated host galaxies (larger n).The narrower PSF model follow the same trends as the broader PSF model except for the relation of n and galaxy size (no meaningful trend due to no successful measurements for hosts with R e =0. 12 and 1. 92).
For the measurement success rate, the broader PSF model generally produce similar results as the fiducial PSF model.The only obvious difference is for the host galaxies with small R e .The measurement success rate is much lower for the narrower PSF model.For example, for host galaxies with AGN-to-host ratio=10, almost no successful measurements can be obtained for R e at m Sérsic 21.5 mag (∼ 5.5 mag above the point source depth).
To better facilitate a detailed look at the effect of the AGNto-host ratio, Figure 11 shows the differences of the parameters as a function of AGN-to-host ratio in mock AGNs with m Sérsic = 24 mag.Results from the broader PSF model are similar to those from the fiducial PSF model, except for n where gradual offset and scatter start to appear at AGN-tohost ratio 1.However, the narrower PSF model would lead to noticeable offsets and scatters in all five parameters.The AGN-to-host ratio tends to be underestimated in AGNdominated systems (AGN-to-host ratio 1), leading to a higher host fraction.Almost no successful measurements of R e and n can be obtained at AGN-to-host ratio 1. m Sérsic can be significantly underestimated (hence host flux overestimated) at AGN-to-host ratio 1.
The goodness of fit can be assessed with the reduced χ 2 (χ 2 ν ). Figure 12a presents the median χ 2 ν from galfit decomposition of our mock images with successful m Sérsic measurements.For systems with AGN-to-host ratio=0.1 or 1, χ 2 ν from galfit is generally 1 across the full m Sérsic range (m Sérsic > 19 mag).Although results with the fiducial PSF model often have the lowest absolute values of χ 2 ν , we are unable to tell which PSF model performs the best, given the overall good performance of all three PSF models.On the other hand, when the host and the AGN are both well above the flux limit (e.g., m Sérsic < 21 mag or ∼ 6 mag above the PS depth, and AGN-to-host ratio=10), the χ 2 ν metric is a good parameter to select the best PSF model, since the S/N of both the AGN and host components is extremely high.
However, as the whole system approaches the imaging depth, χ 2 ν values from the broader and narrower PSF models start to become comparable to that from the fiducial PSF model.At m Sérsic 23 mag the χ 2 ν values from three PSF models are all statistically acceptable for the fit.In fact, in some cases, the mismatched PSF models produce smaller χ 2 ν values than the fiducial PSF model.We find that the distributions of χ 2 ν of three PSF models in systems with 24 ≤ (m Sérsic /mag) ≤ 27.5 are quite similar (Figure 12b), albeit with differences in the number of objects with successful m Sérsic measurements for different AGN-to-host ratios.Real systems have light profiles that are more complicated than a PSF + a single Sérsic due to the presence of substructures (such as bulge, bar, lens, dust lane, and star-forming clumps), and the discriminating power of χ 2 ν may be further weakened.These tests on χ 2 ν suggest that unless both the AGN and the host are much brighter than the imaging depth, a smaller χ 2 ν does not necessarily mean the PSF model is correct.In fact, for most imaging applications where the host (or the full system) is often near the detection limit, PSF mismatch is a realistic concern and cannot be easily identified with the χ 2 ν statistic.The systematic biases in the derived parameters due to PSF mismatches, as we demonstrated here, must be taken into account in the total error budget.
To summarize, PSF mismatch has significant effects on the derived host (and AGN) parameters from AGN-host image decomposition.Adopting a broader PSF model would lead to overestimated host fluxes with more extended host structures.Adopting a narrower PSF model could also lead to overestimated host fluxes but with more concentrated and smaller host structures.For a given object, adopting a narrower PSF model would introduce much worse systematic biases (e.g., underestimate m Sérsic by ∼ 1.6 mag) and lower measurement success rate compared to adopting a broader PSF model.The systematic biases caused by adopting a narrower PSF model are noticeable even in host dominated systems.

Reliability of AGN-Host Offset Measurement
Figure 13 compares the differences of center measurements versus the mean surface brightness of the Sérsic component (µ S érsic ) between the fiducial fitting scenario and scenario 4, where in the latter case the centers of the AGN and host components are untied.µ S érsic is defined as m Sérsic +2.5 log(2qπR 2 e ).The factor of 2 is to account for the fact that half of the total emission is within the effective radius.
We find that one can easily measure an artificial offset between the AGN and its host galaxy centroid, mainly due to incorrect center measurements of host galaxies.The artificial offset generally follows the direction of the position angle of the host (i.e., 45 • ), randomly distributed on either side.The amount of offset increases rapidly with decreasing surface brightness and hits the boundary of 3 × √ 2 pixel when the surface brightness of the host drops below the SB detection limit.The trends are similar in host galaxies with different structures.When considering all the sources with 19 ≤m Sérsic ≤ 27.5 mag, we find that offset measurements in 59-64%, 28-40%, and 11-21% of objects are larger than 1×, 2×, and 3× the uncertainty reported by galfit, respectively.These percentage ranges encompass the AGN-tohost ratio range in our mock data (0.1-10).Excluding objects with µ S érsic below the SB detection limit (associated with very large measurement uncertainty), offsets in 76-82% (15-29%) of the objects are larger than the reported 1σ (3σ) uncertainty.This fraction decreases toward larger AGN-tohost ratio due to larger position measurement uncertainties for both the AGN and its host.
We note that the results above are only for idealized mock AGNs.With the presence of galactic substructures (including both axial symmetric and non-axial symmetric such as bulge, nuclear bar, spiral arms, lens, and dust lane) at the center of the galaxy, we expect more complicated trends in the measured artificial offset.The offsets can be as large as ∼ 80%, 26%, and 7% of R e in objects with R e =0. 12 (4 pixels), 0. 48 (16 pixels), and 1. 92 (64 pixels) when approaching the SB detection limit, respectively.
For completeness, no obvious differences in the global trends and scatter of the four parameters (m Sérsic , m AGN , R e , and n) are found whether the centers of AGN and host are tied in the fit or not.

APPLICATIONS OF OUR METHODOLOGY
As a demonstration of applying our method of optimal PSF modeling to study the properties of host galaxies of broadline AGNs, we select a broad-line AGN SDSS1420+5300A (R.A. = 215.02333• , Decl = 53.01020• ) at z = 1.646, covered by Module A of NIRCam imaging in pointing 1 of the CEERS fields.This AGN has a BH mass (M BH ) of M BH = 10 9.00±0.16M from reverberation mapping (Grier et al. 2019).

PSF Model Construction
CEERS NIRCam observations have seven filters: F115W, F150W, F200W, F277W, F356W, F410M, and F444W.Six wide filters are considered here for the AGN-host decomposition.The early NIRCam observations are taken on June 21-22, 2022 in a coordinated parallel mode with MIRI imaging.We adopt the version 0.5 NIRCam imaging data product from CEERS team (Bagley et al. 2022), which covers all four pointings (1, 2, 3, and 6) of June 2022 observations.Two types of observing configurations are adopted due to the use of different MIRI filters.For pointings 1 and 2, dither patterns of NIRCam are tied to the corresponding MIRI filters that are observed simultaneously, while for pointings 3 and 6, the same 3-POINT-MIRIF770W-WITHNIRCam dither pattern is adopted.Therefore, we construct PSF models for this object using only point sources in the Module A of pointings 1 and 2. To model the spatial variations of the PSF, only a first order polynomial (linear) dependence provides usable PSF models due to the very limited number (∼ 20) of point sources in the field.We construct three PSF models to investigate the differences of host properties: (1) a local PSF model constructed at the location of the AGN in the image; (2) a global PSF model; (3) a PSF model from (2) but broadened by a circular Gaussian kernel, so that the FWHM of the output PSF model is 1σ broader.Here σ is the RMS fractional variation presented in Col. ( 7) of Table 1.
There are several reasons for testing three PSF models: (1) the linear dependence of PSF model constructed locally at the AGN may not be accurate enough; (2) the global PSF model is representative of the median PSF model across the FoV of NIRCam; (3) in case the global PSF model underestimates the true PSF of the target, we slightly broaden it to mitigate the biases in the recovery of host properties, since a narrower fitting PSF can have severe consequences in the host recovery (Section 4.2).We find that the local PSF models are generally narrower compared to the global ones except in the F150W filter.The FWHMs of our PSF models are slightly smaller than those reported in Finkelstein et al. (2023), which are constructed by stacking stars across all four pointings.We make our global PSF models for CEERS fields publicly available to the community in github.

AGN-host Image Decomposition
We employ galfitm (Häußler et al. 2013;Vika et al. 2013) to perform the AGN-host image decomposition.galfitm is a multiwavelength version of galfit, which can model images from multiple filters simultaneously, taking into account the wavelength-dependent galaxy structure.Works on nearby galaxies find that the Sérsic index n and R e tend to vary systematically with wavelength, with increasing n and decreasing R e from rest-frame 0.3 to 2 µm (e.g., Kelvin et al. 2012).This change of structure as a function of wavelength could result from variations in the spatial distributions of stellar populations, metallicity, and dust attenuation.
We consider two sets of fitting configurations for the decomposition: (1) A PSF model for the AGN and a Sérsic model for the host (single Sérsic; SS); (2) A PSF model for the AGN, two Sérsic models for the host (double Sérsic; DS).For the DS fitting, the first Sérsic model is for the elongated emission from the bar/lens-like component revealed in the residual image from fitting configuration 1, and the second Sérsic model is for the emission from the disk-like component.Combining three options for PSF models, we have six sets of fits, denoted as SS1, SS2, SS3, DS1, DS2, DS3.We note that adding an additional bulge component is not evident.Our tests with a PSF+bulge+bar/lens+disk model suggests a bulge with R e ≈ 1 pixel, whose emission is blended and dominated by the central point source.
We use a cutout size of 201 × 201 pixels (∼6 ), which is large enough to enclose most of the emission from our target.A close extended companion that is blended with our target is also fitted simultaneously with a Sérsic model.All other sources in the cutout are masked.The centers of the AGN and the host are tied together and allowed to vary independently in different filters to account for residual astrometry mismatches ( 0.2 pixel).For the SS fitting configuration, we use the same values for the ellipticity and position angle parameters across the wavelength, set m Sérsic and m AGN free to vary (fifth-order polynomial), and allow the n and R e to vary quadratically with wavelength, following Häußler et al. (2013).n is constrained to be in the range of 0.3 to 7. For the DS fitting configuration, we do not allow n and R e to vary with wavelength for both Sérsic models following Häußler et al. (2022), with the assumption that no significant color gradient is present within each component.
We set initial guesses of R e = 10 pixel, n = 3, and m Sérsic and m AGN to 1 mag fainter than the total magnitude for the SS configuration.For the DS configuration, we set R e = 8 pixel, n = 2 for the first Sérsic model, R e = 20 pixel, n = 1 for the second Sérsic model as initial guesses.The two Sérsic models have the same initial m Sérsic , which is 0.75 mag larger than the best-fit m Sérsic from the SS fit.We perform the decomposition twice for both SS and DS configurations with the best-fit parameters of the first run as the initial guess for the second run.
As mentioned in Section 4, true uncertainties of the parameters are often larger than the formal uncertainties reported by the fitting code.We follow the procedures described above to generate mock AGNs using the best-fit parameters.For the sky background, we select a "blank" region in the mosaic near our target in each filter.The background cutout has the same size as the AGN and does not contain any visible source.We further remove residual sky pedestal to ensure that the 3σ-clipped median is exactly zero.For each set, we generate 50 mock AGNs.The final uncertainty of each parameter is the quadrature sum of the uncertainty from the fit to observed data and the standard deviation from the fit to mock AGNs.
Due to the bright magnitude of the host (19.2-21.3mag, ∼ 8 mag brighter than the image depth) and moderate AGNto-host ratio (0.7-2.6), the uncertainties in magnitudes are very small ( 0.01 mag).There are also minor systematic biases between results from mock AGNs and input values.Systematic biases in host magnitude are ∼0.05mag for SS2 and SS3 in the F444W filter, and 0.06 mag for DS1, DS2, and DS3 in all filters.Systematic biases in n are slightly larger, with ∼ 9% for SS2 and SS3 in the F444W filter, 5%, 5%, and 0% for n1 and 17%, 13%, and 14% for n2 in DS1, DS2, and DS3, respectively.Almost no systematic bias ( 1%) is found for R e in any set of fits.The low systematic biases for parameters in SS1 are likely due to the large n, which hits the upper limit (7) in the F444W filter.The results from our six sets of multiwavelength simultaneous AGN+host image decomposition are summarized in Table 3. Systematic biases from the mock AGNs are not applied, as they are smaller than the differences due to adopting different PSF models or adopting different host models.For example, the differences of m Sérsic from using different PSF models are on the order of ∼ 0.1 mag, while the differences from adopting different host models are on the order of ∼ 0.3 mag.These latter systematic uncertainties are substantially larger than formal measurement uncertainties reported by the fit and scatter from mock AGNs.An example of decomposition with the single Sérsic model and global PSF model (SS2) is presented in Figure 14.Thanks to the excellent PSF of NIRCam, the total observed light at radius 0. 2 in the SW filters is dominated by resolved emission from the host galaxy.We find prominent features of a central elongated structure, likely a bar or a lens, and galactic spiral arms in the data−nucleus images.The relatively high global Sérsic index n from three SS sets (3.2-4.0) in the F277W filter, which corresponds to rest-frame ∼ 1µm, is primarily due to the combination of a compact high-surfacebrightness bar/lens and an extended low-surface-brightness disk.On the other hand, n of the second Sérsic component in the three DS sets is ∼ 1.0, consistent with the value for exponential disks (n = 1).Our results highlight the danger of morphological classifications based solely on n from single Sérsic model fits with the presence of resolved substructures in the host galaxy, in particular high surface brightness compact structures such as a bar/lens.

Stellar Mass Estimation
We correct the foreground galactic extinction adopting the extinction curve from Cardelli et al. (1989) with R V = 3.1 and the extinction map from Schlegel et al. (1998).As mentioned in Bagley et al. (2022), the relative and absolute photometric calibration accuracy of the NIRCam images in their data product were still at the level of ∼ 5%.We assume a conservative 10% uncertainty for host fluxes.We then derive the stellar mass by feeding AGN contamination-subtracted host galaxy fluxes to the spectral energy distribution (SED) fitting code CIGALE v2022.1 (Boquien et al. 2019;Yang et al. 2022).The SED fitting is performed on all decomposition model sets (Table 3).
For the models of CIGALE, we adopt a single stellar population model from Bruzual & Charlot (2003) with solar metallicity, a Chabrier (2003) IMF, and a "delayed" star formation history.We also include a nebular emission line model with ionization parameter U = 10 −2 , and a dustatt modified starburst model with extinction curve adapted from Calzetti et al. (2000).We adopt a wide range of stellar ages: 0.50-3.85Gyr with a step of 0.05 Gyr; an e-folding time of the main stellar population (Gyr): 0.001, 0.05, 0.1, 0.25, 0.5, 0.75, 1.0, 1. 25, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20; and color excess of the nebular lines [E(B − V ) (mag)]: 0, 0.001, 0.005, 0.01, 0.03, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7.All other parameters are fixed with their default values.Since the longest filter F444W only covers up to rest-frame 1.7 µm, we do not consider dust emission models.An example of SED fitting to host photometries from DS2 is presented in Figure 15.We adopt a Bayesian estimate of the stellar mass, which is based on the probability density distribution of the likelihood of all templates.Uncertainties of the stellar masses are 0.08-0.09dex (∼ 18 − 21%) for all six sets of fits.This level of uncertainty is larger than that of the input flux (10%), suggesting that the uncertainty is dominated by the degeneracy of the templates.However, we caution that the systematic uncertainties due to the assumptions in the fitting, such as the IMF and star formation history, can be significantly larger than the formal stellar mass uncertainties reported by the fitting (e.g., Conroy 2013;Lower et al. 2020).Stellar mass measurements are presented in Table 3.The bright emission, large size and relatively high contribution of the host galaxy in SDSS1420+5300A enable us to investigate the effects of PSF mismatch and host model mismatch on the derived host properties.

The Effect of PSF and Model Mismatch
As mentioned in Section 5.1, the FWHM of PSF models is generally in ascending order of SS1-SS2-SS3 (DS1-DS2-DS3).For the three SS sets, we find clear trends of increasing host flux, decreasing AGN flux, and increasing n toward narrower PSF models.The difference among the three SS sets is ∼ 0.1 mag in m Sérsic , 0.02-0.07mag in m AGN , 9-16% in AGN-to-host ratio, 1-9% in R e , and ∼ 25% for n.The overall trends of these parameters are consistent with the mock analysis presented in Section 4. Due to the complication of two Sérsic models, the trends of parameters in the DS sets are not as clear as for the SS sets.The parameter differences are larger for the first Sérsic model, likely a result from the larger AGN-to-host ratio and smaller host size.A comparison of the best-fit parameters from six sets of decomposition configurations is shown in Figure 16.
Since host galaxies are fitted with two Sérsic models in the DS sets, we could only compare the total AGN magnitude, total host magnitude and stellar mass between SS and DS sets.We find that host magnitudes tend to be larger (lower flux) while AGN magnitudes tend to be smaller (higher flux) in DS models than in SS models.These differences lead to ∼ 15 − 60% higher AGN-to-host ratio from F115W to F444W in DS sets compared to SS sets.The stellar masses are ∼ 0.11 dex (29%) lower in the DS sets compared to the SS sets.This offset is larger than the uncertainty of individual values (18-21%) and the differences among the three PSF models for the same model configuration (0.04 dex for SS and 0.07 dex for DS).χ 2 ν is smaller (better fits) for DS sets.The improvement is expected given the addition of an extra component and more free parameters.Fits from all three PSF models have similar χ 2 ν , with broadened global PSF models performing slightly better.Due to the overall similar χ 2 ν , we cannot determine which PSF model gives the best results.However, it is likely that the broadened global PSF model provides less biased results according to the mock analysis in Section 4.2.
For stellar mass, we find that the effect of model mismatch (0.11 dex on average) is larger than the effect of PSF mismatch (0.04 dex for SS and 0.07 dex for DS).We would expect more significant impact of PSF mismatch in systems with more AGN-dominated light.

Comparison with Earlier Work
The NIRCam imaging data of SDSS1420+5300A were also analyzed by Ding et al. (2022) using galight with a single Sérsic model for the host emission.There are a few key differences in the PSF construction and image decomposition settings.Ding et al. (2022) construct a PSF library by sorting FWHM of bright point sources and selecting isolated and sharp ones in the FoV of four pointings.They then use the top 3, 5, and 8 PSF models with good performance in terms of χ 2 ν of the fit and construct a combined average PSF model by stacking them.The final reported measurements are performed using either individual PSF models or Comparison of MBH versus M * relation of SDSS1420+5300A (cyan star) with local relations from inactive galaxies (Greene et al. 2020).Red and blue dashed lines represent local relations derived from early-type and late-type galaxies, respectively, with shaded region indicating the intrinsic scatter.Open circles represent X-ray selected 1 < z < 2 broad-line AGNs from Ding et al. (2020).We recalculate their MBH using Equation 5 from Woo et al. (2015) to match the choice of virial factor adopted for SDSS1420+5300A in Grier et al. (2019).Filled triangles represent SDSS-RM quasars at 0.2 z 0.8 with MBH from reverberation mapping in Li et al. (2023).Typical uncertainties are indicated at the lower-right corner.
combined PSF models based on the smallest χ 2 ν value, with uncertainties derived from the dispersion of results using the top-5 PSF models.They perform image decomposition independently for each filter, so that they can use a finer pixel scale of 0. 015 in SW fitlers.
Since the photometric calibration reference file in Ding et al. (2022) is different with the one from Bagley et al. (2022), we only compare their AGN-to-host ratio, R e , and n with those from our SS configurations.Their R e is ∼ 20% smaller than ours in four NIRCam filters listed in their paper (F150W, F200W, F356W, and F444W).Their n is approximately constant around 2 in all four filters.This is inconsistent with the n values from our SS sets, which show a clear positive wavelength dependence (e.g., 2.7 to 5.8 for SS2), as expected since older stellar populations have more concentrated distributions.We also find that their AGN-to-host ratios are systematic higher than ours in all four filters.Despite the fainter host magnitudes, their derived stellar mass is on average ∼ 0.12 and 0.25 dex higher than that from our SS and DS sets, respectively.Since we adopt the same cosmological parameter and IMF, it is likely that at least part of the systematic difference is due to the different SED fitting codes used.5.6.SDSS1420+5300A on the M BH − M * Plane Previous work on the M BH − M * relation of z ≈ 2 AGNs primarily relied on HST imaging for AGN-host decomposition (e.g., Jahnke et al. 2009;Bennert et al. 2011;Schramm & Silverman 2013;Ding et al. 2020), with a few studies relying on SED decomposition (e.g., Bongiorno et al. 2014;Suh et al. 2020).The lack of rest-frame NIR coverage (> 0.5µm) and ∼ 1.6 kpc resolution of HST make the decomposition and stellar mass measurements less robust.Also, the M BH estimates in previous studies are mostly based on the singleepoch virial method, which may suffer from complicated systematics resulting from factors such as the poorly understood size-luminosity of the broad-line region at high-z, statistical biases due to sample selection, and different lines adopted in the BH mass estimation (e.g., Mg II and C IV) other than Hβ (e.g., Shen 2013;Shen et al. 2015).
With a more robust M BH estimate from the Sloan Digital Sky Survey Reverberation Mapping project (Shen et al. 2015;Grier et al. 2019), we can now place SDSS1420+5300A on the M BH − M * plane.We adopt the average M * estimate from the DS sets, M * = 10 11.1±0.1 M , given their better modeling of the host structure.Comparing with the local M BH − M * relation from inactive galaxies in Greene et al. (2020), we find that this object lies ∼ 0.6 dex above the upper envelope (best-fit relation + 1× intrinsic scatter) of local late-type (spiral) galaxies, but consistent with early-type (elliptical and S0) galaxies (Figure 17).The position of this object is in broad agreement with previous studies with host fluxes derived from image decomposition using X-ray selected broad-line AGNs at similar redshifts (e.g., Jahnke et al. 2009;Bennert et al. 2011;Schramm & Silverman 2013;Ding et al. 2020) and SDSS-RM quasars at 0.2 z 0.8 (Li et al. 2023) with M BH from reverberation mapping.An increase of at least 0.4 dex in M * is required to reconcile this object with the local relation for late-type galaxies.This amount of growth is feasible as local Sab-Scd galaxies are found to double their mass within 8 Gyr (Bellstedt et al. 2020).Alternatively, gas-rich major merger, which is more common at high redshift (e.g., Lotz et al. 2011), could change the morphology of the object and bring it back to the local relation.A more comprehensive analysis of the M BH − M * relation at z ≈ 2 requires a larger sample with well understood selection functions, reliable BH masses from reverberation mapping and stellar masses from JWST imaging.

SUMMARY
In this work, we characterize the properties of the NIR-Cam PSF and their spatial variations in eight filters (F070W, F115W, F150W, F200W, F277W, F356W, F444W, and F480M), using public NIRCam imaging data.Using simulated AGN+host images, we investigate the recoverability of host properties with imaging decomposition, and the impact of mismatched PSF models.As an illustration for our methodology, we applied the optimal PSF construction method to NIRCam images in the CEERS field, and measured host properties for the broad-line AGN SDSS1420+5300A at z = 1.646.Our main results are as follows: • By comparing three commonly-used methods for PSF construction, we find that PSF models constructed with PSFEx have the the best quality.PSFEx surpasses the other two methods by producing PSF models with the highest SNR and the best performance in modeling the light profiles of point sources in the image.
• We find clear PSF spatial variations across the FoV of NIRCam, with maximum and RMS variations of the FWHM decreasing from ∼ 20% and 5% in the F070W filter to ∼ 3% and 0.6% in the F444W filter.These spatial variations of NIRCam PSFs can have significant consequences on the imaging decomposition of broad-line AGN+host systems.
• We recommend the use of PSFEx to construct the PSF model, as well as modeling its spatial variations.
If a sufficient number of point sources are available, it is recommended to construct region or local PSF models.Otherwise, adopting the global PSF model still provides satisfactory results under most circumstances, with an average 0.02 mag systematic offset and 0.05 mag random scatter.
• We generate simulated images of mock AGNs spanning a wide range of parameters to investigate the reliability and biases of parameter recovery from AGNhost image decomposition using perfect, narrower and broader PSF models.The narrower or broader PSF models are representative of the maximum spatial variations (e.g., 3-20% fractional changes in the FWHM) of the PSF found in this study.PSF mismatch has a profound impact on the recovered host properties, leading to larger scatter, greater systematic biases, and fewer objects with successful measurements compared with the results from using a perfect PSF model (which is unavailable in realistic situations).
Fluxes of host galaxies are more likely overestimated using either a broader or a narrower PSF model.The broader PSF model tends to produce flatter host galaxies (smaller n), while the narrower PSF model tends to produce more concentrated and smaller host galaxies (larger n and smaller R e ).Narrower PSF models on average have a much worse impact on the recovery of host properties than broader PSF models.The systematic biases in the recovered host parameters generally increase as the AGN-to-host ratio increases.In cases where the AGN light dominates, the host flux can be overestimated by as much as a factor of few using PSF model narrower than the actual PSF at the location of the target.
• Given the higher measurement success rate and lower biases for recovered host parameters with broader PSF models than with narrower PSF models, a practical suggestion is to slightly broaden the PSF model by the typical level of spatial variations across the FoV when very few point sources are available for accurate modeling of the PSF spatial variation.This strategy can mitigate the adverse effects and reduce significant biases from inadvertently adopting a narrower PSF model in the AGN+host decomposition.
• The deviations of recovered parameters from the input values are larger than the formal measurement uncertainties for the vast majority of cases where successful measurements are obtained (measurement/uncertainty> 3).For example, median |deviation|/uncertainty of m Sérsic increases from ∼ 1.2 for AGN-to-host ratio=0.1 to 5.0 for AGN-to-host ra-tio=10, with extreme values as high as ∼ 100.It is strongly recommended to generate object-tailored mock AGNs with input parameters from the best fit to evaluate systematic biases and to obtain more realistic measurement uncertainties.
• When the centers of the AGN and its host galaxy are not tied together, one can easily measure an artificial offset between them, dominated by the offset of the fitted host center.This artificial offset scales with the host surface brightness and is approximately aligned with the position angle of the host galaxy.The artificial offsets can be as large as ∼ 80%, 26%, and 7% of R e in objects with R e =0. 12, 0. 48, and 1. 92 when approaching the SB detection limit, respectively.
• We find that the effects of PSF mismatch on the recovered properties of real AGN hosts (e.g., SDSS1420+5300A) are similar to those found in the mock data.The presence of small-scale high-surfacebrightness substructures (e.g., bar and lens) would lead to high Sérsic indices n from the PSF+single Sérsic fitting.Morphological classifications solely based on n can lead to incorrect results.The impact of model mismatch (i.e., using incorrect model components for the host galaxy) may be greater than that of PSF mismatch in systems where the light contribution from the host galaxy is comparable to or larger than that from the AGN, particularly in the presence of substructures.

Figure 1 .
Figure 1.Final F200W mosaic image with pixel scale of 0. 03 pixel −1 .Labels, vertical and horizontal dashed lines indicate the 36 regions for the characterization of the spatial variation of the PSF.

Figure 2 .
Figure 2. The region PSF model images in different regions of the F070W filter constructed using Swarp (top-left), photutils (top-right), and PSFEx (bottom).The number of point sources used to construct the PSF model, axis ratio (minor-axis/major-axis), and FWHM along the major-axis are shown in the lower-left corner of each panel.

Figure 3 .
Figure 3. Normalized radial surface brightness profiles of PSF models in the F070W filter at different regions generated using (a) Swarp, (b) photutils, and (c) PSFEx.Half-width-at-half-maximum of the global PSF is indicated with vertical black dashed lines.

Figure 4 .
Figure 4. Simultaneous AGN+host decomposition in six-band NIRCam images of the broad-line quasar SDSS1420+5300A at z = 1.646.Images of data, model (AGN+host), data−nucleus, and residual (data−model) are shown from left to right.The best-fit parameters for the host galaxy (assuming a single Sersic profile for simplicity) are shown at the lower-left corner in the model panel.The right-most column shows the radial surface brightness profiles of data and models.The spiral arm features in the residual panel are faint and not included in the model.

Figure 5 .
Figure 5. Normalized probability density distributions of reduced χ 2 (χ 2 ν ) of fitting to point sources using PSF models constructed from different methods, estimated from Gaussian kernels using scipy.stats.gaussiankde.Maroon, red, orange, gold, green, lime, and blue lines represent results from SWarp global, SWarp region, photutils global, photutils region, PSFEx global, PSFEx global, and PSFEx local, respectively.PSFEx local represents the PSF models constructed at the exact location of the source in the image.Columns (a-d) indicate results for (a) 10 brightest sources, sources with magnitude between (b) 0th and 5th percentiles, (c) 5th and 16th percentiles, and (d) 16th to 100th percentiles, respectively.Rows from top to bottom present results in F070W to F444W filters.The median and upper and lower 1σ (84th and 16th percentiles − median) of distributions of χ 2 ν are shown at the upper-right corner.

Figure 6 .
Figure 6.Magnitude differences (∆ = PSF model − PSFEx global) of point sources measured using different PSF construction methods versus magnitudes (mPSF) measured using PSFEx global.Colors are the same as Figure 5. Median differences and standard deviations of ∆mPSF are shown at the top-left corner.Uncertainties of mPSF are 0.01 mag and are thus not shown.Due to < 5 point sources available in any region of the F480M image, we only show comparison for global PSF models in that filter.Combined with Figure 7, our results indicate that (1) the quality of the PSF models from SWarp and photutils can be significantly affected by background; (2) PSF variation does not significantly affect flux measurements for point sources.See Section 3.2.2 for details.ral variation of the PSF.Moreover, since most deep fields by design targeting at high galactic latitude, where few bright stars are available for PSF model construction, the effect of the number of usable point sources on the quality and properties of the PSF also merits further investigation.Last but not least, NIRCam employs Hawaii-2RG HgCdTe detectors, which are found to have a brighter-fatter effect (e.g.,Plazas et al. 2018).A comprehensive and quantitative analysis of this effect on the FWHM and profile of observed NIRCam PSFs would be important for studies in crowded fields such as globular clusters.

Figure 7 .
Figure 7. Same as Figure 6 but with the minuend in the ∆mPSF substituted by mPSF obtained from PSF models with half its original size, to reduce the contamination from imperfect background subtraction.

Figure 8 .Figure 9 .
Figure8.Differences (∆ = output − input) between the best-fit results and the input values for m Sérsic (first row), mAGN (second row), Re (third row), and Sérsic index n (fourth row) versus m Sérsic using the perfect PSF model.The left column shows mock AGNs with fixed Re and q but different n, while the right column shows mock AGNs with fixed n and q but different Re.Red, orange, and blue colors represent objects with AGN-to-host ratio = 0.1, 1, and 10, respectively.Open symbols represent objects with measurement < 3× its reported uncertainty.The vertical cyan dashed line indicates the 5σ point source depth.Error bars represent uncertainties for individual objects as reported by galfit.

Figure 10 .
Figure 10.Same as Figure 8 but for the results from the narrower PSF model.

Figure 11 .
Figure 11.Same as Figure 8 but for differences of properties versus AGN-to-host ratio in mock AGNs with m Sérsic = 24 mag.Symbols are the same as Figure 8 with red, orange, and blue colors representing results from the fiducial, broader, and narrower PSF models, respectively.

Figure 13 .
Figure 13.Offsets of host center (∆d [Sérsic]; first row) and AGN center (∆d [AGN]; middle row) between measurements in mock AGNs and input values, and offsets between center of AGN and host measured in mock AGNs (∆d [Sérsic−AGN]; bottom row) versus versus mean surface brightness of the Sérsic model (µ S érsic ).Mock AGNs are fitted with the fiducial PSF model but with centers of PSF and Sérsic model untied.Color and symbols are the same as Figure 8, but with open symbols indicating objects with offset larger than its uncertainty.

Figure 14 .
Figure 14.Multiwavelength simultaneous AGN+host decomposition in six-band NIRCam images of the broad-line quasar SDSS1420+5300A at z = 1.646 with the single Sérsic component and the global PSF model (SS2) using galfitm.Images of data, model (AGN+host), data−nucleus, and residual (data−model) are shown from left to right.The best-fit parameters for the host galaxy and host-to-total flux fraction f host are shown at the lower-left corner in the model panel.The right-most column shows the radial surface brightness (µ) profiles of data and models, with the bottom panel showing the difference between data and model (∆µ = µ data − µ model ).The apparent spiral arm features in the residual panels indicate the presence of a stellar disk.

Figure 15 .
Figure 15.SED fitting to the AGN-subtracted host photometry in the six NRICam filters from decomposition configuration DS2 using CIGALE.The observed fluxes are shown in purple open circles.The best-fit total model spectrum is shown in the black solid curve, with its components of stellar attenuated emission and nebular emission shown in the yellow and green solid curves.Red filled circles represent the model fluxes and the dashed blue curve represents the stellar emission before attenuation.Relative residuals (Data−Model)/Data are shown in the lower panel.

Figure 16 .
Figure16.Comparison of best-fit parameters (a) Re, (b) n, and (c) AGN-to-host ratio versus total host magnitude (m host ) from six sets of decomposition configurations.Black, red, orange, cyan, green, and blue symbols represent results from SS1, SS2, SS3, DS1, DS2, and DS3, respectively.Open symbols represent results of Sérsic1, and filled symbols represent results of Sérsic2.Circle, triangle, square, cross, star, and diamond symbols represent results in F115W, F150W, F200W, F277W, F356W, and F444W, respectively.Individual uncertainties are smaller than the size of symbols and thus are not shown.
Figure 17.Comparison of MBH versus M * relation of SDSS1420+5300A (cyan star) with local relations from inactive galaxies(Greene et al. 2020).Red and blue dashed lines represent local relations derived from early-type and late-type galaxies, respectively, with shaded region indicating the intrinsic scatter.Open circles represent X-ray selected 1 < z < 2 broad-line AGNs fromDing et al. (2020).We recalculate their MBH using Equation5fromWoo et al. (2015) to match the choice of virial factor adopted for SDSS1420+5300A inGrier et al. (2019).Filled triangles represent SDSS-RM quasars at 0.2 z 0.8 with MBH from reverberation mapping inLi et al. (2023).Typical uncertainties are indicated at the lower-right corner.

Table 1 .
Properties of PSF Models

Table 2 .
Properties of PSF Models Filter name.Col. (2-3): 5σ point source (PS) depth and 3σ surface brightness (SB) depth, where σ is the standard deviation of 50 flux measurements randomly distributed in the background image within an aperture of radius = PSF FWHM and 0. 5, respectively.Measurement for PS is corrected using value in Col. (4) for fluxes outside the aperture.