Discovery of the Ultrahigh-energy Gamma-Ray Source LHAASO J2108+5157

Located at 4410 m above sea level in the Sichuan province of China, LHAASO is a new-generation complex extensive air shower (EAS) array (Cao 2010). It consists of three subarrays: the KM2A, the Water Cherenkov Detector Array (WCDA), and the Wide-Field Air Cherenkov Telescope Array (WFCTA) (He 2018). As a major part of LHAASO, KM2A is composed of 5195 electromagnetic particle detectors (EDs) and 1188 muon detectors (MDs), which are distributed in an area of 1.3 km². Each ED (Lv et al. 2018) consists of four plastic scintillation tiles covered by a 0.5 cm thick lead plate to convert the gamma rays to electron–positron pairs and improve the angular resolution of the array. The EDs detect the electromagnetic particles in the shower that are used to reconstruct information on the air shower, such as the primary direction, core location, and energy. Each MD includes a cylindrical water tank with a diameter of 6.8 m and a height of 1.2 m. The tank is buried under 2.5 m of soil to shield against the high-energy electrons/positrons and photons of the showers. The MDs are designed to detect the muon component of showers, which is used to discriminate between gamma ray– and hadron-induced showers.

Half of the KM2A array, equipped with 2365 EDs and 578 MDs, has been in service since 2019 December. For the half-array, a trigger is generated when 20 EDs are fired within a 400 ns window. This results in a 1 kHz event trigger rate. For each triggered event, the parameters of the air shower, like the direction, core, and gamma/hadron separation variables, are derived from the recorded hit time and charge. The core resolution and angular resolution (68% containment) of the half-array are energy- and zenith-dependent (Aharonian et al. 2020). The core resolution ranges from 4−9 m for events at 20 TeV to 2−4 m for events at 100 TeV, and the angular resolution ranges from 05−08 to 02−03. The energy resolution is about 24% at 20 TeV and 13% at 100 TeV for showers with a zenith angle less than 20°. Thanks to the good measurement of the muon component, KM2A has reached a very high rejection power of the hadron-induced showers. The rejection power is about 10³ at 25 TeV and greater than 4 × 10³ at energies above 100 TeV.

The same simulation data presented in Aharonian et al. (2020) are adopted here to simulate the detector response by reweighting the shower zenith angle distribution to trace the trajectory of LHAASO J2108+5157. For this sample, CORSIKA (version 7.6400; Heck et al. 1998) is used to simulate air showers, and a specific software, G4KM2A (Chen et al. 2019), is used to simulate the detector response. The energy of the gamma rays is sampled from 1 TeV to 10 PeV, and the zenith angle is sampled from 0° to 70°.

The pipeline of KM2A data analysis presented in Aharonian et al. (2020) is designed for surveying the whole sky in the range of decl. from −15° to 75° and the corresponding measurements for the source morphology and energy spectrum. The same pipeline is directly adopted in this analysis. The LHAASO-KM2A data were collected from 2019 December 27 to 2020 November 24. The final live time used for the analysis is 308.33 days, corresponding to a 94% duty cycle. The effective observation time on source LHAASO J2108+5157 is 2525.8 hr. After the pipeline cuts presented in Aharonian et al. (2020), the number of events used in this analysis is 1.2 × 10⁸.

The sky map in celestial coordinates (R.A. and decl.) is divided into a grid of 01 × 01 filled with the number of detected events according to their reconstructed arrival directions (event map). The "direct integration method" (Fleysher et al. 2004) is adopted to estimate the number of CR background events in each grid (background map). The background map is then subtracted from the event map to obtain the source map. The significance levels of sources are computed using a likelihood analysis given a specific source geometry. It takes into account a given source model, the data and background maps, and the expected detector response and calculates a binned Poisson log-likelihood value. A likelihood ratio test is built to compute the test statistic (TS):

$$\begin{eqnarray}&&{\rm{TS}}=2{ln}\displaystyle \frac{{L}_{s+b}}{{L}_{b}}.\end{eqnarray} \tag{ 1 }$$

