The First LHAASO Catalog of Gamma-Ray Sources

We present the first catalog of very-high energy and ultra-high energy gamma-ray sources detected by the Large High Altitude Air Shower Observatory (LHAASO). The catalog was compiled using 508 days of data collected by the Water Cherenkov Detector Array (WCDA) from March 2021 to September 2022 and 933 days of data recorded by

the Kilometer Squared Array (KM2A) from January 2020 to September 2022.This catalog represents the main result from the most sensitive large coverage gamma-ray survey of the sky above 1 TeV, covering declination from −20 • to 80 • .In total, the catalog contains 90 sources with an extended size smaller than 2 • and a significance of detection at > 5σ.Based on our source association criteria, 32 new TeV sources are proposed in this study.Among the 90 sources, 43 sources are detected with ultra-high energy (E > 100 TeV) emission at > 4σ significance level.We provide the position, extension, and spectral characteristics of all the sources in this catalog.

INTRODUCTION
Gamma rays, located at the highest energy band of the cosmic electromagnetic radiation, serve as a powerful probe for astrophysics and fundamental physics under extreme conditions.Most gamma rays are produced through the acceleration and propagation of relativistic particles, such as protons or electrons, in astrophysical sources.Relativistic electrons can scatter low-energy photons, including star light, dust scattered light and the Cosmic Microwave Background (CMB), to the gamma-ray band via the inverse Compton process.The whole universe is filled with low energy photons, especially the CMB.Therefore, gamma rays are a unique tool to probe the relativistic electrons when the surrounding magnetic field is weak and the synchrotron radiation is undetectable.Relativistic protons interact with the surrounding medium to create hadronic cascades, which include secondary π 0 mesons that quickly decay into gamma-rays.Hence, gamma rays are also an important tool to study the origin, acceleration and propagation of cosmic rays (CRs).
Thanks to the advancements in space-based and ground-based gamma-ray detectors, our knowledge about the high energy gamma-ray universe has made impressive progress over the past two decades.At high energy (HE, E > 0.1 GeV), the Fermi Large Area Telescope (Fermi-LAT, Atwood et al. 2009) has been surveying the whole sky since 2008 and detected 6658 Galactic and extragalactic gamma-ray sources using the first twelve years of observations (Abdollahi et al. 2022).Compared to its predecessor, EGRET (Hartman et al. 1999), the number of detected sources has increased by a factor of more than 20, and some new categories of gamma-ray emitters have been revealed.Important evidence about the acceleration of GeV-TeV CRs were found in two ancient supernova remnants (Ackermann et al. 2013).The Dark Matter Particle Explorer (DAMPE) collaboration also reported the detection of more than 200 gamma-ray sources above GeV (Duan et al. 2021).
At Very High Energy (VHE, E > 0.1 TeV), ground-based gamma-ray detectors are necessary due to their large area requirements.The successful operation of the second generation Imaging Atmospheric Cherenkov Telescopes (IACTs), such as H.E.S.S. (Aharonian et al. 2006), MAGIC (Aleksić et al. 2016), and VERITAS (Meagher & VERITAS Collaboration 2015), has significantly increased the number of detected VHE gamma-ray sources from about 10 to over 200 since 2004.Several categories of VHE gamma-ray emitters have been firmly established, including Active Galactic Nuclei (AGNs), pulsars and Pulsar Wind Nebulae (PWNe), Supernova Remnants (SNRs), binaries, starburst galaxies, Gamma-Gay Bursts (GRBs) and others.However, due to their narrow Field Of View (FOV, 3 • − 5 • ) and low duty cycle (10% −15%), IACTs are not suitable for performing long-term comprehensive sky surveys.Most VHE sources are detected while searching for counterparts of sources observed at lower energies, and only certain parts of galactic plane have been surveyed by the IACTs, such as H.E.S.S. (H.E. S. S. Collaboration 2018) and VERITAS (Staszak & VERITAS Collaboration 2015).
To achieve an overall view of the VHE universe, a roughly unbiased sky survey is needed, similar to that carried out by Fermi-LAT in the GeV band.Ground-based Extensive Air Shower (EAS) arrays, with their large field of view and high duty cycle, are ideal detectors for this scientific goal.Several VHE gamma-ray surveys with positive results have been conducted to date, including those by Tibet ASγ (Amenomori et al. 2005), Milagro (Abdo et al. 2007), and ARGO-YBJ (Bartoli et al. 2013).However, due to the limitations of detector sensitivity, only a handful of sources have been observed.Nonetheless, impressive progress has been made in observing some typical extended sources (Bartoli et al. 2014) and variable AGNs (Bartoli et al. 2016(Bartoli et al. , 2012)), highlighting the invaluable role of EAS arrays in VHE gamma-ray observation.Recently, the sensitivity of EAS arrays has greatly improved thanks to the successful operation of the new generation arrays, such as HAWC and LHAASO.The HAWC Observatory has detected 65 VHE sources including 20 new ones using five years of data in their third catalog (3HWC; Albert et al. 2020).In particular, HAWC has revealed a new source category,the TeV pulsar halo (Abeysekara et al. 2017), which is a useful tool to probe the CR diffusion in the interstellar medium (ISM) near PWNe.
Another crucial characteristic of the EAS array is that it can extend the observation to the Ultra-High Energy (UHE, E > 0.1 PeV) range due to its large detector area and long duty cycle.The Tibet ASγ collaboration first reported a UHE gamma-ray source, the Crab Nebula (Amenomori et al. 2019), followed by another three sources reported by the HAWC collaboration (Abeysekara et al. 2020).Recently, the LHAASO collaboration has made exciting progress by reporting the detection of 12 UHE gamma-ray sources with a statistical significance of over 7σ and the maximum energy up to 1.4 PeV (Cao et al. 2021a).Additionally, some sources were observed with VHE emission for the first time (Cao et al. 2021b,c).These observations provide crucial candidates to explore the origin of PeV CRs within the Galaxy.LHAASO also detected PeV gamma-ray emission from the Crab Nebula, revealing an extreme electron accelerator with unprecedented high accelerating efficiency (LHAASO Collaboration 2021).These observation involving the highest gamma-ray energy also shed light on exploring the Lorentz Invariance Violation (Cao et al. 2022a) and dark matter (Cao et al. 2022b).It is worth noting that these results were achieved using only 1 year of data and half of the LHAASO array prior to completion of its construction.
This paper is structured as follows.Section 2 briefly describes the LHAASO detector and the data set.Additionally, it presents the background estimation and significance sky maps.In Section 3, the methods of searching for sources and constructing the catalog are introduced.The characteristics of WCDA and KM2A source components are briefly described and compared.In Section 4, a final source catalog is compiled, including relevant positional, spatial, and spectral information.The results of some typical source populations are discussed in Section 5. Section 6 provides a summary.

