Yubin Li ,Jirong Mao ,Jianbo Qin ,Xianzhong Zheng ,Fengshan Liu ,Yinghe Zhao ,and Xiaohong Zhao

1 Yunnan Observatories,Chinese Academy of Sciences,Kunming 650011,China;jirongmao@mail.ynao.ac.cn

2 Center for Astronomical Mega-Science,Chinese Academy of Sciences,Beijing 100101,China

3 University of Chinese Academy of Sciences,Beijing 100049,China

4 Key Laboratory for the Structure and Evolution of Celestial Objects,Chinese Academy of Sciences,Kunming 650011,China

5 Purple Mountain Observatory,Chinese Academy of Sciences,Nanjing 210023,China

6 National Astronomical Observatories,Chinese Academy of Sciences,Beijing 100101,China

Abstract GW170817 is the unique gravitational-wave (GW) event associated with the electromagnetic (EM) counterpart GRB 170 817A.NGC 4993 is identified as the host galaxy of GW170817/GRB 170 817A.In this paper,we focus on the spatially resolved properties of NGC 4993.We present the photometric results from the comprehensive data analysis of the high spatial-resolution images in the different optical bands.The morphological analysis reveals that NGC 4993 is a typical early-type galaxy without significant remnants of a major galaxy merger.The spatially resolved stellar population properties of NGC 4993 suggest that the galaxy center has passive evolution with the outskirt formed by gas accretion.We derive the merging rate of the compact object per galaxy by a co-evolution scenario of a supermassive black hole and its host galaxy.If the galaxy formation is at redshift 1.0,the merging rate per galaxy is from 3.2×10−4 to 7.7×10−5 within the merging decay time from 1.0 to 5.0 Gyr.The results provide vital information for ongoing GW EM counterpart detections.The Hubble space telescope data analysis presented in this paper can be also applied to Chinese Space Station Telescope research in the future.

Key words: galaxies: evolution–(stars:) binaries: general–gravitational waves

1.Introduction

High-frequency gravitational waves (GWs) originate from compact object mergers.Electromagnetic (EM) radiation is accompanied by gravitational-wave(GW)release(Nakar 2020).Short-duration gamma-ray bursts (GRBs) are produced by compact object mergers(Berger 2014).Thus,GRBs are usually considered as a souce of GWs.GW170817 is the unique source of a GW that has a confirmed EM counterpart,and GRB 170 817A accompanied by GW170817 has multi-wavelength observation.It has been confirmed by both GW and multiwavelength EM data that GW170817/GRB 170 817A arose from a neutron star merger (Abbott et al.2017).It is also important to note that the environment of the compact object merger plays a vital role in both the occurrence and the evolution of a GW source with an EM counterpart.Fortunately,GW170817 occurred only 41 Mpc from Earth.Being so close,the host galaxy of the source,named NGC 4993,was clearly recognized (Hjorth et al.2017).

It has been found that many host galaxies of short GRBs usually follow the sequence of the star-forming galaxies(D’Avanzo et al.2009;Leibler &Berger 2010;Fong et al.2013).However,NGC 4993 is an early-type galaxy and its shape is symmetric,making it appropriate to classify it as S0 style (Palmese et al.2017).The galaxy has a stellar mass of log(M/M⊙)=10.49 and its star-formation rate is as low as 0.003M⊙per year(Pan et al.2017).This is a very unusual case in the catalog of short GRB host galaxies(Nugent et al.2022).Moreover,the compact merger system provides a kind of nucleosynthesis process called the rapid neutron-capture process (r-process).By this process,elements heavier than iron can be produced during the merging time.Ther-process elements were successfully examined at the merging site in NGC 4993 (Pian et al.2017).However,ther-process source may have a delayed timescale of larger than 4 Gyr and binary merger may not be the only way that ther-process is produced in NGC 4993(Skúladóttir&Salvadori 2020).It is necessary to comprehensively perform photometric analysis and obtain the global properties of the GW170817 host galaxy.In order to further investigate the physical properties in the neighborhood of the merging source compared to those of the whole galaxy,spatially resolved photometric measurements of the GW170817 host galaxy are required.Furthermore,for the binary merger case in NGC 4993,spectral analysis has been well performed both at the merger region and in the galactic center (Blanchard et al.2017;Levan et al.2017;Pian et al.2017).However,in principle,physical information for each location inside NGC 4993 needs to be provided.Although an integral field unit can be adopted to obtain the two-dimensional spectrum of a galaxy,photometric measurements are still required in order to study faint GW host galaxies.One may perform a direct GW EM counterpart search in the high-energy band,but the GW location has a huge error circle.The optical band appears to be a more favorable choice for identifying potential host galaxies in the search for GW EM counterparts.The ranked host candidates can be selected by semi-analytic methods or full simulations from galaxy evolution models(Mandhai et al.2022;Perna et al.2022).GRB 170 817A,as the EM counterpart of GW170817,is the only GW EM counterpart that has been identified so far and detailed investigations of the photometric properties of NGC 4993 can provide important information for GW EM follow-up observations in the future.For example,some strategies to select GW EM counterparts among many celestial candidates were established by accurate photometric observations to obtain galaxy properties (Ducoin et al.2020).