Here L_b is the maximum likelihood of the only background hypothesis, and L_s+b is the maximum likelihood of the signal-plus-background hypothesis. Taking Wilks's theorem (Wilks 1938) into account, the TS is distributed as the χ² distribution with the number of degrees of freedom equal to the difference of the number of free parameters between the hypotheses. The usual source discovery test includes the null hypothesis as "there is no source" and the alternative hypothesis as "there is a source with a flux normalization X." There is only one degree of freedom. The significance simplifies as $\sqrt{{\rm{TS}}}$.

We estimate the spectral energy distribution (SED) of LHAASO J2108+5157 with the forward-folding method described in Aharonian et al. (2020). The SED of this source is assumed to follow a power-law spectrum dN/dE = ϕ₀(E/20TeV)^−α . The best-fit values of ϕ₀ and α are obtained by the least-squares fitting method.

The significance maps around LHAASO J2108+5157 in both energy ranges of 25−100 and >100 TeV are shown in Figure 1, smoothed using the point-spread function (PSF) of KM2A in the corresponding energy range. The source is detected with a statistical significance of 9.6σ and 8.5σ, respectively. Given the maximum trials of 3.24 × 10⁶ in the whole sky of −15° < decl. <75°, the posttrial significance at >100 TeV is 6.4σ. The position of LHAASO J2108+5157 is determined by fitting a 2D symmetrical Gaussian model, taking into account the KM2A PSF using events with E_rec > 25 TeV. The centroid of the Gaussian corresponding to the location of LHAASO J2108+5157 is found to be R.A. = 31722 ± 007_stat, decl. = 51.95° ± 005_stat (J2000), which is coincident with the location of events with E_rec > 100 TeV (Cao et al. 2021).

Zoom In Zoom Out Reset image size

To study the morphology of the source, a 2D symmetrical Gaussian template convolved with the KM2A PSF is used to fit the data with E_rec > 25 TeV using Equation (1). The source is found to be pointlike (TS = 118.79), but a slightly extended morphology (x = 121.48) cannot be ruled out due to the limited statistics and uncertainty on the KM2A PSF. An upper limit on the extension of the source is calculated to be 026 at a 95% confidence level (CL). Figure 2 shows the measured angular distribution based on KM2A data related to LHAASO J2108+5157 using events with E_rec > 25 TeV. The distribution is generally consistent with the PSF obtained using Monte Carlo (MC) simulations (χ²/ndf = 9.1/10).

Zoom In Zoom Out Reset image size

The data sets are divided into eight energy bins (Table 1 in Appendix A). The total number of photon-like events detected above 25 and 100 TeV is 140 and 18, respectively. There is one photon-like event with an energy of 434 ± 49 TeV in the last bin. Assuming a single power-law form, the differential energy spectrum of gamma-ray emission from LHAASO J2108+5157 is derived (Figure 3). It can be described by a single power law from 20 to 500 TeV as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{dN}}{{dE}}=(1.59\pm {0.35}_{\mathrm{stat}})\\ \times \ {10}^{-15}{\left(\displaystyle \frac{E}{20\mathrm{TeV}}\right)}^{-2.83\pm 0.18}({\mathrm{TeV}}^{-1}\ {\mathrm{cm}}^{-2}\ {{\rm{s}}}^{-1}).\end{array}\end{eqnarray} \tag{ 2 }$$

The χ² test for goodness of fit (χ²/ndf) equals 4.26/5, while it equals 3.8/4 in a power law with an exponential cutoff hypothesis. The improvement is not significant, with only about 0.7σ compared to fitting by a pure power law. The systematic uncertainty comes largely from the atmospheric model used in the simulation. According to the variation of the event rate during the operation, the overall systematic uncertainty of the flux is about 7%, and that of the spectral index is 0.02 (Aharonian et al. 2020).

