HybPSF: Hybrid PSF reconstruction for the observed JWST NIRCam image

The JWST provides us an unprecedented chance to gaze our universe. With its powerful observation capability in near infrared wavelengths, more details about distant objects could be observed and variety properties of the universe at high-redshift can be investigated more precisely.

We are interested in galaxy clusters which are the largest gravitationally bounded structures in the universe. By studying the mass distribution of galaxy clusters, we can learn more about the dark matter clustering properties (Bradač et al., 2008; Kneib & Natarajan, 2011; Shan et al., 2012; Sereno et al., 2018; Limousin et al., 2022), the investigations of star formation efficiency in cluster contexts by comparing the distribution of the ratio of stellar mass to total mass within galaxy clusters(Bahcall et al., 1995; Jauzac et al., 2012; Bahcall & Kulier, 2014; Hoag et al., 2016; Jauzac et al., 2016; Finney et al., 2018). Studying the high redshift universe with the help of magnification effect of clusters(Strait et al., 2018; Welch et al., 2022; Morishita & Stiavelli, 2023; Sharon et al., 2023). And the nature of dark matter(such as the self-interaction cross-section (Markevitch et al., 2004; Randall et al., 2008; Wittman et al., 2018; Harvey et al., 2019; Gilman et al., 2021; Andrade et al., 2022)) could also be investigated from the spatial offset of the baryonic and dark-matter profiles in clusters.

Accurate mass distribution measurements of galaxy clusters require combining both strong and weak lensing signals in mass modeling methods (Bradac et al., 2004; Bradač et al., 2005, 2006; Oguri et al., 2012; Umetsu et al., 2016; Wong et al., 2017; Congdon & Keeton, 2018). In the weak lensing regime, the precision of shear measurements determines the accuracy of mass distributions, and one of the most significant systematic errors comes from inaccurate Point-Spread-Function (PSF) modeling (Schneider, 2005; Schneider et al., 2006; Mandelbaum, 2018). The PSF describes how point sources would appear in the image and smears the intrinsic galaxy image, ultimately affecting the accuracy of galaxy shape measurements. The PSF is mainly affected by atmospheric turbulence (for ground-based telescopes), telescope optics, and detector effects. The accuracy of shear measurements is heavily affected by imprecise PSF modeling, and reliable mass modeling for galaxy clusters also requires credible PSF modeling. Additionally, astrometry and photometry of point sources cannot be done with high precision without PSF-fitting (Anderson & King, 2000; Mighell, 2005; Nardiello et al., 2022). Therefore, accurate PSF modeling is crucial for various astronomical measurements and analyses, including mass distribution measurements of galaxy clusters, shear measurements, astrometry, and photometry.