In general,studies of the compact merger system can be put in the framework of galaxy formation and evolution in the universe (Gehrels et al.2016;Toffano et al.2019;Adhikari et al.2020).Some host galaxy properties of merging objects,such as star-formation rate and stellar mass,are related to the merging rate per galaxy (Artale et al.2019).The delay time from binary formation to binary merger is dependent on the host galaxy properties (Mapelli et al.2018;Safarzadeh &Berger 2019;McCarthy et al.2020).The star-formation history of host galaxies can be adopted to estimate the merger formation rate (Rose et al.2021).Thus,a compact merger system is naturally linked to the formation and evolution of its host galaxy.For the host galaxy of GW170817,a shell-like structure has been identified in NGC 4993,and this indicates that NGC 4993 was formed by a galactic merger from 400 Myr ago(Ebrová et al.2020).In addition,we note that the accretion activity in the center of NGC 4993 shows some typical features of a low-luminousity active galactic nucleus(Contini 2018;Wu et al.2018).External gas accretion may supply a S0 galaxy(Raimundo 2021).Moreover,many globular clusters have been identified in NGC 4993 (Lee et al.2018).As the globular clusters have a long lifetime,we may expect that the budge of NGC 4993 may exempt it from major galactic mergers in the past evolution.If NGC 4993 has a monolithic evolution over a long cosmic time,compared to the galaxy merger evolution,the monolithic evolution mode may have different effects on the occurrence of the compact merger system.

In this paper,we use data from the Hubble space telescope(HST) observation to comprehensively analyze the properties of NGC 4993.Although some works on the global properties of NGC 4993 have been performed using HST data (Palmese et al.2017;Ebrová et al.2020;Kilpatrick et al.2022),we focus on the spatially identified images of NGC 4993 in this paper.The stellar properties can be effectively derived from the multiband images.Levan et al.(2017) obtained the properties of NGC 4993 from integral field spectroscopy observations.However,we note that space/ground-based spectroscopy observation is only suitable for bright sources.When attempting to identify faint sources during GW EM counterpart searching,photometric observation is almost the only way to investigate the host galaxy properties of the GW EM counterpart candidates,as spectroscopy observation is hard to perform for faint targets.In this paper,we independently perform image analysis of the Hubble observation.The method can be applied to future GW EM counterpart searches.We expect that the research on the HST data presented in this paper is also helpful for China Space Station Telescope (CSST) observations in the future.The data from Pan-STARRS are also considered as a comparison to the data from HST.The near-infrared (IR) data are also provided in this paper.To avoid duplication of some works from VLT and Gemini observations,we collect the photometric data from the 2MASS survey.Although the 2MASS survey is relatively shallow,it is deep enough to provide us with IR information for NGC 4993.

Spatially resolved properties,such as color diagrams and stellar populations,are presented in Section 2.The binary merging rate is estimated in Section 3.We draw a simple conclusion in Section 4.Details on the spatially resolved data analysis are presented in the Appendix.

2.Data Analysis

2.1.Profile Fitting