The LHAASO Detector and Data
LHAASO is a multi-purpose and comprehensive EAS array, designed for the study of CRs and gamma rays across a wide energy range, from sub-TeV to beyond 1 PeV (Ma et al. 2022).It consists of three interconnected detector arrays, a 1.3 km 2 array (KM2A) for gamma-ray detection above 10 TeV, a 78,000 m 2 Water Cherenkov Detector Array (WCDA) for TeV gamma-ray detection, and a Wide Field-of-view Cherenkov Telescopes Array (WFCTA) mainly for CR physics.When a highenergy extraterrestrial particle, gamma ray or CR, enters Earth's atmosphere, it initiates a cascade consisting of secondary hadrons, muons, leptons, and photons known as an air shower.The WCDA and KM2A detectors record different components of these air showers, which are used to reconstruct the types, energies, and arrival directions of the primary particles.The WCDA consists of three ponds: WCDA-1 and WCDA-2, both measuring 150 m × 150 m , and WCDA-3 measuring 300 m × 110 m.The total area of the array is 78,000 m 2 .WCDA-1, WCDA-2, WCDA-3 are composed of 900, 900, 1320 detector units, respectively.Each detector unit is 5 m × 5 m separated by nonreflecting black plastic curtains and is equipped with two upward-facing PMTs (8-inch and 1.5 inch PMT combination) on the bottom at the center of the unit.To further lower the threshold energy, WCDA-2 and WCDA-3 employ a 20-inch and 3-inch PMT combination.Each pond is filled with purified water up to 4 m above the PMT photo-cathodes.A closed recycling system is implemented to maintain water purity and achieve an attenuation length for near-ultraviolet light longer than 15 meters.
The results presented here for WCDA were obtained using the full-array configuration from 2021 March 5 to 2022 September 30.For each PMT, a hit is formed with the threshold of 1/3 photoelectron (PE) for an 8-inch PMT, and 1 PE for a 20-inch MCP-PMT.A trigger algorithm was implemented to record CR air showers by requiring at least 30 PMTs fired among 12 × 12 PMT arrays simultaneously within a window of 250 ns.To ensure a reliable data sample, data quality checks were performed based on the DAQ data status and reconstruction quality.The total effective live time used in the data analysis is around 508 days with a trigger rate around 35 kHz.The number of gammalike events recorded by WCDA is around 1.29×10 9 events after certain data selection and gammaray/background discrimination cuts.More details about the array and the reconstruction can be found in Aharonian et al. (2021a).
KM2A is composed of 5195 electromagnetic particle detectors (EDs) and 1188 muon detectors (MDs) distributed over an area of 1.3 km 2 .Each ED consists of 1-m 2 plastic scintillator covered by a 0.5-cm thick lead plate and a 1.5-inch photomultiplier tube (PMT).The typical detection efficiency is about 98% and time resolution is about 2 ns.The response time and signal amplitude of each ED is calibrated using offline methods (Lv et al. 2018;Aharonian et al. 2022).A trigger is generated when 20 EDs are fired within a 400 ns window.The signals recorded by EDs are used to determine the energy and arrival direction of the primary particles.Each MD consists of a cylindrical water tank, with a diameter of 6.8 m and a height of 1.2 m, and a 8-inch PMT, which is filled with pure water and buried under 2.5 m of soil.The MDs are designed to detect the muon component of showers, which is used to discriminate between gamma-ray and hadron-induced showers.The performance of KM2A for the gamma-rays with energies from 10 TeV to 1.6 PeV has been thoroughly tested using detector simulations and observations of the Crab Nebula (Aharonian et al. 2021b).The resolution is energy-and zenith-dependent.For showers with a zenith angle less than 20°, the angular resolution ranges from 0.5 • at 20 TeV to 0.2 • at 100 TeV.The energy resolution is about 24% at 20 TeV and 13% at 100 TeV.The rejection power of the hadron-induced showers is about 1000 at 20 TeV and greater than 4000 at energies above 100 TeV.
The KM2A detectors were constructed and merged into the data acquisition system (DAQ) in stages.The first half of the KM2A consisting of 2365 EDs and 578 MDs started operating in December 2019.This partial array was expanded to a 3/4 array, comprising 3978 EDs and 917 MDs, in December 2020.The whole KM2A was completed and operated starting 2021 July 19.The KM2A data collected from 2019 December 27 to 2022 July 31 were used for the analysis in this work.The total live time is 933 days, corresponding to 730 days of complete KM2A detector.The pipeline of KM2A data analysis presented in Aharonian et al. (2021b) is directly adopted in this analysis with the same event selection conditions.After the pipeline cuts, the number of events used in this analysis is ∼ 1.35 × 10 7 with reconstructed energies above 25 TeV.

Data Binning and Background Estimation
We use only events with zenith angles less than 50 • in this analysis, corresponding to the survey region in the declination band from −20 • to 80 • .For the WCDA data, events are divided into five groups according to the number of hits (N hit ), i.e., 100-200, 200-300, 300-500, 500-800, ≥800.For Crab-like sources, the corresponding energies roughly range from 1 TeV to 25 TeV.It should be noted that events in the same bin for a source with a harder spectrum, or at a larger declination, will tend to have a larger average energy.The PSF width of WCDA data depends on the shower size which is closely related to N hit .Therefore, it is insensitive to depend on declination or spectrum.The Crab observation provides a clear measurement of the angular resolution (denoted as ϕ 68 ), which is 0.73 • ,0.46 • ,0.37 • ,0.29 • and 0.25 • in the five bins, respectively.For the KM2A data, events are divided into five groups per decade according to reconstructed energy (E rec ), i.e, 25-40 TeV, 40-63 TeV, 63-100 TeV, 100-160 TeV, 160-250 TeV, 250-400 TeV, 400-630 TeV, 630-1000TeV, 1000-1600 TeV, >1600 TeV.The median energy and the angular resolution in each E rec bin slightly vary with the declination and spectrum of the source.The properties of each bin are shown in Table 1.
The sky maps in celestial coordinates (right ascension and declination) are divided into 0.1 • × 0.1 • pixels, and each pixel is filled with the number of the detected events according to their arrival direction.To obtain the excess of gamma-ray induced showers in each pixel, the "direct integration method" (Fleysher et al. 2004) is adopted to estimate the number of background events.This method uses events with the same direction in local coordinates but different arrival times to estimate the background.In this work, the integrated time is 10 hours and the events within the regions of the Galactic plane (|b| < 10 • ) and VHE gamma-ray sources (with space angle less than 5 • ) are excluded to estimate the background which is dominated by the residual CRs.The isotropic diffuse gamma rays and electrons may contribute slightly to the background.
The diffuse Galactic gamma-ray emission (GDE) resulting from the interaction of CRs with the interstellar medium (ISM) and background photons is an essential component of the gamma-ray sky.In Galactic plane, the VHE and UHE GDE has already been clearly detected (Abramowski et al. 2014;Amenomori et al. 2021;Cao et al. 2023).Therefore, the GDE is also an essential background to take into account for detecting and characterizing gamma-ray sources.In VHE and UHE band, the GDE have been interpreted to be mainly pion decay photons generated from hadronic interactions of CR with ISM.Approximately 99% of the ISM mass is gas.Modeling the VHE and UHE GDE requires knowledge of CR intensities and spectra, along with the distributions of gas, throughout the Galaxy.However, due to the variation of gas density and CR density across the Galaxy, it is not possible to completely disentangle the GDE.To simplify the analysis, we ignore the variation of CR density.The flux morphology of GDE is assumed to follow the spatial distribution of gas, which coexist with dust grains.Observations have shown that the gas-to-dust ratio leads to a mass ratio of M gas /M dust ∼ 100.The dust column density can therefore provide a template for the GDE morphology.In this work, we use the GDE template derived by Planck maps of dust optical depth (Planck Collaboration et al. 2014, 2016).The spectral function is assumed to be identical in all regions.To avoid the influence of the gamma-ray sources, the spectrum of GDE is estimated using the region of Galactic latitude Note-ϕ 68 is the 68% containment radius of angular resolution.E M C γ is median energy in each bin.We consider a reference source with a broken power-law spectrum of an index of 2.5 in 1-25 TeV and an index of 3.5 above 25 TeV, at a declination of 30 • .5 • < |b| < 10 • and then extrapolating it to other Galactic plane regions (|b| < 5 • ).Note that the region of several VHE sources in the region of 5 • < |b| < 10 • is also masked for the measurement of GDE, effectively excluding the contribution of both Galactic and extragalactic sources.The GDE spectrum for this region is found to be well reproduced by the diffuse Galactic gamma-ray emission model, which takes into account the local CR spectrum and the gas column density.We add the extrapolated GDE maps from the measurements in the region of 5 • < |b| < 10 • to the background maps for subsequent analysis.

Significance Map
We can calculate the significance of excess centered at each pixel using the event and background maps as described in previous section, taking the detector responses into consideration.The WCDA and KM2A data are separately analyzed.To combine the data bins of each instrument, the 3Dlikelihood method is used.A likelihood ratio test is performed between the background-only model and the one point source model.The likelihood calculation assumes that the number of counts in each pixel is distributed according to a Poisson distribution, with the mean given by the number of background counts for the background-only model and the number of background plus expected gamma-ray counts for the source model (Stewart 2009).The test statistic (TS), defined as TS = 2 ln(L/L 0 ), where L 0 is the maximum likelihood value for the null hypothesis and L is the maximum value for the source hypothesis, is used to estimate the significance.In this work, a power-law spectrum is assumed with an index of 2.6 for WCDA data in the energy range 1−25 TeV and 3.0 for KM2A data at energies E > 25 TeV as initial conditions.This leaves only one free parameter for the likelihood calculation.According to Wilks' Theorem, the TS is distributed as χ 2 with one degree of freedom (dof), and the significance can be estimated with S = √ TS. Figure 1 shows the significance maps obtained in the energy bands 1 TeV < E < 25 TeV and E > 25 TeV in Galactic coordinates.The signals are clearly visible.However, most sources in the Galactic plane are nearby and overlapping.Hence, further analysis is needed to derive each source separately.

