A Hierarchical PSF Reconstruction Method
Pedro Alonso Jun Zhang Cong Liu



Jun Zhang
betajzhang@sjtu.edu.cn


Department of Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China
Shanghai Key Laboratory for Particle Physics and Cosmology, Shanghai 200240, China


Department of Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China
Shanghai Key Laboratory for Particle Physics and Cosmology, Shanghai 200240, China


Department of Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China
Shanghai Key Laboratory for Particle Physics and Cosmology, Shanghai 200240, China


Reconstruction of the point spread function (PSF) plays an important role in many areas of astronomy, including photometry, astrometry, galaxy morphology, and shear measurement. The atmospheric and instrumental effects are the two main contributors to the PSF, both of which may exhibit complex spatial features. Current PSF reconstruction schemes typically rely on individual exposures, and its ability of reproducing the complicated features of the PSF distribution is therefore limited by the number of stars. Interestingly, in conventional methods, after stacking the model residuals of the PSF ellipticities and (relative) sizes from a large number of exposures, one can often observe some stable and nontrivial spatial patterns on the entire focal plane, which could be quite detrimental to, e.g., weak lensing measurements. These PSF residual patterns are caused by instrumental effects as they consistently appear in different exposures. Taking this as an advantage, we propose a multi-layer PSF reconstruction method to remove such PSF residuals, the second and third layers of which make use of all available exposures together. We test our method on the i-band data of the second release of Hyper Suprime-Cam. Our method successfully eliminates most of the PSF residuals. Using the Fourier_Quad shear measurement method, we further test the performance of the resulting PSF fields on shear recovery using the field distortion effect. The PSF residuals have strong correlations with the shear residuals, and our new multi-layer PSF reconstruction method can remove most of such systematic errors related to PSF, leading to much smaller shear biases.






§ INTRODUCTION
 
The point spread function (PSF) measures the diffraction of light in optical systems, which makes point-like sources appear extended. In astronomical images, the PSF effect distorts the shape and size of all celestial objects. For ground-based telescopes, the properties of the PSF are mainly determined by the telescope optics <cit.> as well as the atmospheric turbulence <cit.>. A minor contributor to the PSF is the pixelation effect, resulting from the finite pixel size of the CCD images <cit.>.

PSF plays a particularly important role in weak lensing studies, since it convolves the lensed galaxies and distorts their shapes akin to the cosmic shear. This makes the PSF effect the most important source of systematics in shear measurements <cit.>. To ensure the reliability of the measurements, precise modeling of the PSF is essential.

When modeling the PSF, it is convenient to describe it as a combination of time-variant and time-invariant features. Time-variant features vary from exposure to exposure and are produced by the atmospheric turbulence and some instrument-related issues such as misalignments in the optical components of the telescope, mechanical deformations caused by gravitational or thermal effects, tracking errors, or instrumental instabilities. On the other hand, time-invariant features are quite stable from exposure to exposure.

Optical designs, mechanical elements of the telescope, as well as characteristics of the imaging sensor, focal length, or aperture size can all induce time-invariant features on the PSF.

Multiple methods have been proposed to model the PSF effect for shear measurement <cit.>. With time-variant features complicating the PSF modeling, most of these methods rely solely on individual exposures to build the PSF model, thus being limited by the number of available stars. An effective way to deal with a limited number of stars is to assume a specific functional form for the distribution of the PSF, e.g., polynomial functions. Several methods that follow this approach have been successful in capturing most of the PSF features on the exposure level. Other methods, based on Principal Component Analysis (PCA) <cit.> or machine learning <cit.>, have also shown impressive results in capturing the spatial features of the PSF. The optimal PSF reconstruction method generally depends on the particular shear measurement being considered, and the particular dataset.

In current methods, however, when stacking the PSF ellipticity residuals of the PSF models of many exposures, some complex features appear <cit.>. These features are time-invariant, and therefore related to instrumentation. In this work, we present a novel PSF reconstruction method that performs three interpolations in a hierarchical manner, with the aim of modelling the high-frequency yet stable spatial features of the PSF distribution. Our method builds upon the PSF reconstruction method of Liu et al. 2024 (in preparation), constructing a model of the PSF power spectrum. We apply our method to the i-band of HSC pDR2 data, in which systematic PSF residuals have been previously reported. Our method successfully removes most of the systematic PSF residuals found in HSC DR2, significantly improving the PSF model.