To quantitatively characterize the photometric structure of NGC 4993,the Sérsic profile of the galaxy can be measured by GALFIT (Peng et al.2002).The results can be obtained in bands spanning optical to the near-infrared wavelength.We take the images of theg,r,i,zandybands from the Pan-STARRS survey.The F606W image of the Hubble legacy observation is also selected as a reference.It is shown that NGC 4993 has a Sérsic index of 4 fitted by the Sérsic profile and it is proved that the galaxy is an early-type galaxy dominated by a bulge component.We also see that the Sérsic index is increased as the observational wavelength is increased.The detailed results are listed in Table 1.The results of our measurements are consistent to those of other works (Palmese et al.2017).We further identify that the effective radius is aboutwhich corresponds to a linear distance of 3 kpc.According to the stellar mass of(Pan et al.2017),NGC 4993 is slightly above but consistent with the mass–size relation of local early-type galaxies.The image of NGC 4993 can be subtracted by the fitting model.As an example,Figure 1 shows both the original image and the residual image after fitting in therband.Some substructures,such as dust lanes and shells,are shown in the residual image.Here,we suggest that the bulge has early star formation during passive evolution and that the outskirt has late star formation due to external gas accretion.Although the substructures indicate that some possible merger activities happened in the past,we may consider the substructures as evidence of gas accretion.

Table 1The magnitudes of NGC 4993 given by GALFIT fitting

Figure 1.The image of NGC 4993 in the r band(left panel)and the residual image of NGC 4993 after image subtraction by GALFIT(right panel).The symbol “+” represents the center of the galaxy and “×” represents the position of GW170817.

Figure 2.The two-dimensional distribution of the optical and near-infrared colors obtained by Pan-STARRS and 2MASS survey for NGC 4993.The symbol “+” represents the center of the galaxy and “×” represents the position of GW170817.

Figure 4.The one-dimensional optical and near-infrared color profiles of NGC 4993 derived using images from the Pan-STARRS and 2MASS surveys.The black line represents the color profile of NGC 4993 and the red point represents the color at the position of GW170817.

Figure 5.The one-dimensional optical and near-infrared color profiles of NGC 4993 derived using images obtained by the HST.The black line represents the color profile of NGC 4993 and the red point represents the color at the position of GW170817.

Figure 6.The two-dimensional distributions of age,SFR,mass,sSFR,Av and metallicity in NGC 4993.The asymmetrical distribution of metallicity is shown in panel(f).In our work,we have derived the metallicity error of about 0.01.This error is as large as the metallicity difference between the center and the outskirt,and it effects the metallicity distribution.

Figure 7.The one-dimensional profiles of age,SFR,mass,sSFR, Av and metallicity in NGC 4993.

2.2.Color Diagrams

Detailed investigations by photometric analysis are crucial to understand stellar population properties in a galaxy.Here,we use Pan-STARRS,2MASS and HST data to obtain the twodimensional color distribution of NGC 4993.The galactic attenuation for the colormaps is corrected (Schlafly &Finkbeiner 2011).The results are shown in Figures 2 and 3.It is clearly seen that NGC 4993 has a red center and a blue outskirt.Furthermore,the color at the position of GW170817 seems very similar to that in the adjacent regions.It means that the stellar population underlying the position where GW170817 occurs is identical to the stellar population of the GW170817 neighborhood.The local environment of GW170817 represents the common properties of the outskirt in NGC 4993.

We can further identify the one-dimensional color profiles of NGC 4993.We utilize elliptical annulus photometry to get the surface brightness profile in each band.The geometry parameters,such as center,ellipticity and radius of each annulus,depend on the structure and the depth of the images.When producing the color profile,the definition for both the center and the ellipticity of the annuli should be clarified.The center and the related ellipticity of the galaxy are not completely identical in different bands.Thus,when producing the color profile,the center and the ellipticity (0.166) are all fixed for the GALFIT fitting using the measurement result of the F160W band.The radius along the major axis has the range ofin the HST images and ofin the Pan-STARRS images.When producing the one-dimensional surface brightness profile,the radii of the annuli are logarithmic increased instead of linear.As the radius increases,although the signal-to-noise of a single pixel decreases,more pixels are integrated,leading to the similarity of the signal-to-noise from the inner regions to the outskirts.Thus,the signal-to-noise ratio of the photometry is similar for each annulus at different radius.After the elliptical annulus photometry is performed in each band,we obtain the surface brightness profile and the onedimensional color profile.The results are shown in Table 2 and Figures 4 and 5.It can be seen that,generally,NGC 4993 has negative color gradients and it is confirmed that NGC 4993 has a red center and a blue outskirt.Moreover,the color gradients become shallow as the radius is increased.It is indicated that the formation process for the galactic core is different to that of the galactic outskirt.Here,we identify that the inner region of the galaxy is within 0.5Reand the outskirt region of the galaxy is beyond 0.5Re,whereReis the effective radius of the galaxy.