CONSTRUCTION OF THE CATALOG
The identification of point-like gamma-ray sources and their corresponding significance can be roughly derived from Figure 1.However, it is important to note that the significance may be overestimated due to the overlap with nearby sources.Conversely, in the search for point sources, a significant portion of the sources may actually be extended, resulting in an underestimation of their significance.To improve source detection, the significance of a given source is reassessed by coupling the fitting of localization, extension and spectrum, and new potential sources are also explored.In the first step, the WCDA and KM2A data are analyzed separately to achieve two source component catalogs, one for energies ranging from 1−25 TeV and another for energies above 25 TeV.These two source component catalogs are merged into the final source catalog following a specific procedure.The sources with clear E > 100 TeV emission will also be identified in the final catalog.

All Sky Fitting and Source Component Detection
The WCDA data and KM2A data are analyzed following the same strategy to fit the entire sky region and yield the source components with the characteristics of localization, extension and spectrum.

Determination of Seeds and ROIs
The signals shown in Figure 1 are potential point sources, denoted as "seeds".To identify these seeds in the maps, all local maxima within a 0.5 • region centered around the candidate, with a significant above 4σ, are selected.Most of the seeds are concentrated in the inner Galactic plane, while in the outer Galactic plane, the seeds are clustering in the region near Geminga.In other areas of sky region, a few individual seeds stand out.
To facilitate the fitting process, all sky data are split into different regions of interest (ROI) based on the seeds and their clustering characteristics.For the inner Galactic plane region, a sliding window 20 • × 20 • ROI is adopted along Galactic longitude with a step length of 10 • .For the Geminga region, a ROI with the radius of 12 • centered at right ascension (α 2000 ) 100.2 • and declination (δ 2000 ) 16.2 • is used.For isolated seeds, the ROIs are circular regions with a 3 • radius, and the overlapping ROIs are merged into one.

Source Component Detection Procedure
The seeds obtained in Section 3.1.1may not correspond to a real source component because the spectrum and extension characteristics are not fully considered.Moreover, the position can be influenced by nearby sources with small angular separations.Thus, an analysis pipeline based on a three-dimensional maximum likelihood algorithm is developed for further analysis.Using this analysis pipeline, the spectral information and spatial morphology and significance for individual source components can be determined simultaneously.Some technical details are described below: • The first step of the 3D maximum likelihood algorithm is to build a source model characterized by a 2D-Gaussian morphology based on the initial positions of the seeds.The source spectrum is assumed to follow a power-law shape dN/dE = N 0 (E/E 0 ) −Γ , where N 0 is the differential flux at E 0 , E 0 is the reference energy which is set to 3 TeV for WCDA data and to 50 TeV for KM2A data, respectively, and Γ is the spectral index.
• Starting with one Gaussian component, we add Gaussian components successively to the source model of each ROI.A likelihood ratio test is used to compare between two models with N and N + 1 Gaussian components.We define TS = 2 ln(L N +1 /L N ), where L N and L N +1 represents the maximum likelihood of the model with N and N + 1 source components, respectively.During this iterative process, an additional source characterized by 5 free parameters (α 2000 , δ 2000 , extension, differential flux and spectral index) is added to the model if TS > 25 (3.8σ for 5 dof).
• The iteration process is terminated once no significant residual is left in the TS map for each ROI.Here, no significant residual means that no pixel is with TS > 16 (4σ for 1 dof) in the residual TS map, which is calculated by adding the current sources to the background model and using the method similar to that in Sec.2.3.
• During the fitting process, all of the parameters are optimized simultaneously.The Spectral parameter of the GDE is always keep to fixed.
• After completing the previous steps, the WCDA and KM2A source component lists are formed, respectively, by merging all the list determined in each ROI.Due to the overlap ROIs, the same source component may appear in two or more ROIs.In such cases, the one that was closest to the center of the ROI is selected.
This analysis pipeline yields a preliminary source component list, in which most components have good estimates of parameters.However two possible biases may exist for two cases: 1) the target Gaussian component is close to the ROI edge or a bright source; 2) the iteration fitting scheme described above could converge toward a local maximum.Additional checks have been implemented for these two cases using a new ROI centered on the position of each component and re-analyzing the target source component with free parameters for all components in the ROI.If new parameters of the target source component do not significantly differ from the previous iteration fitting, the new parameters of this source component are considered reliable and are adopted.The significant difference refers to the situation where the variance of a single parameter is greater than the corresponding statistical error.If significant variations are observed, a manual fitting process is conducted.This involves adjusting the starting values of the parameters and/or the size of the ROI until the parameters of the target source component converge.This convergence is determined through a visual inspection of likelihood profiles.Once the parameters are obtained, the TS value of each source component is evaluated within each ROI.In this evaluation, the positions and extensions of the other components are fixed, while the spectral parameters are left free for optimization.

False Positive Expectation
The trial factor can be roughly estimated as f t = Ω/[2π(1 − cos 0.5 • )] ∼ 3.5 × 10 4 , where Ω is the solid angle of the LHAASO survey sky within the declination ranging from −20 • to 80 • , and 0.5 • is the minimum location searching radius for seeds.Since the TS value for the above components follows the χ 2 distribution with 5 degrees of freedom, TS = 25 corresponds to the pvalue of 6.9 × 10 −5 .Based on this, we can expect that 2.4 sources in our survey are false detections due to background fluctuations.Thus, we set a higher threshold of TS = 37 (5σ for 5 dof) for reported sources, corresponding to an expectation of 0.01 false-positive sources.For sources with both WCDA and KM2A components, the summed TS value of 50 (5.1σfor 10 dof) corresponds to the comparable false-detection to that of the single component of TS = 37.Hence, we also report the sources of which two components are with 25 < TS < 37. Furthermore, due to the mismodeling GDE background, the sources in the gas-rich region with a low significance level might still be potential false detection.Therefore, further study is required to determine whether or not they are real sources.

Localization, Extension, and Energy Spectrum
The localization of each source component, represented by α 2000 and δ 2000 , is given by the procedure in Sec 3.1.2.The statistical errors of the localization (σ α 2000 ,stat , σ δ 2000 ,stat ) are also estimated by analyzing the shape of the likelihood function with the HESSE tool in the Minuit2 package.It is worth noting that the localization of a faint or soft source component is more sensitive to the influence of diffuse emission and nearby sources.As a result, we may find asymmetric position error circles for these components, where σ α 2000 ,stat cos δ 2000 ̸ = σ δ 2000 ,stat .To account for this, we conservatively consider the bigger error as the one-dimensional 1σ position error, i.e., x stat = max[σ δ 2000 ,stat , σ α 2000 ,stat cos δ 2000 ].In this work, we report the position uncertainties at the 95% confidence level as σ p,95,stat = 2.45×x stat , where the 2.45 is the factor −2 ln(1 − 0.95) by which the 95% confidence level position uncertainty is related to x stat for a two-dimensional axisymmetric Gaussian.Figure 2 illustrates the position uncertainties as a function of the TS values.The relatively large dispersion seen at a given TS is due to variations in local conditions , the source extension and the source spectrum.Considering the systematic error, the error radius (95% confidence) of the LHAASO sources ranges from ∼ 0.04 • to ∼ 0.8 • .The extension size of each source component is described by a two-dimensional Gaussian σ (namely r 39 in this work) corresponding to 39% of the source flux.A statistical error for the measured extension (σ r 39 ,stat ) at the 1σ level can also be determined.The measurement of r 39 is invalid for the point-like components unresolved by LHAASO.Thus, for all source components, the extension Test Statistic (TS ext ) is performed.Here TS ext = 2 ln(L ext /L ps ), where L ps and L ext are the maximum-likelihood values assuming the source spatial model is a point-like and 2D-Gaussian model, respectively.When the TS ext is lower than 9 corresponding to a significance smaller than 3σ level, we change the 2D-Gaussian spatial model to a point-like model and re-optimize the position and spectral parameters.The upper limits on the extension (r 39,ul,stat ) at the 95% confidence level for these point-like components are derived from the shape of the likelihood function.Figure 3 shows the extension size of each source component.Roughly 1/3 are unresolved by LHAASO and modeled by a point-like morphology, and roughly 2/3 are modeled by a 2D-Gaussian morphology with an extension size ranging from ∼ 0.2 • to > 2 • .The spectral energy distribution (SED) of each source component is modeled by a simple power law (PL).Figure 4 shows the distribution of the SED parameters of the source components.The differential flux has a mean of ∼ 3 × 10 −12 TeV cm −2 s −1 for WCDA components at 3 TeV and ∼ 4 × 10 −13 TeV cm −2 s −1 for KM2A components at 50 TeV.The photon index distributions are obviously different, with an average of ∼ 2.5 for WCDA components and ∼ 3.5 for KM2A components, implying that a single power law shape cannot describe the overall SED for the majority of 1LHAASO catalog sources.
We compare the flux measurements with the sensitivity of source derived by our used data set.The flux sensitivity is defined as the flux normalization required to have 50% probability of detecting a source at 5σ level.As shown in Figure 5, the sensitivity dependents on the declination, spectral index and source size.Sources passing directly overhead, with a declination of 30 • for LHAASO (20 • for HAWC), exhibit the highest sensitivity.The sensitivity for a source with a softer spectrum or larger size, or at a larger declination, will tend to decrease.Thus, about 10 WCDA sources, which have fluxes above the 3HWC point source sensitivity, have not been included into the 3HWC catalog due to their large size and/or soft spectrum.In addition, it is possible that the source searching strategy employed in the 3HWC catalog results in the failure to detect several of these sources.
Figure 6 shows the 1LHAASO source fluxes related to source sizes or source photon spectral index.In the flux-size panel, we illustrate the approximate flux sensitivity limit of the 1LHAASO source as a function of source size.It can be observed that the sensitivity worsens as the source size increases, as mentioned earlier.For extended sources with sizes larger than 1 • , we may find a possible lack of sources with flux F 3 > 4 × 10 −13 TeV −1 cm −2 s −1 or F 50 > 4 × 10 −16 TeV −1 cm −2 s −1 .In the index-flux panel, the photon index in the WCDA band appears to be uncorrelated with the flux F 3 .However, there might be a broken power-law correlation between the photon index in the KM2A  2).Left: the integral sensitivity is shown at 3 TeV for three hypotheses: 1) point source with spectrum of E −2.0 ; 2) point source with spectrum of E −2.5 ; 3) 0.5 • gaussian source with spectrum of E −2.5 .The 3HWC sensitivity corresponds to the point-source search with the 3HWC data set (Albert et al. 2020).Right: the integral sensitivity is shown at 50 TeV for three hypotheses: 1) point source with spectrum of E −3.0 ; 2) point source with spectrum of E −3.5 ; 3) 0.5 • gaussian source with spectrum of E −3.5 .
band and the logarithm of F 50 , with a break occurring at a flux of F 50 ∼ 4 × 10 −17 TeV −1 cm −2 s −1 .It is not yet clear whether this phenomenon is caused by an instrument selection effect or a physical relation.