In this analysis, LHAASO J2108+5157 is a pointlike source with a 95% position uncertainty of 014. We first searched for the gamma-ray counterpart within 014 of the center of LHAASO J2108+5157 in the VHE catalog TeVCat. ³³ No VHE counterpart of LHASSO J2108+5157 is found, even when enlarging the searching radius to 05. In the fourth Fermi-LAT Source Catalog (4FGL; Abdollahi et al. 2020), the high-energy point source 4FGL J2108.0+5155 is spatially coincident with LHAASO J2108+5157 at an angular distance of ∼013. The coincidence possibility of a high-energy point source being found in the region of LHAASO J2108+5157 is about 0.006. We dedicated an analysis to 4FGL J2108.0+5155 using ∼12.2 yr of Fermi-LAT data (see Appendix B for details). At a high significance level (7.8σ), 4FGL J2108.0+5155 is spatially extended (namely, 4FGL J2108.0+5155e), and the extension of a 2D Gaussian model is ∼048. The extrapolation of the spectrum of 4FGL J2108.0+5155e predicts a differential flux of 4.4 × 10⁻¹³ erg cm⁻² s⁻¹ at 10 TeV, a factor of 10 lower than that of LHAASO J2108+5157. Considering that the angular size of 4FGL J2108.0+5155e is about two times larger than the 95% upper limit extension (UL_ext,95% = 026) of LHAASO J2108+5157, the physical coincidence between the two sources cannot be clearly identified. Therefore, we derived the upper limit flux with the same spatial template as LHAASO J2108+5157 (σ = 026), centered at the position of LHAASO J2108+5157 above 10 GeV, which is used to limit the 10 GeV–1 TeV emission associated with LHAASO J2108+5157 (as shown in Figure 5).

Zoom In Zoom Out Reset image size

For the X-ray observation, Swift-XRT surveyed this region with an exposure of 4.7 ks (Stroh & Falcone 2013). However, no X-ray counterparts are found within 026 of the center of LHAASO J2108+5157. The closest X-ray source is the eclipsing binary RX J2107.3+5202, with a separation of 03. At radio wavelengths, Canadian Galactic Plane Survey carried out a high-resolution survey of the 408 and 1420 MHz continue emission covering the LHAASO J2108+5157 region. As shown in Figure 4, an extended radio source is within the region of UL_ext,95% of LHAASO J2108+5157, so it seems plausible that they are associated with the star-forming region nearby. The radio upper limits are derived from the 026 region centered at the position of LHAASO J2108+5157.

Zoom In Zoom Out Reset image size

We use the 2.6 mm CO line survey (Dame et al. 2001) to search for molecular cloud clumps in the direction of LHAASO J2108+5157. Two peaks at ∼−13 and ∼−2 km s^–1 are identified as the molecular clouds [MML2017]4607 and [MML2017]2870 (Miville-Deschênes et al. 2017), respectively (see Appendix C for details). As shown in Figure 4, LHAASO J2108+5157 lies near the center of the molecular cloud [MML2017]4607 (Miville-Deschênes et al. 2017), and this cloud is within the upper limit of the extension of LHAASO J2108+5157. The average angular radius of the cloud [MML2017]4607 is 0236, with a mass of 8469 M_⊙ at a distance of ∼3.28 kpc. The number density is estimated to be ∼30 cm⁻³. The coincidence possibility of a molecular cloud being found in the region of LHAASO J2108+5157 is about 0.04.

Most possible candidates to accelerate particles up to hundreds of TeV include SNRs, pulsar wind nebulae (PWNe), and young stellar clusters. Based on those catalogs collected by SIMBAD, ³⁴ we searched for the possible accelerators within 08 from the center of the LHAASO source. In spite of missing SNRs and PWNe, we found two young star clusters, Kronberger 82 (Kronberger et al. 2006) and Kronberger 80 (Kharchenko et al. 2016).