Generally, PSFs should be first measured at star positions and then properly interpolated to galaxy positions (Mandelbaum, 2018; Jarvis et al., 2021; Zuntz et al., 2018), and a certain number of stars are usually required for PSF modeling. However, for space telescopes like HST and JWST, star images are often undersampled (Adorf, 1989; Lauer, 1999; Wylezalek et al., 2022; Smith et al., 2009), and the number of PSF stars in a single exposure might not be enough to construct the PSF model for classical PSF modeling methods (van der Marel et al., 2002; Rhodes et al., 2007; Gillis et al., 2020). Therefore, forward modeling methods (e.g., Tiny Tim for HST (Krist et al., 2011) and WebbPSF¹¹1https://webbpsf.readthedocs.io/en/latest/intro.html (Perrin et al., 2012, 2014, 2015) for JWST) are essential for obtaining the PSF of observed images.Forward modeling methods can model the PSF by simulating PSFs for a telescope’s instruments in imaging modes. However, there may be discrepancies between the model and data due to differences in instrument configurations between simulations and practical telescope status (Hoffmann & Anderson, 2018). To obtain more accurate PSF models, some improved methods have been developed (Rhodes et al., 2007; van der Wel et al., 2012; Brooks et al., 2019; Liaudat et al., 2021, 2023). For example, WaveDiff (Liaudat et al., 2023) reconstructs the best fit wavefront based on a differentiable optical forward model according to the data to obtain better PSF models. WebbPSF calculates the PSF by using the results of in-flight wavefront sensing measurements (Optical Path Difference files, OPD), which are produced by the JWST Wavefront Sensing Subsystem (WSS) roughly every two days in flight (Brooks et al., 2019). However, there is still a prominent discrepancy at the inner region of the PSF between data and WebbPSF models (Ding et al., 2022; Ono et al., 2022).

In this paper, we introduce a hybrid PSF modeling method, HybPSF, for JWST NIRCam images that combines simulated PSFs from WebbPSF with observed data. We analyze the residuals obtained by subtracting the WebbPSF PSF model from the data using Principal-Component-Analysis (PCA) method and construct the residual model in the images. PCA is a statistical non-parametric procedure that converts a set of observations of possible variables into a set of linearly uncorrelated principal components and is popular for data compression and information extraction (Pearson, 1901; Shlens, 2014; Li et al., 2016; Nie et al., 2021a). The final PSF model is a composition of the WebbPSF model and the constructed residual model. SMACS J0723 is a massive galaxy cluster (z$\sim$0.388) (Ebeling et al., 2001; Repp & Ebeling, 2018) previously observed by the Hubble Space Telescope (HST) in the HST/ACS and WFC3 cameras (Coe et al., 2019). JWST’s First Deep Field has released observed data of the massive galaxy cluster SMACS J0723 (Pontoppidan et al., 2022), which were performed with NIRCAM filters F090W, F150W, F200W, F277W, F356W, and F444W. We will use the SMACS J0723 NIRCam imaging data to demonstrate our method. The advantage of SMACS J0723 is that there are more stars, which are from LMC (Summers et al., 2023), that can be used in PSF modeling and facilitate accurate shear measurement.

This paper is structured as follows: We introduce the data reduction for star image selection and the discrepancy between the model and data in Section 2 and 3 respectively. The hybrid PSF model and its comparisons with WebbPSF model is described in Section 4. We provide a summary of this work in Section 5.

In this work, we use the JWST first released NIRCam Imaging data of SMACS0723, including F090W, F150W, F200W in the short wavelength channel and F277W, F356W, and F444W in the long wavelength channel, to develop and validate our PSF modeling method. We started our data processing based on the ”stage two” products, which are processed by JWST²²2https://jwst-pipeline.readthedocs.io/en/stable/jwst/introduction.html version 1.8.1 in the context of $jwst\_1017.pmap$ (Bushouse et al., 2017, 2019, 2022). The ”1/f” noise (Livingston, 1965; Rest, 2015) is also removed using a method developed by Wang et al. (in prep.)³³3https://github.com/chriswillott/jwst. We mark cosmic rays or bad pixels using an improved statistical algorithm in multiple exposures coaddition to reduce their impact on PSF reconstruction. To minimize the effect of cosmic rays or bad pixels removal on PSF reconstruction, we fill normal flux values to the marked pixels generated by the statistics over the normal coadded pixels on the same position. Finally, we extract our PSF star images from the above processed data and convert pixel values into counts (Nardiello et al., 2022).

We select our PSF star image candidates according to the “FLUX$\_$RADIUS”, “FWHM$\_$IMAGE” and “FLUX$\_$AUTO” parameters measured from SExtractor (Bertin & Arnouts, 1996, 2010). In the parameter space, stars appear as a locus of points with constant size and are separate from the larger cluster of galaxies (Jarvis et al., 2021). Figure 1 shows an example of the selected PSF stars from channel F444W, “$jw02736001001\_02105\_00009\_nrcblong\_cal.fits$” in the parameter space where the star locus appears separate from the other points. We use the Friends-of-Friends (FOF) method to analyze the images and select stamps that have multiple sources as blended stamps (Zhang et al., 2019). Each source is identified as a number of connected pixels that are above 2$\sigma$ by using the FOF method, and the number of pixels in each group should be larger than 7. We found that over $20\%$ and $80\%$ of the stamps in short and long wavelength channels are blended, respectively. Next, we select unblended stamps as PSF star image candidates. For the blended stamps, we choose those stamps with the disconnected pixels of multi objects as candidates. To avoid the effect of bright pixels from non-star objects, the pixels of objects located at the edge of the stamps are masked and set to zero. Meanwhile we construct a mask image $\mathcal{M}$ for each stamp in which the masked pixels have zero value and the others are set to 1 for each selected blended image.

Besides, faint star images generally are inappropriate for PSF modelings (Paulin-Henriksson et al., 2009; Zhang et al., 2022). In order to use as many stars as possible for PSF modeling, we choose stamps that have SNRs (Signal-to-Noise Ratios) larger than 10 as PSF star candidates. Our star stamps are arranged into $80\times 80$ pixels ($\sim 2.48^{\prime\prime}\times 2.48^{\prime\prime}$ for short wavelength channel and $\sim 5.04^{\prime\prime}\times 5.04^{\prime\prime}$ for long wavelength channel), and with the result that $\sim 40\%$ of the star images have SNR $<40$. Despite the low SNR being inadequate for typical PSF modeling, the central region of these images still offers valuable information for HybPSF due to its emphasis on discrepancies between data and WebbPSF. Finally, we further preserve candidates with ellipticity within 3$\sigma$ to exclude the potential contaminated images. Figure 2 shows examples of star stamp candidates and the final selected star images of $jw02736001001\_02101\_00002\_nrcalong\_cal.fits$, which is the second exposure of Module A in long wavelength channel “F277W”. Processed images show that the contamination at the margin of the stamps is masked effectively when compared with the candidate images, and significant diffraction structures are visible in the star images.

Figure 3 shows the number of selected stars on each chip from all of the NIRCam Imaging data of SMACS0723 used in this work. The histogram of star numbers shows two peaks, corresponding to the short and long wavelength filters, respectively. On average, approximately 12 stars per chip are detected from the short wavelength band data, which is significantly lower than the typical number ( 100) required for precise PSF modeling (Paulin-Henriksson et al., 2009, 2008), and we will discuss later on the consequences of this low number on the analysis.

WebbPSF provides an easy-to-use way to generate accurate PSF for JWST NIRCam images. We first check the consistency of the generated PSFs with the observed data by generating PSFs with WebbPSF using the corresponding positions, detectors, and filters. In addition, the telescope Optical Path Difference (OPD) files are used, and the OPD information is loaded through the “$load\_wss\_opd\_by\_date$” method⁴⁴4https://webbpsf.readthedocs.io/en/latest/jwst_measured_opds.html in WebbPSF, where the argument of “DATE-BEG” can be read from the header of the corresponding fits image.

We compared the profiles of WebbPSF PSFs and star images. The center of the WebbPSF model is generated at the center of a pixel, while the center of the star images is randomly sampled within a pixel. To make a fair comparison, we first generated the PSFs of corresponding star images with an ”oversample” factor of 2. Then, we aligned the centers of the oversampled PSFs with the star images using cubic-spline interpolation. The centers of the star images were estimated using the fast shape estimation algorithm (Li et al., 2012). Next, we downsampled the oversampled PSF images to the same resolution as the data. Finally, we masked both the star images and the corresponding PSFs to have the same masked area, as mentioned in section 2, and normalized their fluxes.

Figure 4 shows the residual images obtained by subtracting the corresponding PSFs from the star images, which were selected from $jw02736001001\_02101\_00002\_nrcalong\_cal.fits$. The residuals exhibit similar patterns, with almost all of them having negative values in the central region. This indicates that the profile of the simulated PSF by WebbPSF is sharper than that of the star images. This discrepancy has also been reported by Ding et al. (2022) and Ono et al. (2022). The difference in value of some residuals could be larger than 0.02. Such a remarkable difference between the PSF model and the actual PSFs could introduce systematic errors into the final measurements of astrometry, photometry, and cluster mass reconstruction.

The general way to construct the PSF model for a galaxy is to first reconstruct the PSF using the observed star images, which represent the exact PSFs at the star positions, and then interpolate the PSF to the target positions in the same image, due to the PSF variations in the field. The number of stars required depends on the different interpolation methods (Bergé et al., 2012; Gentile et al., 2013; Schmitz et al., 2020). Usually, the PSF reconstruction is applied on a single Charge-Coupled-Device (CCD) chip due to the potential discontinuities of the PSF field across the chips on single exposures (Dahle et al., 2002; Jee et al., 2009; Jee & Tyson, 2011; Jee et al., 2013; Finner et al., 2021; Nardiello et al., 2022). This limits the number of stars available for PSF modeling. Moreover, the star images have very complex structures, as shown in Figure 2, which increases the complexity (Paulin-Henriksson et al., 2008) of the PSF model. Constructing a PSF model with high complexity requires numerous stars (Paulin-Henriksson et al., 2008, 2009). However, as shown in Figure 3, the number of selected stars is fairly low, especially for the short wavelength band, which may not be enough to construct a high-precision PSF model.

Principal component analysis (PCA), a non-parametric data compression method, is widely used in PSF modeling for HST (Jee et al., 2007; Nakajima et al., 2009; Schrabback et al., 2010; Finner et al., 2017, 2023c, 2023b). However, the star images are too noisy to construct accurate PSF modeling using the PCA method. Expectation-Maximization-PCA (EMPCA) (Bailey, 2012), an improved PCA method, provides an approach to compose Principal-components (PC) with noisy data-sets or with missing values. Li et al. (2016) used EMPCA to reconstruct the PSF with high accuracy, and Maturi (2017) modeled optical images of galaxies with high precision using the EMPCA method. This suggests that EMPCA might also be suitable for our star images with masked pixels and low SNRs. In the next step, we use the long wavelength data of channel ”F277W” to test the PSF reconstruction efficiency of the EMPCA method. To evaluate the reconstruction effectiveness of EMPCA, we assess the size $R$ and ellipticity $e$ (Paulin-Henriksson et al., 2008, 2009) of the PSFs. These parameters are used to estimate the potential errors in the PSF model constructed by EMPCA. $e$ and $R$ are usually defined based on the central second brightness moments $Q_{ij}$ of the PSF profile (Hoekstra et al., 1998; Bartelmann & Schneider, 2001):

To suppress the noise at the outskirts of the stamps for the $e$ and $R$ measurements, we use a circular Gaussian function with a full-width at half-maximum (FWHM) of approximately 9 pixels.

Figure 5 shows the shape residuals between the data and reconstructed PSFs from EMPCA. The discrepancy between the model and the data could be 0.4 or more, which is mainly caused by the low SNR and limited number of star images used for PSF modeling. The results are far from the general reconstruction accuracy ($\sim$0.01) required for cluster mass by weak lensing (Jee et al., 2005; Holhjem et al., 2009). Therefore, the customary PSF modeling methods (such as EMPCA) might be inadequate for this situation. It is worth mentioning that Finner et al. (2023a) obtained a very excellent PSF modeling accuracy based on the coadded image of channel F200W. The primary difference between our results and theirs is that they utilized star images with a relatively high SNR of 80$\sim$1000 to construct the inner region (with a stamp size of approximately $0.62^{\prime\prime}\times 0.62^{\prime\prime}$ in the channel of F200W) of the PSFs.

Although there are some differences between the PSF model from WebbPSF and the observed data, the models still display excellent agreement with the data. Furthermore, we can see the prominent diffraction spike structures (Pontoppidan et al., 2022), which can also be found in the WebbPSF PSF model (Perrin et al., 2012), in the observed data. However, this structure is hard to reconstruct from the PSF star images because star images are fairly faint, and the spikes are totally drowned in noise. Thus, WebbPSF provides us invaluable information about the PSF outer region profile for JWST NIRCam imaging. And the similarity of the residuals gives us a potential possibility to improve the WebbPSF PSF model. We will introduce our supplementary method for improving the WebbPSF model in the next section.

In this paper, we attempt to improve the WebbPSF PSF model by addressing the residuals. As discussed earlier, the key question now is how to parameterize the residual images. Given the similar structures of the residual images, we can use the PCA method to extract the most dominant structures of residuals. In the PCA scenario, the approach can be expressed as follows:

where $I_{i}$ represents the $i$th star image in the field. $WebbPSF_{i}$ is the corresponding WebbPSF PSF model, and $m_{i}$ is provided for WebbPSF to capture the potential variations in the field. $Res_{i}$ is the corresponding residual image, $P_{l}$ represents the $l$th principal component, and $C_{il}$ is the coefficient of $P_{l}$. Finally, the hybrid-PSF(HybPSF) is defined as:

And the number of principal components $l$ used in this work is $l_{max}=10$, which provides the best trade-off between time consumption and reconstruction performance for the entire residual dataset.

In the classical PCA scenario, the star images should be center-aligned, which is usually achieved through interpolation (Jee et al., 2007; Jee & Tyson, 2011). However, interpolation can cause potential problems such as aliasing effects and smoothing. To avoid these issues, we use a modified iSPCA method (Nie et al., 2021a) to construct the principal components of the residual images. iSPCA uses a set of oversampled basis functions to align the center of observed images without modifying the data. The principal components are then composed of these basis functions: $P_{l}=\Sigma_{m=1}^{m_{\rm{max}}}D_{m}B_{m}$, where $B_{m}$ is the basis function and $D_{m}$ is the coefficients of the basis functions. In Nie et al. (2021a), the Moffatlets basis function showed flexibility in capturing complex structures, so we also use the Moffatlets basis function in this work. Because we want to construct the principal components for the residual images rather than the star images, the $\chi^{2}$ equation of iSPCA can be rewritten as:

In equation LABEL:eq:chi, $W_{ik}$ represents the weight of the $ik$th pixel, where the subscript $l$ and $k$ represent the order of basis functions and pixel index, respectively. In our approach, $W_{ik}$ is the inverse of the estimated noise variance, defined as $\frac{1}{\sigma_{ik}^{2}}$ (Nie et al., 2021a). However, in our case, the star images have masked pixels that need to be considered. Therefore, we set ${W_{ik}}=\frac{\mathcal{M}}{\sigma_{ik}^{2}}$, where $\mathcal{M}$ is the mask image mentioned in section 2.

The basis function in equation (LABEL:eq:chi) can be sampled at any resolution theoretically. However, in this work, we consider the time consumption and the oversampling factor of the WebbPSF PSF model. Therefore, we generate the basis functions with an oversampling factor of $``oversample^{\prime\prime}\times 2$. This means that the constructed principal components are also oversampled. We then downsample the principal components to the same resolution as the WebbPSF PSF model to construct the HybPSF.

It is worth noting that although principal components are used to construct the residuals, the basis functions must be aligned with the center of star images. Additionally, the residuals exhibit the most significant discrepancy in the central region of the images, so we can only calculate the $\chi^{2}$ in the inner region of the stamps. For this work, we used only the central $44\times 44$ pixels of each image to calculate $\chi^{2}$, rather than the entire stamp area, which also helped to reduce time consumption. After calculating the coefficients of $m_{i}$, $C_{il}$ and $D_{im}$ from equation LABEL:eq:chi, we reconstructed the PSF of each star image by combining the $WebbPSF_{i}$ and $Res_{i}$ models from iSPCA, respectively. In order to capture the spatial variations of the PSFs, we used bivariate polynomials to fit the coefficients of the principal components and WebbPSF in the field for each CCD chip (Jee et al., 2007). We found that polynomials below second order were unable to fit the residuals field well, so we only adopted chips with PSF star numbers larger than 6 for PCA. Figure 3 shows the distribution of star numbers on the accepted chips on the right side of the pink dashed line, while chips on the left were not used in this work.

Now we summarize the model process of HybPSF as follows :

1. $Extracting\ star\ images$: In this step, we extract the PSF star images and construct the masked areas and weight images and then normalise each images.

2. $Constructing\ residual\ images$: Generating the corresponding WebbPSF, and subtracting the centrally aligned and flux normalised model from the star images to get the residual images.

3. $Calculating\ \chi^{2}\ equation$: Aligning the WebbPSF to the center of each data respectively, and use the data, WebbPSF, and masked images to calculate the $\chi^{2}$ equation.

4. $Constructing\ the\ {\rm HybPSF}\ model$: Applying the polynomial fitting to the coefficients $C_{il}$ and $m_{i}$ of each chip respectively, and combining the WebbPSF and $Resi$ model at target positions from WebbPSF and fitted polynomial respectively. The maximum order of polynomials depends on the number of PSF stars used in the PSF reconstruction, for example: when the number of PSF stars is greater than 21, the maximum order is 5.

Figure 6 displays the residual image examples of the star images after subtracting the HybPSF, which correspond to $jw02736001001\_02101\_00002\_nrcalong\_cal.fits$. We observe that the maximum difference between HybPSF and the data is less than 0.006, which is nearly one order of magnitude smaller than that shown in Figure 4. The pattern of the residual images shows less similarity, indicating that HybPSF has successfully extracted the most prominent structure in the residuals. However, the residuals still exhibit some structures, which may result from the small number of stars used for PCA decomposition. Overall, almost all of the images show considerable improvement. Figure 7 provides an example of ellipticity comparison between the data, WebbPSF, and HybPSF. The primary features of the ellipticity distributions in the field are captured by HybPSF and WebbPSF, and the distribution of whiskers in the residual panel appears random.

Figure 8 shows distributions of the shape parameter residuals (${\delta R^{2}}/{R^{2}}$,$\delta e_{1}$ and $\delta e_{2}$) for the JWST NIRCam imaging data used in this work, respectively. And $\delta e$ and $\delta R^{2}/R^{2}$ values are obtained by subtracting the $e$ and $R$ of the PSF model from the measurement of the star images. The light-colored distributions in Figure 8 represent the results from all of the samples, including images with low SNRs, which are used as a contrast. As shown in Figure 8, the ${\delta R^{2}}/{R^{2}}$ distributions in channels F150W, F200W, F277W, F356W, and F444W, and the $\delta e_{1}$ distributions in F090W appear separated for the two methods. However, all of the residuals are spread over a wide range for both methods. To avoid the noise effect in shape estimations, we use star images with high SNRs (SNR$>$100) as a contrast to obtain more explicit comparisons. The results, shown in Figure 8 with a heavier color, demonstrate that the ${\delta R^{2}}/{R^{2}}$ and $\delta e_{1}$ distributions for the two methods mentioned above are clearly separated. The peaks of the HybPSF distribution (shown in red) are closer to 0 compared to those of the WebbPSF distribution (shown in blue). The results from high SNR data indicate more compact distributions.

Furthermore, we present the statistical results of $\delta e$ and $\delta R^{2}/R^{2}$ for the two methods in Figure 9. As with Figure 8, the results are obtained for two types of datasets (all samples and high SNR samples). The scatter plots from high SNR data show much smaller ranges than those from all samples. The results of $\delta R^{2}/R^{2}$ are clearly separated, suggesting that the size precision of the PSF model is improved significantly, and the ellipticity precision is also improved to some extent. Table 1 shows the mean shape residuals for WebbPSF and HybPSF in each channel from samples with high SNRs, respectively. The results of HybPSF are closer to zero than that of WebbpSF. Specifically, the mean residuals of $\delta R^{2}/R^{2}$ in the F150W, F200W, F277W, F356W, and F444W channels decreased by an order of magnitude (from a few percent to parts per thousand). Additionally, the values of the ellipticity component $\delta e_{1}$ in channel F090W also show a decline. The 1$\sigma$ uncertainty of the residuals is constrained to few percent, indicating that the model may be acceptable for cluster mass analysis using weak lensing (Jee et al., 2005; Holhjem et al., 2009).

Based on Figures 8 and 9, the mean residual values of HybPSF from all samples and high SNRs samples are statistically consistent. This suggests that low SNRs images could also provide effective constraints for JWST NIRCam imaging PSF reconstruction in our method.

The HybPSF model can be regarded as the central profile calibrated WebbPSF, and the ellipticity has more dependence on the peripheral profile of the PSF according to its definition. Hence, the improvement for PSF size is more evident than that for ellipticity. We also observe that the error of $\delta e$ is smaller than that of $\delta R^{2}/R^{2}$ in Figure 9. Our method can also reconstruct the PSFs with high precision for the ”stage three” data products of JWST, which may be used more widely. To demonstrate this effect, we use the Drizzle coadded images (Fruchter & Hook, 2002; Wang & Li, 2017) from one chip in the long wavelength channels for PSF reconstruction using HybPSF, since they would expect to have more PSF stars with high SNRs. In the coadded images, we choose star images with SNR$>$40 for PSF modeling and images with SNR$>$100 for contrast. Additionally, we restrict the comparison region of the PSFs to the inner $30\times 30$ pixels to avoid noise effects. Generally, the stacked image would not be the average of multiple exposures due to rotation, geometric distortion, etc. (Fruchter & Hook, 2002). For simplicity, we construct the WebbPSF for coadded images by averaging the WebbPSF corresponding to multiple exposures obtained for each star.

Figure 10 displays the shape residual results of the coadded images, and the width of the residual distributions is much smaller than those in Figure 9. The 1$\sigma$ uncertainty of these points is markedly reduced to $\sim$0.01 for $\delta R^{2}/R^{2}$ and parts per thousand for $\delta e$, as shown in Table 2.

Based on the residuals and shape parameter comparisons, HybPSF appears to have better performance than WebbPSF. HybPSF constructs a residual model based on PCA and provides supplements to the WebbPSF PSF model. We expect this could help us to make better measurements for cluster mass, astrometry, photometry, etc. Furthermore, HybPSF can provide a PSF model with a considerable area that is difficult to achieve using only the empirical method based on star images. This advantage is promising for the upsampling deconvolution method (Wang et al., 2022). However, we also observe some structures left in the residuals, as seen in Figure 4. Therefore, the PSF is not reconstructed perfectly by our method, which may be caused by the small number of star images in the single exposure. Additionally, we did not consider the spectral energy distribution (SED) dependence of each stellar object in WebbPSF, which may help to calculate more accurate PSF models.

The shape comparisons above mainly focus on the inner region of PSFs to avoid noise effects. One of the features of HybPSF is that it includes outer profiles, such as the diffraction spikes. We cut a large stamp ($6.3^{\prime\prime}\times 6.3^{\prime\prime}$) with prominent diffraction spikes from the long wavelength channel F356W to test the performance of HybPSF. In this test, we set the coefficient $m$ of WebbPSF to be 1. This is because that the coefficient $m_{i}$ is estimated from faint star images which show no distinct spikes structures in the very outer regions, and the constraints for considerably bright stars might be inappropriate. Figure 11 shows the spike correction of WebbPSF and HybPSF, respectively. The star flux is calibrated using the flux of the inner circular region with a radius of 10 pixels from the center of images. The outer region of HybPSF is constructed from WebbPSF, so they have very similar performance, as shown in Figure 11. This suggests that HybPSF can also capture the spike structures at very large scales of the JWST NIRCam image PSFs. We then estimate the AB magnitude of the star based on the two PSF models, and the difference in AB magnitude estimated from the two PSF models is $\sim 0.01$ magnitude.

With the powerful observation capability of JWST, we are entering a new era of exploring our universe. However, accurate modeling of its PSF is also required for high-precision astrophysical measurements. In this paper, we introduce a PSF modeling method called HybPSF for the NIRCam images, which combines observed data and simulated PSFs from WebbPSF.

To extract proper star images from the observed data, we began with the ”stage two” calibrated file and then filled up the contaminated pixels and performed deblending for each star image in later steps. However, comparing the observed data and the simulated PSF from WebbPSF, we found that the residuals still had significant structures, suggesting that the PSF model could be improved potentially.

We then attempted to extract the remaining signals from the residuals using the iSPCA method, and supplemented the constructed residual model back to WebbPSF to obtain a PSF model (HybPSF) that would be more consistent with the observed data. The value range of the corresponding residuals from HybPSF reveals the improvements of our method. Finally, we checked the model precision on the PSF shape parameters: $e$ and $R$. The $e$ and $R$ of the HybPSF models showed better precision than WebbPSF to some extent, especially on $R$. This suggests that HybPSF is more consistent with the data than WebbPSF. Additionally, our method can provide oversampled PSF images, which could be promising for up-sampling PSF deconvolution (Wang et al. in prep). Additionally, this method is also applicable to i2d.fits images. Given our scientific interest in accurately measuring shear, we aimed to avoid uncertainties in the PSF resulting from image stacking. Therefore, we required accurate PSF information from individual images (cal.fits). The PSF obtained in this work will also be utilized for PSF deconvolution and image stacking (Wang in prep).

However, there are still some discrepancies between our model and the data, which can be seen in the residual images and shape residuals. This may be caused by the limited number of stars used for PSF modeling, and some chips were dropped due to the limited number of stars. This issue may be alleviated by considering all star images in the entire field, as done in Liaudat et al. (2020); Finner et al. (2023a), and galaxy images may also provide some help in solving this problem (Nie et al., 2021b). Additionally, SED information of the stellar objects can help WebbPSF construct more accurate PSF models, which we will explore in future work. The brighter-fatter effect is another important systematic in PSF modeling. However, Finner et al. (2023a) mentioned that this effect would not significantly impact cluster mass reconstruction. We will check for these influences on PSF modeling in the future. Overall, we expect that this improved PSF model will help us obtain better astrophysical measurements, such as cluster mass, photometry, and more.

This work is supported by National Key R$\&$D Program of China No. 2022YFF0503403 and the Ministry of Science and Technology of China (grant Nos. 2020SKA0110100). HYS acknowledges the support from NSFC of China under grant 11973070, Key Research Program of Frontier Sciences, CAS, Grant No. ZDBS-LY-7013 and Program of Shanghai Academic/Technology Research Leader. We acknowledge the support from the science research grants from the China Manned Space Project with NO. CMS-CSST-2021-A01, CMS-CSST-2021-A04.

The Early Release Observations and associated materials were developed, executed, and compiled by the ERO production team: Hannah Braun, Claire Blome, Matthew Brown, Margaret Carruthers, Dan Coe, Joseph DePasquale, Nestor Espinoza, Macarena Garcia Marin, Karl Gordon, Alaina Henry, Leah Hustak, Andi James, Ann Jenkins, Anton Koekemoer, Stephanie LaMassa, David Law, Alexandra Lockwood, Amaya Moro-Martin, Susan Mullally, Alyssa Pagan, Dani Player, Klaus Pontoppidan, Charles Proffitt, Christine Pulliam, Leah Ramsay, Swara Ravindranath, Neill Reid, Massimo Robberto, Elena Sabbi, Leonardo Ubeda. The EROs were also made possible by the foundational efforts and support from the JWST instruments, STScI planning and scheduling, and Data Management teams.