Source Detection Above 100 TeV
An exciting characteristic of LHAASO is that it can extend the observation to the UHE regime owing to the large detector area.The UHE sources are crucial to identify the emission mechanism and explore the maximum acceleration energy within the sources, and they are also very important

Merging of the Source Components
A WCDA component and a KM2A component are considered to be the same source according to the agreement in position, i.e., r 2 d < (σ p,95,W ) 2 + (σ p,95,K ) 2 + δ 2 p , where r d is the positional offset between the WCDA component and the KM2A component, and σ p,95,W and σ p,95,K are the position uncertainties at the 95% confidence level for the WCDA and the KM2A components, respectively.δ p is the position-offset correction of the components, possibly due to the physical position variation or the changes in morphology at different energy bands.Additionally, incomplete descriptions of the GDE and/or of nearby sources in the analysis can also contribute to the positional offset.It is difficult to determine the value of the δ p for each source.We expect that δ p is related to the source extension.As a simple criterion, we adopt a loose condition with δ p = (r 39,W + r 39,K ), where r 39,W and r 39,K are the extensions of the WCDA and KM2A components, respectively.We merge the closest pairs of the WCDA and KM2A components that satisfy this condition.Ultimately, 54 pairs of components are merged.

Known TeV Source Association
We have conducted a preliminary association of LHAASO sources with known TeV sources based on their position.A searching radius is defined by r sr = σ 2 p,95 + r 2 39 + (0.3) 2 , where σ p,95 and r 39 are the position error and extension size of the source component.For the LHAASO sources with KM2A and WCDA components, the position and extension of the component with higher significance is used.For the point-like source components, r 39 is set to the extension upper limit.The average value of the position errors for the known TeV sources detected by EAS arrays is about 0.3 • .While for the sources detected by IACTs, the average position error is around 0.05 • .Roughly, we used the factor of 0.3 degree for the searching radius.If r sr > 1 • , a maximum searching radius of 1 • is applied.
The vast majority of the known TeV sources from ground based observatories are included in the TeV online catalog, i.e., TeVCat1 (Wakely & Horan 2008).It is frequently updated, containing 252 TeV source entries to date.When associating LHAASO sources with known TeV sources, we used the canonical name and position reported in TeVCat.We exclude the TeV pulsars for the associations, since the pulsed TeV emission is expected to be null in our source searching procedure.If no TeV source in TeVCat is found for association, we further compare 1LHAASO sources with 3HWC sources, which include some sources not yet included in TeVCat so far.We list the closest TeV counterparts within the searching radius for each LHAASO source, as shown in table 2.

Systematic Uncertainties and Caveats
Based on the measurement of the standard candle, the Crab Nebula, with a position of α 2000 = 83.633• and δ 2000 = 22.015 • for the systematic error estimation, a systematic pointing error can be evaluated as 0.03 • for KM2A data.For the WCDA data, an investigation of the pointing systematic error has been conducted using three point sources: Crab Nebula, Mrk 421, and Mrk 501.The maximum error is from Mrk 501 at about ∼ 0.04 • , which is taken as the pointing error of the WCDA array in this work.It is important to note that no clear zenith angle dependence of the systematic error is found according to the observation of the Crab Nebula using events at different zenith angles, ranging over 0 • -15 • , 15 • -30 • , and 30 • -50 • .Furthermore, no systematic location bias related to declination has been identified based on the observation of the CR Moon shadow at different declination bands.
The systematic uncertainty of the size of 1LHAASO sources could be contributed by the uncertainties of the Point Spread Function (PSF).For KM2A data, we obtained a systematic bias for the ϕ 68 at the order of ∼ 0.08 • by comparing the Crab measurement with our simulated data.For WCDA data, we compared the observed events profile among Crab, Mrk 421, and Mrk 501.The uncertainties can also yield a systematic bias of ∼ 0.05 • .Thus, we can conservatively estimate the systematic error to be at the order of σ r 39 ,sys ∼ 0.05 • for WCDA and σ r 39 ,sys ∼ 0.08 • for KM2A component.
The systematic errors affecting the spectrum have been investigated in Aharonian et al. (2021a) and in Aharonian et al. (2021b).The main systematic error is contributed by the atmospheric model in the Monte Carlo simulations.The total systematic uncertainty is estimated to be 7% on the flux and 0.02 on the spectral index for KM2A SED measurement.In the case of WCDA data, the overall systematic uncertainty can be as large as +8% −5% on the flux, which is estimated by the same method in Aharonian et al. (2021a).The power-law spectral shape can adequately describe the WCDA components, while it is not suitable for about 1/3 KM2A components which shown an evident curved shape at energies above 25 TeV.A detailed spectrum study needs to be carried out in the future.
Mismodeling of GDE could affect the fitted locations, extensions, and SEDs of several source components, especially for those with lower fluxes and larger sizes.A rough assessment of the impact of the GDE was performed by excluding it from the background maps for the source components in the Galactic plane region.As a result, 11 KM2A source components and 10 WCDA source components exhibited changes (in terms of location, flux, index, or extension) exceeding three times the statistical errors.Source components that were significantly affected by the GDE were tentatively labeled as such.The GDE test is not thorough because we cannot determine all the sources affected by the GDE mismodeling.As a conservative approach, source components with an extended size greater than 2 • were excluded due to the growing impact of the GDE with increasing component size.Further understanding of very extended sources or obvious GDE-impacted sources requires deeper study of the GDE model, which is beyond the scope of this work.
As shown in Figure 7, at 95% confidence level, approximately 40% of the merged sources exhibit a noticeable bias in position or differences in extension.The shift could be an indication of an energydependent morphology of these sources.However, it should be noted that due to the poor angular resolution of LHAASO, the offset in position and extension could also be influenced by emission from nearby unresolved sources or mismodeling of GDE.There is a possibility of a false merged source, which could be a combination of WCDA and KM2A components that are physically unrelated.To further examine this, we have roughly conducted individual investigations on the surface brightness distribution and spectrum of each merged source.During this process, conflicts were found in 11 merged sources, which are tentatively labeled as dubious mergers.More comprehensive studies are needed to confirm the physical association between the two components of the dubious merged sources.