Table 2The color gradients of NGC 4993

2.3.Stellar Population

Usually the HST has a better spatial resolution than groundbased optical telescopes.We obtain the spectral energy distribution (SED) composed of fluxes in the HST F606W,F814W,F110W,F140W and F160W bands in each pixel of NGC 4993.we fit these SEDs using the Code Investigating GALaxy Emission (CIGALE Noll et al.2009;Boquien et al.2019) to obtain the galaxy properties on a resolved scale in NGC 4993.CIGALE combines a library of single stellar populations and variable attenuation curves with Star Formation History (SFH) models to generate a large number of grid SED models to fit the observed data.To build a stellar composition,we use the BC03 stellar population synthesis model(Bruzual&Charlot 2003)with the Chabrier initial mass function (Chabrier 2003).The metallicity number that we take is from 0.2 to 2.5Z⊙(forZ⊙=0.02).A delayed star-formation historyis adopted,where t is the stellar age(varies from 1 to 13 Gyr)and τ is the e-folding time(varies from 0.1 to 11 Gyr).For the dust attenuation,we adopt a fixed Calzetti attenuation curve (Calzetti et al.2000) with E(B-V) varying from 0.0 to 0.3 mag.The nebular emission is also included in our SED fitting.All the modules and parameters are summarized in Table 3.CIGALE makes use of flat priors.The best-fitting parameters and the corresponding uncertainties are the likelihood-weighted mean and the standard deviation of all models,based on the probability distribution functions generated with Markov chain Monte Carlo sampling.We finally obtain the distribution of the age,the Star Formation Rate(SFR),the sSFR,the extinction and the metallicity of NGC 4993.The two-dimensional distribution and the one-dimensional profiles are shown in Figures 6 and 7,respectively.Finally,the spatially resolved properties of NGC 4993 are clearly shown.

Table 3Modules and Input Parameters with CIGALE for Generating our Model SEDs

In order to illustrate the properties of NGC 4993 compared to those of other normal galaxies,we utilize the data set of the Sloan Digital Sky Survey (SDSS).We adopt the DR10 value-added catalog including the spectroscopic line measurements as well as the galaxy parameters derived from theu,g,r,iandzimaging and spectroscopic data.We limit our sample to the redshift range of 0 ＜z＜0.02,as NGC 4993 is at redshift of about 0.01.A lower limit ofM*≥108M⊙is adopted to ensure the mass completeness in this redshift range.We also removed sources with insecure redshift measurements as well as galaxies spectroscopically classified as quasi-stellar objects.There were 8629 sources left after applying these selection criteria.Figure 8 shows theg−rcolor as a function of the stellar mass.NGC 4993 follows the normal properties of the red and early-type galaxies in the figure.It seems that providing the ranked GW host galaxies through galaxy property investigations is an important aspect of GW EM counterpart detection.

Figure 8.Background gray contour shows the g−r color as a function of the stellar mass for the galaxies in the SDSS sample with z ＜0.02.NGC 4993 is labeled as the blue star.The dashed line separates the galaxy sample into the red data set and the blue data set..

3.Modeling Estimation

The merging rate of compact objects in a galaxy is the result of the convolution of the star-formation rate in the galaxy with a certain delay time.The delay time is defined as the duration from binary formation to merger occurrence.The merging rate of compact objects per galaxy is calculated by