Because LHAASO J2108+5157 and the molecular cloud [MML2017]4607 have a spatial coincidence, there is a preference for hadronic origin. The UHE emission could be produced by protons accelerated up to PeV colliding with the ambient dense gas. We use the NAIMA package (Zabalza 2015) to estimate the parent particle spectrum to best reproduce the observed gamma-ray spectrum from the LHAASO J2108+5157 region. The observed gamma rays are attributed to the decay of π⁰ mesons produced in inelastic collisions between accelerated protons and target gas in [MML2017]4607. For the energy distribution of the parent particles, we assume an exponential cutoff power-law form. The index is fixed to 2, which is predicted in the standard diffusive shock acceleration. We obtain a cutoff energy of about 600 TeV. The total energy of the CR protons is $2\times {10}^{48}{\left(\tfrac{n}{30{\mathrm{cm}}^{-3}}\right)}^{-1}{\left(\tfrac{D}{3.28\mathrm{kpc}}\right)}^{2}$ erg, which is less than 10% of the typical kinetic energy of supernova explosions and consistent with the energy content of escaping CRs from SNRs (Cui et al. 2018), where n is the gas density and D is the distance to the source. Considering the soft GeV spectrum of 4FGL J2108.0+5155e in the lower energy band and the hard TeV spectrum of LHAASO J2108+5157, together with the different spatial extensions, we suggest that the much-extended gamma-ray emission from 4FGL J2108.0+5155e may be associated with an old SNR, like W28. And the pointlike gamma-ray emission from LHAASO J2108+5157 may be produced by the interaction between the escaping CRs from the SNR and the molecular cloud clumps, like the three TeV sources in the south of W28 detected by HESS (Aharonian et al. 2008). In addition, the young massive stars may also operate as proton PeVatrons with a dominant contribution to the Galactic CRs (Aharonian et al. 2019). Two young star clusters, Kronberger 80 and Kronberger 82, indicated in Figure 4, are located near LHAASO J2108+5157 with angular distances of 062 and 045, respectively. Kronberger 80, at a distance of 4.9 kpc (Kharchenko et al. 2016), is about 1.6 kpc away from the cloud [MML2017]4607, implying that Kronberger 80 is unlikely to interact with [MML2017]4607 to radiate UHE gamma-ray photons. Considering the absence of distance measurements for Kronberger 82, we cannot conclude whether or not Kronberger 82 is the source near [MML2017]4607 for explaining the UHE gamma-ray emission.

Nevertheless, UHE gamma rays can also be generated by high-energy electrons, as shown in Figure 5. According to the previous observations, a large fraction of the VHE sources are identified to be PWNe. Recent LHAASO observations also found that most of the sources detected with energy above 100 TeV are spatially coincident with identified pulsars (Cao et al. 2021). A likely scenario for LHAASO J2108+5157 would be that the UHE gamma-ray emission stems from a PWN powered by a yet-unknown pulsar. According to the spectrum of LHAASO J2108+5157, the gamma-ray luminosity at 20−500 TeV is estimated to be $1.4\times {10}^{32}{\left(\tfrac{D}{1\mathrm{kpc}}\right)}^{2}$ erg s⁻¹. It can be reasonably powered by a typical gamma-ray pulsar, which usually has a spin-down luminosity ranging from 10³⁵ to 10³⁹ erg s⁻¹ (Mattana et al. 2009). The PWN scenario for LHAASO J2108+5157 suggests that the UHE gamma-ray emission is produced by the inverse Compton scattering process of electrons. We assumed that the primary electron spectrum follows a power law with an exponential cutoff ${dN}/{dE}\propto {E}^{-{\rm{\Gamma }}}\exp (-E/{E}_{c})$. Considering the absence of multiwavelength data, the spectral index of the electrons is adopted to be 2.2, and the magnetic field strength is adopted to be 3 μG, both of which are reasonable for typical PWNe. The cutoff energy of electrons is set to be 200 TeV to explain the UHE gamma-ray spectrum of LHAASO J2108+5157. It corresponds to a synchrotron energy loss timescale of 10⁴ yr with a magnetic field strength of 3 μG, and the total energy of electrons with energy above 1 GeV is estimated to be $1\times {10}^{46}{\left(\tfrac{D}{1\mathrm{kpc}}\right)}^{2}$ erg. Because of the absence of a pulsar counterpart, the PWN scenario remains uncertain.

Zoom In Zoom Out Reset image size