RESULTS
Following the above procedure, the source catalog of LHAASO has been constructed.Overall, 90 sources with extension < 2 • are found over the whole LHAASO survey sky.Among them, 65 sources exhibit extended morphology with a confidence level greater than 3σ.A total of 54 sources have been simultaneously detected by both WCDA and KM2A.Among all the sources, 43 UHE sources have been detected at > 4σ confidence level when E > 100TeV.
Table 2 presents a comprehensive list of all LHAASO sources obtained through the above procedure, ordered by α 20002 .To illustrate these catalog sources, a detailed view of 82 sources with Galactic latitude |b| < 12    The search for and identification of galactic sources capable of accelerating CRs up to PeV energies, known as PeVatrons, is of great importance.UHE gamma rays serve as a crucial and highly promising tool for achieving this target.With its unprecedented sensitivity at UHE, LHAASO represents the best instrument to survey the PeVatrons.In fact, using approximately 300 days of data collected by the half KM2A, LHAASO has already detected 12 UHE sources with significance above 7σ and maximum energies reaching up to 1.4 PeV (Cao et al. 2021a), providing crucial candidates for hadronic PeVatron exploration.In this paper, 43 sources with significance above 4σ at energy beyond 100 TeV are listed in Table 2. Since these sources are well detected at E > 25 TeV with much higher significance, they are significant enough to be identified as UHE sources.Among the 43 sources, 22 sources are detected with significance above 7σ (TS 100 ≈ 54) at E > 100 TeV, improving by a factor of 2 the previous number of sources in Cao et al. (2021a).It is worth noting that 57% (43 out of 75) E > 25 TeV sources are UHE sources.According to Figure 9, most of these UHE sources have higher significance or harder spectral index than the other E > 25 TeV sources, which indicates that the remaining E > 25 TeV sources may also be detected as UHE sources by LHAASO in the future with further accumulation of data.This provides important evidence for the exciting fact that the Milky Way is full of UHE sources and full of PeV particle accelerators.Further deep analysis focusing on these sources one by one or type by type to identify the emission mechanism is a straightforward task to explore and identify the hadronic PeVatrons.According to Table 2, 76% (33 out of 43) UHE sources are detected at energies 1−25 TeV by WCDA. Figure 9 shows the spectral index versus significance for these sources at energies 1−25 TeV.For comparison, other sources detected at energies 1−25 TeV are also presented in the figure, excluding the extragalactic sources listed in Table 2.The UHE sources have harder spectral index or higher significance than the others in the energy band 1-25 TeV.Up to now, more than 250 VHE sources have been detected, and more than 100 sources are within the Galaxy with a significant fraction being located in the southern sky, which is out the FOV of LHAASO.These Galactic sources would be important UHE candidates and may be revealed with UHE emission in the future observations.
It is worth noting that 10 out of 43 UHE sources were not detected by WCDA in 1-25 TeV band.These sources exhibit dominant gamma-ray emissions at energies around a few tens of TeV or E > 100 TeV.Among them, sources like 1LHAASO J0216+4237u possibly belong to a new population of gamma-ray sources.These would demonstrate the distinctive importance of the UHE window at higher energy, which could explore new phenomena and new extreme celestial bodies of the Universe.With more accumulation of data, LHAASO will discover more sources with such features in the future.These sources should be also important candidates for PeVatron exploration and identification.