where φ(t) is the star-formation rate,dp/dt(t−td) is the delay time distribution,andtdis the minimum delay time.Usually,we have(Mapelli et al.2018;Safarzadeh &Berger 2019;Adhikari et al.2020).Here,we assume that the delay time distribution can be a δ-function asdp/dt=δ(t−td),and we obtainR=λφ(td).The case ofdp/dt=δ(t−td)indicates that binaries have an instant merging after we take the delay time oftd(Adhikari et al.2020).Here,we usethat was given by Safarzadeh &Berger(2019).Some globular clusters in NGC 4993 have been identified by deep photometric measurements(Lee et al.2018),although GW170817 did not have a globular cluster origin(Fong et al.2019).It is indicated that the environment of GW170817 was not strongly affected by the galactic merger during the galaxy evolution time.Therefore,we may consider NGC 4993 as a giant S0 galaxy with a passive evolution.

The calculation is in the framework of the co-evolution between the central black hole and its host galaxy.We utilize the model of the co-evolution between a central black hole and its host galaxy provided by Granato et al.(2004).The semianalytic model can provide the star-formation processes of a galaxy in a given dark matter halo at a certain redshift.The feedback from both supernovae and active galactic nuclei is involved in the model.This model has wide applications.For example,a careful investigation on high-redshift star formation and absorption had been performed (Mao et al.2007).In particular,Mao et al.(2010) utilized the model to investigate the physical properties of long-duration GRB host galaxies.The star-formation process of a galaxy in a certain dark matter halo at a certain redshift is affected by active galactic nucleus feedback and the supernova.Here,we can apply this model to further constrain the merger rate in NGC 4993.The starformation rate in a galaxy is

wherem(0)=18%MH,MHis the dark matter halo mass and we set the parameters=5.The condensation timescale for gas conversion into stars in a dark matter halo at a given redshift can be presented as

whereESNis the energy release from supernova and we take the parameterf=0.3.

We assume that the merging occurs immediately after the delay time oftd.If the host galaxy was formed at redshift 1.0 according to the mass-weighted stellar age of about 8.0 Gyr and it has a passive evolution,we obtain the merging rate per galaxy of 3.2×10−4when we take the delay timetd=1.0 Gyr.The merging rate per galaxy can be about 7.7×10−5when we take the delay timetdto be 5.0 Gyr.It is confirmed that compact object merging in a galaxy is a rare event in the universe.This constraint provides a valuable reference for future GW EM counterpart identification studies.

4.Conclusion

We comprehensively perform photometric analysis of NGC 4993.The spatially resolved properties of the galaxy are clearly presented.Although the shell of NGC 4993 was identified as evidence for galaxy merging,the mass of the shell seems too small to be the product of a major galaxy merger (Kilpatrick et al.2022).We suggest that the galaxy center has passive evolution and the outskirt is formed by gas accretion.We estimate the compact binary merging rate per galaxy as 3.2×10−4to 7.7×10−5within the merging decay time from 1.0 to 5.0 Gyr.The methods used for the spatially resolved data analysis and the physical constraints on the binary merging in a galaxy are very useful for future GW EM counterpart detections.The HST data analysis presented in this paper can be also applied to CSST research in the future.

Acknowledgments

All the HST and Pan-STARRS data used in this paper can be found in MAST: data set 10.17909/T97P46,10.17909/55e7-5x63,and 10.17909/s0zg-jx37.This work is supported by the National Science Foundation of China (NSFC 11673062),the China Manned Space Project (CMS-CSST-2021-A06) and the Yunnan Revitalization Talent Support Program (YunLing Scholar Award).J.Q.acknowledges support from the Jiangsu Funding Program for Excellent Postdoctoral Talent (No.2022ZB473).X.Z.Z.is grateful for support from the NSFC(11773076 and 12073078),the National Key R&D Program of China (2017YFA0402703) and science research grants from the China Manned Space Project(Nos.CMS-CSST-2021-A02,CMS-CSST-2021-A04 and CMS-CSST-2021-A07).F.L.is grateful for support from the NSFC(11733006 and 12273052).Y.Z.is grateful for support from the NSFC(12173079).X.Z is grateful for support from the NSFC (U1831135).

Appendix A Image Selection and Data Reduction