The structure of this paper is as follows. <ref> introduces the HSC dataset. The first layer of interpolation is explained in <ref>. The second and third layers of interpolation are described in <ref>, which also includes our main results. We conclude in <ref>.





§ HSC DATASET


The Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) project <cit.>
is an optical multi-layer imaging survey that covers approximately 1400 deg^2 in five bands (g, r, i, z, y) in its Wide layer (r∼26). In addition, the survey includes two deeper layers, Deep and UltraDeep, covering 27 deg^2 (r∼27) and 3.5 deg^2 (r∼28), respectively. The project utilizes the Hyper Suprime-Cam, a wide field optical camera built on the 8.2m Subaru Telescope, and aims to address some of the most important problems in astrophysics and comsmology, with a focus on weak gravitational lensing, galaxy evolution, supernovae, and galactic structure. Their data is publicly available at their official website[https://hsc-release.mtk.nao.ac.jp/].

In this paper we use the i-band data of the second public data release of the HSC <cit.>. HSC pDR2 covers an area of 300 square degrees in the Wide layer, in all five bands. The data was collected over 174 nights of observation, from March of 2014 to January of 2018. HSC pDR2 includes significant improvements over the previous release (HSC pDR1), which include improved background subtraction, PSF modeling, and object detection procedures. However, we only use the background-removed CCD images, performing our own PSF reconstruction. The galaxies are selected from the official HSC catalog <cit.>. Fig. <ref> shows the distribution of the PSF full width at half maximum (FWHM), in real space, for the i band of HSC pDR2.



§ PSF RECONSTRUCTION





Our method builds upon the PSF reconstruction method of Liu et al. 2024 (in preparation), which is part of the Fourier_Quad shear measurement pipeline <cit.>. Fourier_Quad utilizes the quadrupole moments of the power spectrum of galaxy images to measure the cosmic shear. In line with this approach, we perform PSF interpolation directly on the power spectrum of the stars, building a model of the power spectrum of the PSF (PPS hereafter). 



 §.§ Star selection

Several techniques have been developed for star-galaxy separation. Some methods rely on morphological features <cit.> or color information of the sources <cit.>. Moreover, modern approaches based on machine learning methods such as decision trees <cit.> or deep convolutional neural networks <cit.> have shown promising results in star-galaxy separation.

In this work, we follow the following procedures to select out bright stars for PSF reconstruction (Liu et al. 2024, in preparation):

1. For each exposure, we select sources with SNR≥100 to form our initial set of star candidates. The rest of the steps only make use of their power spectra, which are all normalized so that the power P(k⃗=0) is unity. Note that in Fourier space, point sources should have the most extended profiles. We therefore measure the area of each candidate (defined as the number of pixels above 0.02). The distribution of the area typically has a Gaussian shape. We throw away those candidates that are away from the peak of the distribution by more than 3σ to remove the outliers. 

2. From the remaining star candidates, we form the first star model (power spectrum) pixel-by-pixel.
Each pixel value is determined by sorting the corresponding pixel values from the power of all candidates, and taking the lower bound of the top 25%.

3. The similarity between the candidates and the star model can be quantified by defining a χ^2 as:

    χ^2=2∑_i=1^N(I_i^n-I_i^model)^2/∑_i=1^N(I_i^n+I_i^model)

in which I_i^n refers to the value of the i^th pixel in the power image of the n^th candidate, and I_i^model refers to the corresponding pixel value of the model. N is the total number of pixels involved, which are typically chosen to be those located in the middle part of the stamp. The distribution of the χ^2 forms a Gaussian-like function, and we remove those candidates whose χ^2 are more than 3σ away from the peak of the distribution. 

4. We build the new star model, this time as a function of location, by interpolating the pixel values of all the remaining candidates with polynomial functions of order nine. Note that if there are not enough candidates left for the fitting, we simply stop processing the exposure further. We again use eq.(<ref>) to define the χ^2 between each candidate and the star model at its location. From the distribution of χ^2, we again remove the candidates that are more than 3σ away from the peak. The surviving candidates are treated as stars for our PSF reconstruction. 



 §.§ First layer of interpolation


Our first interpolation is a 2D polynomial fitting of third order on the power spectrum of the stars on each CCD image, following the procedure of Liu et al. 2024 (in preparation). To ensure the reliability of the PSF model, we only include CCDs containing a minimum of 20 stars. Each star power spectrum is centered on a 48x48 stamp, and the interpolation is done pixel-by-pixel.

To evaluate the quality of the interpolation, we calculate the PPS residuals, represented as the ellipticity and relative size residuals at the positions of the stars. The ellipticity components, e_1 and e_2, and the size are defined based on the quadrupole moments, Q_ij, as:


    e_1 =   Q_20-Q_02/Q_20+Q_02
    
    
        e_2 =   2Q_11/Q_20+Q_02
    
    
        size =    Q_20+Q_02/Q_00,

where

    Q_ij = ∑_P(k)>P_0 P(k)k^i_x k^j_y.

We only include pixels above P_0=0.02· P(k=0) in the calculation.

Ellipticity residuals are calculated as the difference between the ellipticities of the original (e_1,true, e_2,true) and the predicted (e_1,pred, e_2,pred) stamps, while the relative size residuals are calculated as (size_true-size_pred)/size_true. Fig. <ref> shows the stacked ellipticity and relative size residuals including all exposures. While the residuals are very small for most part of the exposure, we find some remaining systematic residuals near the boundaries, which are not visible on single exposures but become prominent after stacking a large number of them. Note that similar PSF residual patterns are also observed with the official HSC pipeline <cit.>, as well as in the DES data <cit.> and the DECaLS data <cit.>. In the next layer of interpolation, we build a model for these systematic PSF residuals.



§ IMPROVED PSF RECONSTRUCTION










The systematic PPS residuals observed in fig. <ref> are related to instrumental effects rather than the atmospheric turbulence. Given the strong spatial and temporal dependencies of the atmospheric turbulence, it is unlikely to produce any persistent features on the PSF residuals.
In the following we present our second layer of interpolation, which builds a model for the systematic features of the PPS residuals. Unlike in the first interpolation, the features that we want to model are systematic, common to all exposures. Therefore, instead of performing an interpolation on individual chips, we collect the PPS residual stamps from all exposures and place them into a single exposure. We then perform a single interpolation on each CCD of that exposure. This procedure significantly increases the amount of available data for interpolation.

To collect the PPS residuals from all exposures and place them into a single exposure, we must first take into account the different PSF sizes in different exposures (fig. <ref>). We rescale the PPS residual stamps to a common size, with the re-scaling factor determined as the mean PSF size of the exposure. It is important to emphasize that the above procedure of rescaling the PSF stamps from different exposures and place them into a single exposure is only beneficial when modelling systematic features, common to all exposures. In our first interpolation, the main contributor to the PSF was the atmospheric turbulence, which strongly varies from exposure to exposure, hence the interpolation was performed on individual exposures.
 


 §.§ Rescaling of the PPS residuals





We rescale each PPS residual to a reference PSF FWHM of 1 arcsec in real space. For each PPS residual stamp, we calculate the rescaling factor s as the mean PSF FWHM size, in real space, of all stars in the same CCD in unit of arcsec. To build the rescaled PPS residual stamp, we project the coordinates (u, v) of each pixel back to the original stamp, as (u· s, v· s). In most cases, the projected coordinates (u· s, v· s) lies within four pixels of the original stamp. The value of the (u, v) coordinate in the rescaled stamp is calculated as a weighted average of those four neighboring pixels in the original stamp, with the weights determined by the inverse of the pixel distances to (u· s, v· s).



 §.§ Principal Component Analysis (PCA) of the PPS residuals



Once the PPS residual stamps have been rescaled, we place them into a single exposure, and apply principal component analysis to the stamps on each CCD. Principal component analysis <cit.> is a widely used method for dimensionality reduction in data analysis. It transforms an N-dimensional system into a lower dimensional representation, characterized by Principal Components, or PCs. These PCs form the basis of the new space, representing the main features of the data. They are orthogonal to each other, which helps eliminate the correlation between variables, making PCA particularly useful for large datasets, where numerous correlated variables make data interpretation difficult. The main motivation for applying PCA to the PPS residuals is to capture the main features of the data, eliminating unnecessary information or noise. Additionally, dimensionality reduction drastically decreases the computational time, since instead of performing the interpolation pixel-by-pixel, now we perform an interpolation for each principal component.

We choose 100 as the number of principal components, resulting in each stamp being represented by the 100 coefficients to the PCs. Fig. <ref> and fig. <ref> show the first 30 principal components of the CCDs number 50 and 53, respectively <cit.>. The CCD number 50 is located in the central region of the exposure, whereas the CCD number 53 is located in the boundary. We observe that their first four principal components are almost identical, diverging from the fifth, showing the different features of the PPS residuals in the two regions.



 §.§ Reconstruction of the PPS residuals



We build a model for each principal component coefficient at the CCD level by interpolating the PC coefficients of the PPS residuals within that CCD. We test two different interpolation methods, polynomial and random forest, and compare their performance.



  §.§.§ Polynomial



The features of the PPS residuals are particularly complex near the boundaries of the exposure, thus a single interpolation on a CCD level is unlikely to fully capture them. Since the amount of data per CCD is now very large, we further divide each CCD into four equal parts (2× 2 in the CCD plane), and fit a polynomial of order six to the PC coefficients of the PPS residuals within each part. The predicted PPS residual stamp (48x48) at the position of each star is constructed as the sum of the PCs, weighted by the predicted PC coefficients. Finally, we rescale the current PPS residuals model back to the original PSF size at each position. We follow the same procedure as <ref>, with the scaling factor being the inverse of s.

The improved PPS model is built by adding the PPS residuals model to the original PPS model from the first interpolation. Fig. <ref> shows the new ellipticity and relative size residuals, calculated at the position of the stars, as the difference between the ellipticity components/relative size of the true star and the new PPS model. We observe a significant improvement, with the systematic PPS residuals vanishing almost completely.



  §.§.§ Random Forest


In this section we follow the same procedure as in <ref>, but using a machine learning algorithm called random forest to interpolate the principal component coefficients of the PPS residuals. Random Forest <cit.> is a widely used machine learning algorithm, applicable to both classification and regression problems. It is an ensemble method that combines the predictions of multiple models, called decision trees, to infer a final prediction. Each decision tree is built on a randomly selected subset of the data, and the output of the random forest is the average output of all the decision trees. This significantly reduces the risk of overfitting.

Each decision tree is built through a series of binary splits of the data, starting from the root node. At first, two subnodes of the root node are created, based on a splitting condition. Following the same procedure, each of the subnodes is further divided into two new subnodes, based on new splitting conditions, and so on, building the tree. When a stopping condition is met, we stop splitting that node.

In our case, each decision tree performs a regression on a randomly selected subset of our data, where each datapoint is described by the coordinates (x, y) on the CCD, and the principal component coefficient value, z. Starting from the root node, we perform the first split of the data. The algorithm splits the data into two groups or subnodes (A and B), based on a condition—either on x or y—that better splits the data according to the values of z. It is important to highlight that the algorithm sees x and y as distinct features of the data, rather than as coordinates. As a result, each split is based on a single feature, either x or y. To determine the optimal splitting condition for each node, we use the mean squared error (MSE), a widely adopted metric for splitting in regression tasks. We define the MSE of each possible split as:

    MSE = ∑_(x,y)∈ A(z-z_pred, A)^2+∑_(x,y)∈ B(z-z_pred, B)^2,

where z_pred, A and z_pred, B are the predicted PC coefficients associated to each group, and it is the same for all elements in the group. The split with the lowest MSE value is the optimal split.

We continue splitting each subnode (A and B), following the same procedure as for the root node. We stop splitting a node when it contains less than 500 datapoints. This stopping criterion is found empirically, as further splitting the data results in overfitting.

A decision tree generates predictions for new coordinates (x, y) by following the splitting conditions for each node. Starting at the root node and moving down the tree, the new data reaches a leaf node, i.e., a node without subnodes. The predicted PC coefficient for the new coordinates is the value associated to that leaf node, z_pred.

As in <ref>, we divide each CCD into four parts (2× 2 in the CCD plane), and build models for the PC coefficients of the PPS residual stamps within each part. We build 100 decision trees for each principal component coefficient, and determine the predicted PC coefficient at each CCD position, z_pred, as the average predictions of all decision trees. We build the PPS residual model as the weighted sum of the principal components, with the weights given by the PC coefficients predicted by random forest, z_pred. To obtain the final PPS residual model at each CCD position, we rescale the PPS residual model back to the PSF size of each exposure. Fig. <ref> shows the ellipticity and relative size residuals of our improved PPS model, calculated as the sum of the PPS model from the first interpolation and the PPS residual model. As in the polynomial case, we obtain very small PPS residuals, removing almost completely the systematic PPS residuals of fig. <ref>. Although polynomial and random forest obtain comparable results, the PPS residuals of random forest are slightly smaller. However, this does not necessarily mean that random forest is going to be more accurate in making predictions at the position of galaxies. As a non-parametric algorithm, random forest could in principle capture more complex features on the data. However, they are also more prone to overfitting. To evaluate the predictions at galaxy positions of our new PPS models, we use the field distortion test.



 §.§ Field distortion test


The field distortion <cit.> is an optical aberration characterized by the deviation from global rectilinear projection. It induces a distortion in the shape of galaxies in a similar way as the cosmic shear, which can be directly derived from astrometry parameters. <cit.> proposed a method to use the field distortion to evaluate the accuracy of shear recovery directly on real galaxies, by comparing the measured field distortion induced shear (FDS) and the true field distortion inferred from astrometry. This is known as the field distortion test. We refer to <cit.> for a derivation of the FDS equations.

We use the FD test to evaluate the performance of our PSF model on real galaxies, and compare the results of polynomial and random forest. We evaluate the true FDS, g_1(FD) and g_2(FD), against the FDS recovered from galaxies, g_1(gal) and g_2(gal).

Fig. <ref> compares the FD test results for SNR_F≥ 4 <cit.>, after the first and second interpolation. Note that the FDS signals are removed from the galaxies. The x-axis represents the true FD signal and the y-axis the difference between the recovered FD signal from astrometry and the true one. The results after the first interpolation clearly show the imprints of the PPS residuals of fig. <ref>, and suggest that all galaxies with a FD signal of |g_1,2(FD)|>0.02 should be removed from the shear catalogs, as their PSF models are not reliable. We find a significant improvement after the second interpolation, particularly in the |g_1,2(FD)|>0.02 range. This improvement allows the shear catalog to include a much larger number of galaxies, thus enhancing its statistical power.

Fig. <ref> shows that although random forest obtained slightly smaller PPS residuals compared to the polynomial case, both methods perform equally well in predicting the PPS model at galaxy positions. This results demonstrate the efficacy of the FD test in evaluating the PSF reconstruction accuracy directly on galaxy positions. 



 §.§ Further improvements









After smoothing the PPS residuals shown fig.<ref>, we find an interesting phenomenon: the PPS residuals show a very small but systematic slope, similar to that shown by the field distortion induced shear (FDS). Fig.<ref> displays the FDS (g_1 and g_2), and the ellipticity residuals (e_1 and e_2), which show a clear correlation. This correlation was not observed earlier mainly because this effect is too small and it was covered by the noisy PPS residuals signal.

We perform a third interpolation that aims to correct for this correlation. As we observe in fig.<ref>, the systematic PPS residual is now global, thus we perform a single interpolation on the scale of the whole exposure.

As in the second layer of interpolation, we first rescale the PPS residual stamps—obtained after the second interpolation—to a reference PSF FWHM of 1 arcsec, in real space. We place them into a single exposure and apply principal component analysis (<ref>) to the PPS residual stamps in entire exposure. We use 100 principal components in our analysis. To find this number, we evaluate the quality of our PPS model for different number of PCs, using the field distortion test. Fig. <ref> shows the correlation between the PPS residuals (e_1 and e_2) and the FD signal (g_1 and g_2) for different number of PCs (10, 20, 30, 50, and 100). From 10 to 50 PCs, we observe an clear improvement as the number of PCs increase, with the results stabilizing after 50 PCs, without a noticeable improvement. Based on these results, we conclude that a minimum of 50 PCs is necessary to capture all the features of the PPS residuals. However, it is important to note that these results might vary depending on the specific dataset, thus we recommend maintaining the number of PCs above 50, ideally 100, to ensure optimal performance.

Using all available data, each set of PC coefficients is fit to a polynomial of third order, building a model for the PC coefficient. These models are used to predict the PC coefficients at different positions on the exposure. Next, a model for the PPS residuals is built for each coordinate on the exposure as a weight sum of the PCs, weighted by the predicted PC coefficient at that position. Lastly, we rescale the PPS residual model back to its original PSF size, as in <ref>. This PPS residuals model is added to the PPS model after the second interpolation, building our final PPS model.

Fig. <ref> compares the PPS residuals with the field distortion shears after the second and third interpolations. Although we observe a slight improvement after the third interpolation, this is not quite significant. We have explicitly tested various polynomial orders for our third interpolation, ranging from 2 to 6. The third order demonstrated the best performance. In addition, we have also tested our interpolation using random forest, but in this case it performs significantly worse than polynomial, likely due to overfitting.

We summarize all the steps of our PSF reconstruction scheme in a flowchart shown in fig. <ref>.




§ CONCLUSION
 
Point spread function (PSF) reconstruction is a crucial step towards accurate shear measurements, as it directly affects the observed galaxy morphology. Although current PSF reconstruction methods can capture most of the PSF features, there are often some PSF residuals remaining that limit the accuracy of, e.g., shear measurements. It is therefore essential to develop new techniques that can further improve current PSF reconstruction methods. In this work we introduce a novel PSF reconstruction method composed by three layers of interpolation, following a hierarchical scheme. Using data from the second data release of the Hyper Suprime-Cam (HSC DR2) as an example, our method significantly reduces the systematic PSF residuals, improving the PSF model.

The first layer of interpolation (<ref>) is a 2D polynomial fitting of third order to the power spectrum of the stars, pixel-by-pixel, on each CCD separately. This interpolation captures most of the spatial features of the PSF, building our initial PSF model.

The second layer of interpolation, described in <ref>, is the most important and novel part of this work. In this layer we model the systematic PSF residuals remaining after the first interpolation, using the data of all the exposures simultaneously. Note that to do so, we need to first uniform the image sizes of the PSF residuals. We apply PCA to the PSF residuals and build a model for each principal component coefficient. We then build the PSF residual model as a weighted sum of the principal components, with the weights being the predicted PC coefficients. To make corrections on the PSF model, we rescale the PSF residuals model back to the original PSF size of each exposure and add it to the original PSF model. Fig. <ref> and fig. <ref> show the PSF residuals of our improved PSF model, for polynomial and random forest interpolations, respectively. In both cases our method successfully removes most of the systematic features of the PSF residuals. In addition, we test our model directly on real galaxies, using the field distortion test (<ref>). We obtain significant improvements compared to the results of the first interpolation. Fig. <ref> shows the field distortion test results for the cases of polynomial and random forest interpolations. Both interpolations have comparable performances, leading to very small shear biases.

In <ref> we study a correlation between the PSF residuals and the field distortion signal, which is found after smoothing the PPS residuals after the second interpolation (see fig.<ref>). In this case, the systematic features are global, hence we perform a single interpolation on the exposure level, making use of all available exposures together. As in the second interpolation, we first rescale the PSF residuals to a common size. Then, we apply PCA to the entire data, and build a model for each PC coefficient. Our new PSF residuals model is built as a weighted sum of the principal components, with the weights given by the predicted PC coefficient at each position on the CCD. Lastly, we rescale the PSF residuals model back to the original PSF size of each exposure, and add it to the PSF model. Fig.<ref> presents the PPS residuals versus the field distortion shears, showing a small improvement after our third interpolation, although not quite significant. Nevertheless, this effect is minor and does not significantly impact the PSF model.

Overall, our model introduces a way to model the systematic PSF residuals, successfully removing most of the systematic PSF residuals in HSC DR2. In the current framework, our machine learning approach (random forest) performs similarly to polynomial, in the second interpolation. More exotic machine and deep learning algorithms may further improve the interpolation, reducing the PSF residuals.

In conclusion, our method aims to be a stepping stone towards building better PSF models, reducing the systematic biases induced by the PSF and helping produce more accurate and reliable shear catalogs, which are essential to the understanding of the distribution of dark matter, galaxy evolution, and to constrain cosmological parameters.



§ ACKNOWLEDGEMENTS

The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from the Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University. 

This paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at  http://dm.lsst.org

This paper is based [in part] on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by the Subaru Telescope and Astronomy Data Center (ADC) at National Astronomical Observatory of Japan. Data analysis was in part carried out with the cooperation of Center for Computational Astrophysics (CfCA), National Astronomical Observatory of Japan. The Subaru Telescope is honored and grateful for the opportunity of observing the Universe from Maunakea, which has the cultural, historical and natural significance in Hawaii.

JZ is supported by the National Key Basic Research and Development Program of China (2023YFA1607800, 2023YFA1607802), the NSFC grants (11621303, 11890691, 12073017), and the science research grants from China Manned Space Project (No. CMS-CSST-2021-A01). The computations in this paper were run on the π 2.0 cluster supported by the Center of High Performance Computing at Shanghai Jiaotong University, and the Gravity supercomputer of the Astronomy Department, Shanghai Jiaotong University.



aasjournal

