TeV Sources Discovered by LHAASO
As shown in Table 2, there are 25 sources without any known TeV source association.Additionally, 7 sources have conflicting extensions (with a difference of more than 0.5 • ) compared to the known TeV sources within the searching region.Therefore, tentatively, these 32 sources are announced as new TeV sources detected by LHAASO in this catalog.To determine the physical nature of these sources it is often necessary to analyze their spectral and morphological features, as well as the information about multi-wavelength associations.In this section, we conduct a preliminary analysis of the associations for each new TeV source.The majority of these newly discovered LHAASO sources are expected to be Galactic sources, due to their extended properties or their KM2A component detection.However, several point-like sources detected only by WCDA could potentially have an extra-galactic origin.It is known that there are a number of types of Galactic gamma-ray emitters, such as SNRs, PSRs and their PWNe, massive star clusters (MSCs), star-forming regions(SFRs), superbubbles, binaries, etc.The vast majority of detected Galactic TeV sources are likely to be associated with SNRs and PWNe (H.E. S. S. Collaboration 2018).To search for SNR and PWN associations, we take into account the most complete catalog of SNRs and PWNe to date, SNRcat3 (Ferrand & Safi-Harb 2012).Although we do not expect pulsed TeV emission from a pulsar to be detected by LHAASO, the pulsar is a probe to illustrate the characteristic properties of an association SNR and/or PWN, and is also an indicator of a possible unseen SNR and/or PWN.For the pulsar associations, we used the web-based Australia Telescope National Facility Pulsar catalog (ATNF pulsar catalog4 ).There is less association between low spin-down luminosity pulsars and known TeV sources, hence we apply a cut of Ė > 10 34 erg s −1 .Additionally, we search for associated high-energy (HE) gamma-ray sources in the 4FGL catalog (Abdollahi et al. 2020), which covers the sources at energies ranging from 0.1 GeV to >500 GeV.Note that in this section we claim an association or a counterpart based on a position less than 0.5 • away.For the likely extragalactic sources, the AGN catalog is adopted to search for associations within the position error.The association results are listed in Table 3.
Seven of the new TeV sources do not have any associations.They are tentatively classified as "dark" sources in this work.We list and discuss these "dark" sources briefly here: • 1LHAASO J0007+5659u is a point-like source detected only by KM2A.It is an UHE source with a significance of TS 100 = 43.5 at energies above 100 TeV, implying a Galactic origin.
• 1LHAASO J0206+4302u and 1LHAASO J0212+4254u are the point-like sources only detected by KM2A.Due to both of them being UHE sources, we suggest a Galactic origin, although these two sources are obviously far from the Galactic plane with Galactic latitude b ∼ −17 • .As shown in Figure 13, these two sources and 1LHAASO J0216+4237u are close by and likely to be connected by bridges.Considering that they have a similar spectral shape (with the index of ≈ 2.5), we favor a physical association among these three 1LHAASO sources.
• 1LHAASO J1937+2128 is an extended source with the size of ≈ 1.3 • with the WCDA and KM2A components.Although 3HWC J1935+213 (0.36 • away) and 3HWC J1936+223 (0.85 • away) are found in point searches and are within the 1 • region around the center position of this 1LHAASO source, we report this source as a new detection due to its larger extension size.
• 1LHAASO J1959+1129u is a point-like TeV source detected above 25 TeV with TS = 90.2, and also detected as an UHE source with TS 100 = 59.3.We cannot find any pulsar, SNR or PWN counterparts within a 0.5 • region around its position.Interestingly, it is 0.23 • from the low mass X-ray binary (LXB) 4U 1957+11 which is a black hole candidate with the simplest and cleanest soft, disk-dominated spectra (Nowak et al. 2012) in the X-ray regime.The TeV emission from 11LHAASO J1959+1129u is possibly associated with LXB 4U 1957+11.A dedicated study of 1LHAASO J1959+1129u from our LHAASO group is ongoing.
• 1LHAASO J2200+5643u is an extended TeV source with the size of ∼ 0.5 • detected by WCDA and KM2A simultaneously.It is an UHE source with the significance of TS 100 = 38.The spectral index is ≈ 1.77 at energies ranging from ∼ 10 TeV to 30 TeV as measured by WCDA and is ≈ 3.44 at energies above 25 TeV as detected by KM2A, implying a peak or a break energy at a few tens of TeV.
• 1LHAASO J2229+5927u is a large extended source with the size ∼ 1.8 • detected by WCDA and by KM2A.Within the 1 • region centered at the position of the WCDA component, we cannot find any counterpart.The famous SNR G106.3+2.7, which has been extensively studied as a candidate PeVatron, is located at the edge of the extended region of 1LHAASO J2229+5927u.It is possible that this 1LHAASO source is the product of CRs escaping from SNR G106.3+2.7.More detailed study for this 1LHAASO source is ongoing by our collaboration.
Eight of the new TeV sources have GeV counterparts only.The details of these sources are as follows: • 1LHAASO 0056+6346u is an extended source (r 39 ≈ 0.3 • ) detected by WCDA and KM2A detectors simultaneously.It is also a UHE source with TS 100 = 94.The spectral shape is obviously curved with the index ≈ 2.35 by WCDA and ≈ 3.33 by KM2A.The unidentified point-like GeV source 4FGL J0057.9+6326 is 0.38 • from this 1LHAASO source.
• 1LHAASO J0500+4454 is an extended source with a size of ∼ 0.4 • , only detected by WCDA in this work.The extension characteristic supports a Galactic origin.An unidentified point-like GeV source 4FGL J0501.7+4459 is 0.32 • from this 1LHAASO source.
• 1LHAASO J1858+0330 is an extended source with a size of ≈ 0.5 • detected by WCDA with a significance of TS = 114 and by KM2A with a significance of TS = 299.It is located in the Galactic plane (l = 36.8• and b = 0.1 • ), towards which a complex gas distribution exists.The closest GeV source is 4FGL J1857.9+0313c which is identified as the a blazar candidate of uncertain type by the Fermi-LAT group.We exclude that the 1LHAASO source is associated with this extragalactic GeV source due to the extension characteristic.Another unidentified point-like GeV source 4FGL J1858.0+0354 is 0.41 • from this 1LHAASO source, and an association cannot be ruled out.
• 1LHAASO J2027+3657 is an extended source with r 39 ≈ 0.23 • only detected by KM2A.It is resolved from the previously published LHAASO source LHASSO J2018+3651, which is one of the brightest sources observed by LHAASO.The point-like GeV source 4FGL J2026.5+3718c is located within a 0.5 • region around this 1LHAASO source.
A weak WCDA component with TS ≈ 20 is found but is not included in this catalog.The unidentified point-like GeV source 4FGL J02049.3+4440c is 0.32 • from this 1LHAASO source.
• 1HAASO J0249+6022 is an extended source with r 39 ≈ 0.71 • for the WCDA component and with r 39 ≈ 0.38 • for the KM2A component.We can find that the KM2A component overlaps the extended region of the WCDA component, with a positional offset of ∼ 0.45 • .The WCDA component at energies 1−20 TeV could plausibly include the emission of two TeV sources, resulting in the large position offset compared to the KM2A component.In our searching radius, we just find one pulsar counterpart (PSR J0248+6021, d = 2.0 kpc Ė = 2.13 × 10 35 erg s −1 , τ c = 62.4 kyr, p 0 = 217.1 ms), 0.16 • from the position of the KM2A component.The 1LHAASO source is likely to correspond to an astrophysical system associated with PSR J0248+6021, such as a composite SNR.
• 1LHAASO J0428+5531 is an extended source detected by WCDA and KM2A detectors simultaneously.The WCDA component is with r 39 ≈ 1.18 • .We can find that the shell-type SNR G150.3+04.5, with a size of 3.0 • × 2.5 • (Gao & Han 2014), is 0.26 • from the center of the WCDA component.In the GeV band, an extended source 4FGL J0426.5+5522emodeled as a 1.515 • disk is also associated with SNR G150.3+04.5 (Ajello et al. 2017).The WCDA component of 1LHAASO J0428+5531 is possibly of SNR origin due to the comparable location and extension size in the multi-wavelength observation.The extended KM2A component (r 39 ≈ 0.23 • ) is ≈ 0.97 • from the central position of the WCDA component and ≈ 0.1 • away from the radio west shell of SNR G150.3+04.5 (with a size of 64.1 ′ × 18.8 ′ , also named SNR G150.8+03.8 in Gerbrandt et al. 2014).A point-like GeV source 4FGL J0426.5+5434, which is possibly a pulsar candidate due to the curved SED and the cutoff energy at a few GeV, is 0.06 • from the position of KM2A component.Whether the KM2A component of 1LHAASO J0428+5531 is associated with the SNR G150.8+03.8needs more study due to the obvious offset of the position and the extension size.
• 1LHAASO J0534+3533 is a point-like source detected by WCDA and KM2A in this work, which is near the centroid position of the shell-type SNR G172.8+01.5, 0.3 • away.However, the radio shell size of SNR G172.8+01.5 is estimated as 4.4 • × 3.4 • , which is much larger than the TeV emission size of 1LHAASO J0534+3533.If the radio shell is physically associated with 1LHAASO J0534+3533, the TeV emission is likely to be a product of the central pulsar wind nebula.
• 1LHAASO J1814-1636u is an extended UHE source with the size of ≈ 0.68 • , only detected by KM2A.Only two SNRs, i.e., G14.1-00.1 (0.43 • away) and SNR G014.3+00.1 (0.31 • away) are found in our searching catalog.However, both SNR candidates are with < 0.2 • radio shell size, which is inconsistent with that of gamma-ray emission.An unidentified GeV source 4FGL J1816.2-1654c is found at 0.4 • away from the position of the TeV emission.
• Source 1LHAASO J1831-1028 only has an extended KM2A component and with the size of ≈ 0.94 • .In our association procedure, we just find three SNR candidates, i.e., G021.0-00.4(0.3 • away), G021.5-00.1 (0.36 • away) and G021.8-00.6 (0.49 • away), all of which are with a radio size of < 0.2 • .Due to the larger conflict between the size of TeV emission and that of radio shell, we disfavor that the same population of electrons produces the TeV emission and the radio emission.
• 1LHAASO J1852+0050u is an extended source with a ∼ 0.64 In addition, the SNRcat source G033.6+00.1, which is a small extended radio source with the size of 10 ′ , is suggested to be associated with the GeV source.
• 1LHAASO J1954+3253 is an extended TeV source (r 39 ≈ 0.17 • ), only detected by WCDA with the significance of TS = 144.SNR G069.0+02.7 (also named CTB 80) with the size ∼ 80 ′ overlaps this 1LHAASO source.SNR CTB 80 is an old SNR with a kinematic distance of ∼ 1.5 kpc.At the center of SNR CTA 80, the pulsar PSR B1951+32 ( Ė = 3.74 × 10 36 erg s −1 ,τ c = 107.0kyr, d = 3.00 kpc) is generally regarded as the compact object associated with this old SNR.A X-ray PWN around the pulsar PSR B1951+32 has been identified by ROSAT, which is shown as a 5 arcmin extended nebula.The pulsar and the SNR/PWN are also detected by F ermi-LAT, which are 4FGL J1952.9+3252 and 4FGL J1955.1+3321,respectively.We note that the position of TeV emission is ∼ 0.33 • from that of the pulsar and its X-ray PWN, for which more study is needed to confirm the physical association between the 1LHAASO source and CTB 80.
• 1LHAASO J2002+3244u is a point-like source with TS=74.0 as detected by WCDA and with TS = 43.6 as detected by KM2A.It is also an UHE source with TS 100 = 28.1.4FGL J2002.3+3246 is spatially coincident with this source, which is identified as a potential association with a SNR or PWN by the F ermi-LAT group.A shell type SNR G069.7+01.0 is possibly associated with 1LHAASO 2002+3244u (0.04 • away).We favor that 1LHAASO 2002+3244u has an SNR origin because the position and size of TeV emission agree with the radio shell of SNR G069.7+01.0.
• 1LHAASO J2028+3352 is a large extended source with r 39 ≈ 1.6 • , only detected by KM2A.0.36 • from their centroid position, a middle-aged pulsar PSR J2028+3332 ( Ė = 3.48 × 10 34 erg s −1 ,τ c = 576.0kyr) is found, implying a possible TeV halo identification.Two GeV sources are found within the 0.5 • region around the 1LHAASO source, of which one is the pulsar PSR J2028+3332 and the other is an unidentified point-like source.
• 1LHAASO J2238+5900 is an extended source with the size of ≈ 0.51 • as detected by WCDA with TS = 110.2and with the size of ≈ 0.44 • detected by KM2A with TS= 361.Just one pulsar PSR J2238+5903 ( Ė = 8.99 × 10 35 erg s −1 , d = 2.83 kpc, τ c = 26.6 kyr) is found within 0.5 • region of this source, at 0.07 • from the centroid position of the KM2A component.The young pulsar and the size of TeV emission shrinking with increasing energy support that 1LHAASO J2238+5900 has a PWN origin.
One possible new extragalactic source is as follows: • 1LHAASO J1219+2915 is a point-like source, only detected by WCDA with a significance of TS = 49.2, and located at high Galactic latitude (b ∼ 82.5 • ).We cannot find any pulsar or SNR/PWN counterpart associated with this source.It is a likely extragalactic source due to the high Galactic latitude and null detection by KM2A.The closest AGN counterpart is the LINER-type AGN NGC 4278, 0.05 • from the 1LHAASO source.It is the most possible association even though we can not firmly identify this 1LHAASO source at present.