We use both optical and near-IR images to get the spatially resolved properties of NGC 4993.The images are obtained from the Pan-STARRS survey (Chambers et al.2016),the 2MASS survey (Skrutskie et al.2006) and the HST legacy survey (Alexander et al.2018;Lyman et al.2018;Margutti et al.2018;Lamb et al.2019;Piro et al.2019).The reduction of these ground-based and space-based data are summarized as follows.

The images obtained by the Pan-STARRS survey have five bands (g,r,i,zandy) and cover the wavelength range from 0.45 to 1 μm.The pixel size isand the spatial resolution ranges fromindicating that the 250 pc scale substructure of NGC 4993 can be identified if we accept 41 Mpc as the distance of NGC 4993 (Cantiello et al.2018).

We also obtained 2MASS images of NGC 4993 in the three near-IR bands (namedJ,HandK).The pixel size isand the spatial resolution iscorresponding to a 600 pc scale substructure for NGC 4993.The morphology of NGC 4993 is more representative for the stellar mass distribution in the NIR bands than that in the optical bands.Furthermore,with the NIR observations,the spectral energy distribution (SED) can be extended toward longer wavelength,and this can give us more information on the distribution of the stellar population.

The images obtained by the HST have ultra-high spatial resolution due to the absence of air turbulence.The pixel size ranges fromfor the optical filter WFC3/F814W tofor the NIR filters(WFC3/F110W,WFC3/F140W and WFC3/F160W).The spatial resolution ranges fromtofor the optical bands and fromfor the NIR bands,corresponding to 30–80 pc scale substructure for NGC 4993.The detailed morphology and the distribution of the stellar population for NGC 4993 can be recovered by analysis of the HST images and it benefits from an unprecedented high spatial resolution.To study the underlying stellar populations of NGC 4993,especially the stellar populations at the position of GW170817,the images of NGC 4993 obtained before the GW170817 occurrence are better to be used to avoid contamination from the afterglow of GRB 170 817A.However,before the GW170817 occurrence,NGC 4993 was only observed by ACS/F606W as part of the Schedule Gap Pilot program (PI: Andrea Bellini ID:14840) observed on 28 April 2017 and the exposure time is 696 s.It is challenging to obtain the distribution of the stellar population within the galaxy through a single-band image.After the GW170817 occurrence,NGC 4993 was monitored by the HST at different bands and so multi-band images are available for the purpose of the stellar population synthesis.We select the images obtained at least 50 days after the GW170817 trigger to ensure that the light of NGC 4993 is dominated by the stellar populations.At this stage,the images of NGC 4993 are only contaminated by the afterglow of the structured jet in GRB 170 817A,which is fainter than 26 mag,4 mag dimmer than the kilonova (Fong et al.2019).Therefore,the properties of the stellar population at the position of GW170817 can be precisely inferred.Table A1 lists the information about the images of NGC 4993 observed by the HST that we use in this work,including observing time,exposure time,instruments,filters and proposal ID.Briefly,the images obtained by the HST cover five optical and NIR bands (F606W/F814W/F110W/F140W/F160W),which enables us to get the distribution of the stellar population with ultra-high resolution.We then stack the images obtained in the same band to enhance the signal-to-noise ratio.The stacking process was as follows: (1) remove cosmic rays and bad pixels in the images,which have been bias-subtracted and flatfield-divided.(2) Subtract the background for the images.(3) Match the astrometry of the images using the astrometry information listed in the head file to the same reference image;this reference image was obtained on 2017-12-08 20:33:09 in the F110W filter.(4) Extract the empirical point-spread function (PSF) for all the images by stacking the stars in each corresponding image.(5) Match all the PSFs to the worst/largest PSF and the PSF obtained from the image observed on 2017-12-06 01:45:51 in the F160W filter with an FWHM(full width at half maximum)ofThis ensures that the value of each corresponding pixel grid at the same position is a weighted average of the neighboring pixels produced in the same way.Thus,the same pixel grid represents the same region in the galaxy,such that we can compare the same region in the different images.(6)Combine the astrometry and PSF-matched images obtained in the same band using the exposure time as the weight.The value at each pixel is calculated by

wheretexpis the exposure time,(x,y) is the coordinate in the image and DN is the digital number at a pixel.Since the DN is in the unit of electrons s−1in the calibrated images obtained by HST,the final combined image is still in the unit of electrons s−1,and the photometric zero-point of the combined image is unchanged.Finally,we get the stacked images of NGC 4993 in the five bands(F606W,F814W,F110W,F140W and F160W),respectively.After the stacked multi-band images are recovered,we use the Sérsic profile to describe the morphology of the galaxy in the different bands and we also get the colormap and the color profile of the galaxy from the images in the different bands.

Table A1Summary of the images obtained by the HST

We then focus on reduction of the above HST data.The HST images are already drizzled products that can be downloaded from the STSCI website.Then we shift the different images to a fixed reference coordinate and stack the images at the same band using the exposure time as the weight to increase the signal-to-noise ratio.Before stacking the images at the same band,the PSF was also matched to ensure that the images to be added have identical PSFs.

Appendix B The Extraction and Matching of the PSF

We normally use PSF to describe the response of an optical system for pointed sources.In principle,an image of a celestial source is the convolution of the PSF with the intrinsic intensity distribution.Thus,understand the intrinsic surface brightness profile of a galaxy,we must obtain the PSF in the image.The existence of the PSF also means that the value at each pixel is the weighted average of the pixels around it (including the pixel itself,which usually has the maximal weight).Different images usually have different PSFs due to the variance of the air turbulence and the variance of the filters with the different observing strategies.The air turbulence is exempted if we take images from space-based telescopes.If the astrometry is matched between different images,the pixel at the same position is indeed from different regions,as the adjacent pixels contribute different weights to the same pixel.In order to do photometry at the same pixel in different images directly,we should match both the astrometry and the PSF of the images before stacking the images and producing colormaps or color profiles.To match the PSF,we first extract the PSF in each image.The pointed sources in each image are adopted.When extracting the PSF,the empirical method is used to select the unsaturated pointed sources without contamination from the neighboring sources and to stack the images of the pointed sources weighted by the flux of each pointed source.