Pulsar Wind Nebulae or TeV Halos in the 1LHAASO Catalog
Pulsars left behind from supernova explosions are rapidly spinning neutron stars with extremely strong magnetic fields.The pulsed gamma-ray emission dominates at GeV energies, and only a few cases can extend to VHE.At VHE, the steady emission is from a PWN produced by the termination of the ultra-relativistic wind.PWNe constitute one of the largest VHE source populations within the Galaxy.The Crab nebula PWN has also been identified as an electron PeVatron (LHAASO Collaboration 2021).The diffusion of escaping particles from a PWN would leads to extended VHE gamma-ray emission, which has been extensive discussed as a "TeV halo" after the discoveries of several cases (Abeysekara et al. 2017;Aharonian et al. 2021c).TeV halos are thought to form around middle-aged pulsars with ages of at least several tens of thousands of years.PWNe or TeV halos should also contribute a significant fraction of the 1LHAASO sources, and an effective method to check for this is to search for spatial coincidence with the known pulsars.
To search for pulsars associated with 1LHAASO sources, the ATNF pulsar catalog is adopted.For the hunting of associations in astronomy, an important work is to estimate the possibility of spatially coincidence by accident.A similar method as that used in Mattox et al. (1997) is adopted here.According to a previous empirical result at VHE, the identified PWN or halo type VHE sources are always associated with pulsars with high spindown power.Therefore, for each pulsar within 0.5 • of a 1LHAASO source, the chance probability P c is estimated according to the space angle r and pulsar spindown power Ė using Where r 0 is the characteristic angle between confusing sources, and where ρ( Ė) is the number density of pulsars with Ė not lower than that of the candidate pulsar.With this searching, 65 1LHAASO sources are found with at least one pulsar nearby within 0.5 • .To decrease the fake association by accident, the pulsars with chance probability higher than 1% are excluded.After this filter, 35 1LHAASO sources are found with one associated pulsar each and 2 1LHAASO sources, 1LHAASO J0359+5406 and 1LHAASO J1929+1846, are found with two associated pulsars each.For the source with two associations, the associated pulsar with a lower chance probability is listed.There are also two pairs of 1LHAASO sources, 1LHAASO J1848-0001u vs. 1LHAASO J1850-0004 and 1LHAASO J2020+3638 vs. 1LHAASO J2020+3649u, associated with the same pulsar.Again, the 1LHAASO source with a lower chance probability is listed.Finally, 35 associated pulsars are derived for the 1LHAASO sources.Detailed information about these associations is listed in Table 4.
Figure 10 shows the spindown power versus age of the associated pulsars.For comparison, all the pulsars of the ATNF catalog within LHAASO FOV are also shown in the figure.As expected, the TeV-gamma associated pulsars are among that with the most energetic spindown power ( Ė > 10 34 erg s −1 ).Among these associated pulsars, the ages of 24 pulsars are less than 100k years.The corresponding 1LHAASO sources have a high possibility to be PWNe.The ages of 11 pulsars are older than 100k years, and this marks the corresponding 1LHAASO sources to possibly be PWNe/TeV halos.It is worth noting that some of them have already been identified as PWNe or PWN/TeV halos in the TeVCat, as marker in the Table 4.An exciting result is that one pulsar, i.e., PSR J0218+4232, is a millisecond pulsar, which has been expected to produce VHE emission but there has been a lack of observation evidence.The spatial coincidence between this millisecond pulsar and 1LHAASO J0216+4237u is notable, with a confidence level of 2.9σ.This finding favor the existence of VHE emission around millisecond pulsars, although a non-physical association cannot be definitively ruled out.According to Figure 10, most of the energetic pulsars with Ė higher than 10 36 erg s −1 .within the FOV of LHAASO are associated with 1LHAASO sources.These shows that the PWNe of energetic pulsars are promising as effective to VHE gamma-ray emitters.It is worth noting that no VHE or UHE emission was found from a handful of energetic pulsars with Ė higher than 10 37 erg s −1 .according to Figure 10.This shows that emission of PWNs may also be affected by other parameters of the pulsars or the surrounding environment, besides the spindown power.Among the 35 1LHAASO sources with pulsar associations, twenty-two are labeled as UHE sources.This result suggests that PWNe also contribute a significant portion of the UHE sources and acceleration of electrons close to or up to 1 PeV is common for energetic PWNe.In addition to the Crab Nebula, more electron PeVatrons may be revealed by LHAASO in the future.

SUMMARY
The 1LHAASO catalog, with 90 VHE gamma-ray sources, is the first comprehensive search conducted using data from the LHAASO observatory.The catalog utilizes 508 days of data from the full WCDA and 933 days of data from the full KM2A and its partial array.With 32 sources discovered, this survey is the most sensitive and roughly unbiased investigation of large sky regions in the VHE band, covering a declination range from −20 • to 80 • .Seven new sources do not have counterparts in the GeV gamma-ray, pulsar, and SNR catalogs, and are referred to as "dark sources".Among the 90 1LHAASO sources, sixty-nine are detected with emissions in the energy range 1-25 TeV, while seventy-five are detected at emission energy of E > 25 TeV, both with a significance above 5σ.The catalog provides information on the extension and spectrum of each source, with 65 sources exhibiting clear extension.To be conservative and avoid confusion with diffuse gamma rays from the Galactic plane, only sources with extensions less than 2 • are included in the 1LHAASO catalog.Additionally, 43 sources exhibit clear UHE emission, with a significance above 4σ at E > 100 TeV.This significantly expands the current number of known UHE sources, by a factor of 3, and suggests that most Galactic VHE sources are potential UHE sources, indicating a prevalence of PeV particle accelerators in the Milky Way.Ten out of the forty-three UHE sources are not detected at the 1-25 TeV range, possibly representing a new class of gamma-ray sources dominated by emission above tens of TeV.Thirty-five 1LHAASO sources are associated with energetic pulsars, with a chance coincident probability of less than 1%.Apart from sources already identified as PWNe or TeV halos, the remaining sources associated with energetic ( Ė > 10 36 erg s −1 ) pulsars are likely PWNs or TeV halos.The majority of energetic pulsars within the LHAASO field of view are associated with 1LHAASO sources, and a significant fraction of these sources are also UHE sources.This suggests that PWNe associated with energetic pulsars are effective emitters of VHE and UHE radiation and have a general capability to accelerate electrons close to or up to PeV energies.
With a large FOV and full duty cycle, the full LHAASO detector configuration has been monitoring the sky since July 2021 and will continue operations for the next 20 years.With the the accumulation of data, the sensitivity in both the VHE and UHE bands will improve.Further optimization of data analysis and refining detector response simulation are ongoing, which may also expand its capacity.Deep and multi-wavelength analysis focusing on the 1LHAASO sources one by one is also ongoing.Therefore, LHAASO is expected to unveil more new discoveries and provide a deeper understanding of the VHE and UHE universe in the near future.

ACKNOWLEDGMENTS
We would like to thank all staff members who work at the LHAASO site above 4400 meters above sea level year-round to maintain the detector and keep the water recycling system, electricity power supply and other components of the experiment operating smoothly.We are grateful to Chengdu Management Committee of Tianfu New Area for the constant financial support for research with LHAASO data.This research work is supported by the following grants: The National Key R&D program of China No.2018YFA0404201, No.2018YFA0404202, No.2018YFA0404203, No.2018YFA0404204, National Natural Science Foundation of China No.12022502, No.U1831208, No.12205314, No.12105301, No.12261160362, No.12105294, No.U1931201, No.12005246, No.12173039        Note-The first column lists the LHAASO catalog name.The name designation is 1LHAASO JHHMM+DDMM, where the 1 refers to this being the first LHAASO catalog, JHHMM+DDMM is according to the source location.For a source observed by both WCDA and KM2A, the source name corresponds to the position of that with higher significance level.If the source is an UHE source with TS 100 > 20, we add "u" at the end of the name, i.e., 1LHAASO JHHMM+DDMMu.For dubiously merged sources, we label the names with "*".
Note-The 2nd column lists the component name, represented by the detector name.If the component suffers a significant GDE impact in our test, we label the component with "*".
Note-The 3rd − 5th columns list the position parameters for the source components.α 2000 and δ 2000 are the right ascension and declination at J2000.0 epoch.σ p,95,stat is the 95% statistical uncertainty of the source component position.The units are degrees.
Note-The 6th column lists the extension of the source components.r 39 is the 39% containment radius of the 2D-Gaussian model.For the point-like source components, the 95% statistical upper limits (i.e., < r 39,ul,stat ) are shown.
Note-The 7th column lists the TS values of the source components.It is distributed as χ 2 with 5 dof.
Note-The 8th − 9th columns list the parameters of the power-law SED of the source components.The power-law shape is defined by dN/dE = N 0 (E/E 0 ) −Γ , where N 0 is the differential flux, E 0 is the reference energy which is 3 TeV for the WCDA component and 50 TeV for the KM2A component, and Γ is the photon spectral index.N 0 is in a units of 10 −13 cm −2 s −1 TeV −1 and 10 −16 cm −2 s −1 TeV −1 for the WCDA the KM2A components, respectively.For a source with single component detection, we give the 95% statistical upper limits of N 0 of the non-detected component with the same position, extension and photon index.
Note-The 10th column lists the TS values of the UHE sources at energies E > 100 TeV.
Note-The 11th column lists the closest known TeV source counterpart within the searching region (as described in Sec.3.5).The angular separation between the TeV counterpart and the LHAASO source is shown in parentheses.