For the images obtained by the Pan-STARRS survey,we use the following criteria to select the pointed sources.(1) We select the sources and each source has the photometric difference between PSFMag and KronMag less than 0.05 mag iniband.Here,PSFMag is the magnitude obtained from the PSF profile and KronMag is the magnitude obtained from the Kron radius.This selection condition was also used to distinguish the pointed source and the extended source(Farrow et al.2014).Pointed sources like stars can be well described by a PSF profile and the PSFMag value is close to the KronMag value.For extended sources like galaxies,the surface brightness profile cannot be described by a PSF profile,and the magnitude obtained by PSFMag is dimmer than that obtained by KronMag.(2) We plot all the detected sources on the magnitude-half-light radius diagram as shown in Figure B1.The pointed sources have a similar half-light radius,which is the half-light radius of the PSF profile,while the extended sources have a larger half-light radius.Thus,the pointed sources and the extended sources lie in the different regions on the magnitude-half-light radius diagram,respectively.This is the algorithm used by PSFExtractor (Bertin 2013).We select the sources with half-light radius less than 5 pixels (corresponding toand with magnitude brighter than 20 mag in the Pan-STARRS image.(3) We also use the parameter CLASS_-STAR derived by SExtractor to select the pointed sources(Bertin &Arnouts 1996).We consider sources that meet the condition ofCLASS_STAR≥0.9.Finally,we select the pointed sources that satisfy all the three criteria above to extract PSFs for the images in the different bands.Sources with contamination from neighboring ones and saturated sources are removed.

Figure B1.The sources on the magnitude-half-light radius diagram.The red dots represent the sources with CLASS_STAR ≥0.9 and the blue circles represent the sources with CLASS_STAR ＜0.9.We can see that the pointed sources and the extended sources can be clearly distinguished in the diagram.

For the images obtained by the HST observation,since the field of view is not large enoughto get the statistical properties of the sources in the images,the pointed sources are visually selected.These selected sources are shown in Figure B2.After the sources are selected,we use the DAOPHOT package in IRAF to get the average PSFs in the different bands.Once the PSF is obtained,the best parameterized Sérsic profile convolving with the PSF for the image in each band can be achieved.Moreover,when stacking the images in each band and producing the colormaps and the color profiles among the different bands,the PSFs should be matched to ensure that the corresponding pixel in different images can be compared.We use the PSFMATCH package in IRAF to match all the PSFs to the worst PSF.The algorithm is described as

Where inten represents the intrinsic intensity distribution without the effect of PSF.If PSF2 is worse than PSF1,we find a kernel satisfying the following equation

Figure B2.The selected pointed sources used to extract the PSF of images obtained by the HST.These sources are unsaturated,relatively bright and devoid of contamination by their neighboring sources.

After the kernel is found,image1 will be convolved with this kernel as

After convolving with this kernel,the PSF is identical for both image1 and image2.Following Equation (B3),the method to find the kernel is presented as

Figure B3.The images of the PSFs before (left panels) and after (right panels) the PSF-matching procedure in the F606W and the F160W bands as an example.

whereF is the Fourier transformation andF-1the inverse Fourier transformation.The noise arising in the empirical PSF by stacking the images can influence this PSF-matching procedure severely.The effect of the noise in the matching procedure needs to be reduced.After the Fourier transformation is performed,the random noise signals are shown in the highfrequency range in the spatial-frequency domain.We test several different methods to reduce the high-frequency noise.For the Pan-STARRS images,we fit the low-frequency and high signal-to-noise ratio components of the matching function with a Gaussian model and apply this Gaussian model to replace the entire PSF function,following the algorithm replace with the parameter filter of the PSFMATCH package in IRAF.For the HST images,a cosine bell function is applied to the PSF-matching function in spatial-frequency space,which reduces the weight of the high-frequency component,following the algorithm cosbell with the parameter filter of the PSFMATCH package in IRAF.This algorithm can also match the direction of the asterism and it is suitable for HST image reduction.In practice,the high-frequency component mentioned above contributes to the center of the PSF.When we adopt the procedure mentioned above to reduce the noise,the matched PSF is slightly larger than the original PSF.To solve this problem,we match the worst PSF to the original PSF by using the same algorithm.After this additional process,the PSFs in all images are almost self-consistent.The FWHM numbers before and after PSF-matching procedure are listed in Table B1.As an example,Figure B3 shows the PSFs of images F606W and F160W before and after then PSF-matching procedure.We can see that not only the FWHMs of the PSFs but also the shapes (including the asterism) of the PSFs are almost identical after the PSF-matching procedure.Figures B4 and B5 show the PSF curves before and after the PSF-matching procedure.These panels indicate that the PSF-matching procedure in this work is reliable.

For Pan-STARRS images,before calculating the color gradients and the color disribution for each color index,the PSF was matched to the larger PSF.The PSF isandfor theg,r,i,zandybands,respectively.When producing theg−rcolormap,for example,since the PSF of theg-band image isand the PSF of ther-band image istheg-band image was matched to ther-band image and the g-r colormap has a resolution of

Figure B4.The aperture-dependent PSF before and after PSF-matching procedure of each image in each band.The blue line shows the aperture-dependent PSF before the PSF-matching procedure and the red line shows the aperture-dependent PSF after the PSF-matching procedure.The black dashed line shows the reference PSF to be matched.The consistency of the PSF after PSF-matching procedure and the reference PSF indicates that the PSF-matching procedure is reliable.

For HST images,in the PSF-matching procedure,all the HST images were matched to the HST image that had the worst/largest PSF (the image of F160W obtained at 2017-12-06 01:45:51).The FWHM of the worst PSF isas mentioned in Appendix A.When removing the noise of the PSF,the center of the PSF was also smoothed since it has a high-frequency component in the spatial-frequency domain.Thus,the ultimate resolution of the HST images used in the analysis isslightly larger than the reference PSF.

For 2MASS images,the FWHMs of PSFs atJ,HandKbands were almost identical atand the PSF-matching procedure can be ignored when we produceJ−H,J−KandH−Kcolor maps.

In addition,there are a few stars in the observational field that overlap the galaxy.We simply mask the regions contaminated by stars using the segmentation map provided by the SExtractor.Thus,our results are not affected by the overlapping stars.

Figure B5.Continued Figure B4.

Table B1FWHM comparison of the original image and the image after the PSF-matching procedure in each HST Band

ORCID iDs