Figure 1 .
Figure 1.Significance maps of the region monitored by LHAASO.A point test source with a spectral index of 2.6 for WCDA data and 3.0 for KM2A data is used.

Figure 2 .
Figure 2. Source component location uncertainties (σ p,95 ) as a function of TS values.σ p,95 is defined by σ 2 p,95,stat + σ 2 p,sys , where σ p,sys is the systematic positional error which is 0.03 • for the KM2A component and 0.04 • for the WCDA component, respectively, considering the pointing error only (Detailed in Section 3.6).The dashed line is a trend for reference, considering a point source at a declination of 20 • with the index of 2.5 for WCDA component and of 3.5 for KM2A component.The marker size is scaled by the source radius.The marker color represents the photon spectral index.The source components with the significance of TS > 37 and the extension size of r 39 < 2 • are plotted.

Figure 3 .
Figure 3.The distribution of the extension sizes (r 39 ) of the source components detected by LHAASO.The size of the point-like source component is set to 0. Due to the limitation of statistical level and angular resolution, the measured minimum extension is approximately 0.17 • , leading to a gap of around 0.1 • .

Figure 4 .
Figure 4.The distribution of the SED parameter of the source components.Left: the distribution of the differential flux (E 2 0 N 0 ).The reference energy E 0 is 3 TeV or 50 TeV for the WCDA or KM2A component, respectively.Right: the distribution of the photon spectral index (Γ) for WCDA and KM2A components.

Figure 5 .
Figure5.1LHAASO source differential flux at 3 TeV (F 3 ) or 5 TeV (F 50 ) as a function of declination.The marker size represents the source size.The color indicates whether or not these sources are the new sources (seen in Section 5.2).Left: the integral sensitivity is shown at 3 TeV for three hypotheses: 1) point source with spectrum of E −2.0 ; 2) point source with spectrum of E −2.5 ; 3) 0.5 • gaussian source with spectrum of E −2.5 .The 3HWC sensitivity corresponds to the point-source search with the 3HWC data set(Albert et al. 2020).Right: the integral sensitivity is shown at 50 TeV for three hypotheses: 1) point source with spectrum of E −3.0 ; 2) point source with spectrum of E −3.5 ; 3) 0.5 • gaussian source with spectrum of E −3.5 .

Figure 6 .
Figure 6.Flux vs. source size and flux vs. photon spectral index.Top: 1LHAASO source differential flux (F 3 or F 50 ) as a function of source size.The error bar of r 39 represents the statistical uncertainty at a 95% confidence level.The dashed lines are integral sensitivity shown at 3 TeV or 50 TeV, assuming a point-like source.Bottom: 1LHAASO source photon spectral index as a function of the differential flux (F 3 or F 50 ).The definition of the new sources is seen in Section 5.2.
• are shown in Figures 11-12, in which most sources are concentrated in the inner Galactic plane.Eight individual sources are detected at Galactic latitude |b| > 12 • , as shown in

Figure 7 .
Figure7.The position and extension comparison between the WCDA and KM2A components.The six most significant sources are highlighted in red.Left: the position offsets (r d ) relative to the declinations.In the case where the declination of WCDA component is smaller than that of the KM2A component, we take the opposite value.The error bar is defined as σ 2 p,stat,95,W + σ 2 p,stat,95,K , where σ p,stat,95,W and σ p,stat,95,K are the statistical positional error of WCDA and KM2A components, respectively.The orange represents the systematic error just including the pointing error.Right: the extensions of the WCDA components relative to that of the KM2A components.The error bar represents the statistical uncertainty at a 95% confidence level.The orange band represents the systematic error just considering that of PSF.

Figure 13 .
Figure13.Among them, 1LHAASO J1653+3943 and 1LHAASO J1104+3810, which correspond to Mrk 421 and Mrk 501 respectively, have been significantly detected by the WCDA detector, however, no gamma-ray emission has been found at E>25 TeV by KM2A detector, due to significant absorption from interstellar radiation fields (ISRFs) and CMB.In addition to Mrk 421 and Mrk 501, there are another 34 sources detected by only one detector.This finding is reasonable considering the curved spectral shape, as illustrated in Figure8.The different shapes of the spectra can result in various relationships in the TS value between WCDA and KM2A detections.

Figure 8 .
Figure 8. TS value of KM2A component versus that of WCDA component.The reference dashed lines indicate the expected TS value for each detector.These values are calculated based on a point source which has a broken power-law spectral shape with a break at an energy of 25 TeV.

Figure 9 .
Figure 9. Distribution of the spectral index and the significance for sources at energy range from 1−25 TeV (left) and E > 25 TeV (right), respectively.The red points indicate the sources also with significance above 4σ at E >100 TeV.The extragalactic sources are excluded in these figures.The six most significant sources are labeled by their name.
The ρ( Ė) is counted using the pulsars nearby the candidate pulsar with Galactic latitude |b−b c | < 2.5 • and longitude |l − l c | < 10 • , where (b c , l c ) denotes the Galactic coordinates of the candidate pulsar.

Figure 10 .
Figure10.Pulsar spindown power Ė versus age for all pulsars of the ATNF catalog within the FOV of LHAASO and the 1LHAASO-associated pulsars with the chance probability P c less than 0.01.
, Department of Science and Technology of Sichuan Province, China No.2021YFSY0030, Project for Young Scientists in Basic Research of Chinese Academy of Sciences No.YSBR-061, and in Thailand by the NSRF via the Program Management Unit for Human Resources & Institutional Development, Research and Innovation (No. B37G660015).
. Xi and S.Z.Chen led the drafting of text and performed the data analysis of KM2A, S.C. Hu and M. Zha conducted the data analysis of WCDA.Y.Y.Guo and J.Y.He provided the cross-check.G.M. Xiang illustrated some figures in this article.Zhen Cao, the spokesperson of the LHAASO Collaboration, coordinated the specific working group for this paper involving all corresponding authors.All other authors participated in data analysis, including detector calibration, data processing, event reconstruction, data quality checks, and various simulations, and provided comments on the manuscript.

Figure 11 .
Figure 11.LHAASO significance map within the region 10 • ≤ l ≤ 115 • , |b| ≤ 12 • .Top: WCDA (1 TeV < E < 25 TeV) significance map.Middle: KM2A (E > 25 TeV) significance map.Bottom: KM2A (E > 100 TeV) significance map.In this figure and following, the LHAASO source are represented by gray crosses and white labels.The LHAASO sources for which the WCDA component has higher significance are plotted in the top panel.The LHAASO sources for which the KM2A component has higher significance are plotted in the middle panel.Meanwhile, UHE sources are shown again in the bottom panel.

Table 1 .
Property of each Analysis BinDetectorN hit /E rec (TeV) ϕ 68 ( • ) E M C • WCDA component and ∼ 0.85 • KM2A component.The nearest known TeV source is 2HWC J1852+013* (0.55 • away), which is suffers from GDE impact and is without an extended size measurement.We claim this 1LHAASO source as a new TeV source due to its large extension.The middle-aged pulsar, PSR J1853+0056 (d = 3.84 kpc, Ė = 4.03 × 10 34 erg s −1 , τ c = 204 kyr) is the nearest pulsar, 0.31 • from the position of the KM2A component.The 1LHAASO source is possibly a TeV halo associated with this middle-aged pulsar.On the other hand, an extended GeV source 4FGL J1852.4+0037e(namely Kes 79), modeled by 0.63 • disk shape and identified as the SNR or PWN type by the F ermi-LAT group, is in agreement with the position and size of the 1LHAASO source.
• WCDA component and r 39 ≈ 0.78 • KM2A component.It is resolved from the published LHAASO source LHAASO J1958+2845.At 0.36 • from the position of the WCDA component, a shell type SNR with radio size of 31 ′ × 25 ′ is found.

Table 2 .
1LHAASO source catalog Table 2 continued on next page Table 2 (continued)

Table 2
continued on next page

Table 2
continued on next page