Measurement of the Gamma-Ray Energy Spectrum beyond 100 TeV from the HESS J1843$-$033 Region

HESS J1843$-$033 is a very-high-energy gamma-ray source whose origin remains unidentified. This work presents, for the first time, the energy spectrum of gamma rays beyond $100\, {\rm TeV}$ from the HESS J1843$-$033 region using the data recorded by the Tibet air shower array and its underground muon detector array. A gamma-ray source with an extension of $0.34^{\circ} \pm 0.12^{\circ}$ is successfully detected above $25\, {\rm TeV}$ at $(\alpha,\, \delta) = (281.09^{\circ}\pm 0.10^{\circ},\, -3.76^{\circ}\pm 0.09^{\circ})$ near HESS J1843$-$033 with a statistical significance of $6.2\, \sigma$, and the source is named TASG J1844$-$038. The position of TASG J1844$-$038 is consistent with those of HESS J1843$-$033, eHWC J1842$-$035, and LHAASO J1843$-$0338. The measured gamma-ray energy spectrum in $25\, {\rm TeV}<E<130\, {\rm TeV}$ is described with ${\rm d}N/{\rm d}E = (9.70\pm 1.89)\times 10^{-16} (E/40\, {\rm TeV})^{-3.26\pm 0.30}\, {\rm TeV}^{-1} {\rm cm}^{-2} {\rm s}^{-1}$, and the spectral fit to the combined spectra of HESS J1843$-$033, LHAASO J1843$-$0338, and TASG J1844$-$038 implies the existence of a cutoff at $49.5\pm 9.0\, {\rm TeV}$. Associations of TASG J1844-038 with SNR G28.6$-$0.1 and PSR J1844-0346 are also discussed in detail for the first time.


INTRODUCTION
The cosmic-ray energy spectrum has the so-called knee at about 4 PeV (Kulikov & Kristiansen 1958;Hörandel 2004;Amenomori et al. 2008), and its origin is still unknown. Searches for a PeVatron, an accelerator of such PeV cosmic rays in the Galaxy, are performed through observations of sub-PeV gamma rays (E > 100 TeV) stemming from decays of neutral pions generated by hadronic interactions of cosmic rays with the interstellar medium. So far about a dozen of sources have been discovered in the sub-PeV energy range (for example, see Amenomori et al. (2019); Abeysekara et al. (2020); Cao et al. (2021)), and one of these sources, a supernova remnant (SNR) G106.3+2.7, is a candidate for a PeVatron because of good spatial correlation between the distributions of gamma rays and molecular clouds (Albert et al. 2020;Amenomori et al. 2021a). However, the maximum energy of cosmic-ray protons accelerated by SNR G106.3+2.7 is estimated at ≃ 0.5 PeV lower than the knee energy range, and the robust detection of a PeVatron remains inconclusive. In many cases, gamma rays of leptonic origin also significantly contribute to the total gammaray flux of sources, making the situation confusing. On the other hand, the observation of sub-PeV Galactic diffuse gamma rays provides evidence that PeVatrons exist (and/or existed) in the Galaxy (Amenomori et al. 2021b). The result is consistent with the theory presented by Lipari & Vernetto (2018) arguing that gamma rays of hadronic origin dominate the total Galactic diffuse gamma-ray flux in the sub-PeV range. These facts motivate us to directly specify PeVatrons in the Galaxy.
HESS J1843-033 is one of the unidentified gamma-ray sources discovered by the H.E.S.S. Galactic plane survey (Hoppe 2008; H.E.S.S. collaboration 2018a). The energy spectrum was measured up to 30 TeV and described with a power-law function with an index of ≃ −2.2. HAWC also reported the discovery of 2HWC J1844-032 and eHWC J1842-035 near HESS J1843-033 (Abeysekara et al. 2017b(Abeysekara et al. , 2020, and for the latter, they estimated the integral gamma-ray flux above 56 TeV. Moreover, LHAASO discovered LHAASO J1843-0338 and provided a flux point at 100 TeV (Cao et al. 2021). Theoretically, the nature of the gamma-ray emission has been discussed by Sudoh et al. (2021) and Huang & Li (2022), but their claims are still inconclusive because the gamma-ray energy spectrum is not systematically measured above the several tens of TeV range and its characteristics are thus not revealed yet. This paper presents the first observational result of the energy spectrum of gamma rays beyond 100 TeV from the HESS J1843-033 region and discusses the origin of the gamma-ray emission considering its possible associations with nearby celestial objects.

EXPERIMENT AND DATA ANALYSIS
The Tibet air shower array (AS array) has been operating since 1990 in Yangbajing (90. • 522 E, 30. • 102N, 4, 300 m a.s.l.) in Tibet, China (Amenomori et al. 1992(Amenomori et al. , 1999(Amenomori et al. , 2009 and currently consists of a surface AS array and an under- Table 1. Positions and extensions of TASG J1844-038 and nearby gamma-ray sources. α and δ are the right ascension and declination, respectively, in the J2000 equatorial coordinates. R0.68 denotes the error radius of a source position at the 68% confidence level (see Appendix B). For the source extension of TASG J1844-038, see the text. Numbers enclosed in the parentheses in the Angular distance to TASG J1844-038 column denote the significance of a positional deviation between TASG J1844-038 and a source evaluated with their R0.68's. Data of the nearby sources are taken from H.E.S.S. collaboration (2018a), Abeysekara et al. (2020), andCao et al. (2021). ground muon detector array (MD array, Amenomori et al. (2019)). The AS array comprises 597 plastic scintillation detectors, each with a detection area of 0.5 m 2 , and covers an area of 65, 700 m 2 . The MD array with a total area of 3, 400 m 2 consists of water-Cherenkov type detectors to collect Cherenkov light emitted by penetrating shower muons in the water layer. Each detector of the MD array is made of concrete and is located 2.4 m underground. The soil overburden and the concrete ceiling of the MD detectors correspond to the total thickness of ≃ 550 g cm −2 , allowing only shower muons with E ≳ 1 GeV to reach the water layer of the array for vertical incidence. Using the MD array thus enables us to improve the experiment's sensitivity to sub-PeV gamma rays by more than one order of magnitude by rejecting more than 99.9% of hadronic cosmic rays while keeping 90% of gamma rays in the sub-PeV energy range (Sako et al. 2009;Amenomori et al. 2019). A trigger is issued when any four plastic scintillation detectors detect more than 0.6 particles within a time window of 600 ns. In this work, data taken from 2014 February to 2017 May (719 live days) are analyzed. Reconstruction methods for recorded AS events and selection criteria imposed on the events are the same as Amenomori et al. (2019) except for the following two points. First, the maximum zenith angle of the incoming direction of analyzed events is changed from 40 • to 50 • because the meridian zenith angle of HESS J1843-033 is ≃ 33 • at the site of the experiment. This extension helps us increase the exposure to the HESS J1843-033 region by ≃ 70% and improve statistics of gamma rays from that region. Second, the muon-cut condition optimized for this work is used (see Appendix A). After the selection, events are analyzed, binned in energy with five logarithmically equal bins per decade. For the Monte Carlo simulation performed in this work, see Appendix A.

RESULTS
A gamma-ray source is successfully detected above 25 TeV near HESS J1843-033 with a statistical significance of 6.2 σ calculated from Equation (17) of Li & Ma (1983). Figure 1 shows the significance map of the observed gamma rays. A two-dimensional maximum likelihood analysis is performed to estimate the source center assuming an axisymmetric Gaussian distribution for the distribution of the gamma rays. The resultant center is (α, δ) = (281. • 09 ± 0. • 10, −3. • 76 ± 0. • 09) in the J2000 equatorial coordinates. The source is named TASG J1844-038, and hereafter this name is used accordingly. Table 1 presents the positions of TASG J1844-038 and nearby gamma-ray sources and shows that the position of TASG J1844-038 is statistically consistent with those of HESS J1843-033, eHWC J1842-035, and LHAASO J1843-0338. On the other hand, TASG J1844-038 deviates from HESS J1844-030 and HESS J1846-029 at the 3.2 σ and 4.5 σ levels, respectively, making its associations with these sources unlikely. Figure 2 presents the distribution of events above 25 TeV around the center of TASG J1844-038. The horizontal axis ϕ 2 denotes the square of the angle between the center of TASG J1844-038 and the incoming direction of events. The blue histograms are constructed from background events in OFF regions plus Monte Carlo gamma-ray events from a point source normalized to the number of excess counts in the ON-source region. The source extension is estimated by fitting the following Gaussian function to the observed number of events: where A is a normalization constant, σ ext the extension of TASG J1844-038, σ psf = 0. • 28 the radius of the pointspread function for gamma rays above 25 TeV that follow a power-law energy spectrum with an index of 3.0, and N bg = 29.4 the number of background events. The best-fit result is shown by the black curve in Figure 2, leading to σ ext = 0. • 34 ± 0. • 12 with a reduced chi-square of χ 2 /d.o.f. = 39.5 / 38. Figure 3 shows the gamma-ray energy spectrum measured in this work. A gamma-ray flux is calculated only if the statistical significance of gamma-ray detection in each energy bin exceeds 2 σ, and otherwise, a 95% upper limit is calculated. Our results in 25 TeV < E < 130 TeV can be fitted with a power-law function of dN/dE = (9.70±1.89)×10 −16 (E/40 TeV) −3.26±0.30 TeV −1 cm −2 s −1 with χ 2 /d.o.f. = 2.1/2. As a possible systematic uncertainty, the uncertainty in the source extension σ ext affects the flux normalization by 19% and the spectral index by 0.02. The absolute energy-scale uncertainty of 12% (Amenomori et al. 2009) also affects the flux normalization by 39%. Summing up these two in quadrature, the total systematic uncertainty in the flux normalization is calculated as 43%. The spectrum of TASG J1844-038 smoothly connects to that of HESS J1843-033 and is consistent with the fluxes of eHWC J1842-035 and LHAASO J1843-0338, supporting the common origin for gamma rays from these sources. A fit of a simple power-law function to the combined spectra of HESS J1843-033, LHAASO J1843-0338, and TASG J1844-038 is disfavored at the 5.0 σ level, and the spectrum is better described with the following power-law function with an exponential cutoff: where the best-fit parameters are N 0 = (3.57±0.26)×10 −12 TeV −1 cm −2 s −1 , Γ = 2.02±0.06 and E cut = 49.5±9.0 TeV with χ 2 /d.o.f. = 10.4/8. The flux of eHWC J1842-035 is not included in the combined fit because only the integral flux above 56 TeV is presented by Abeysekara et al. (2020) (the differential flux of eHWC J1842-035 shown in Figure  3 is calculated from the integral flux with a spectral index of −2.7 assumed in their paper). The time variation of the gamma-ray flux is also searched by equally dividing the data into four independent periods and calculating the integral gamma-ray flux in each period. A constant function is fitted to the gamma-ray flux points over the four periods, and the best-fit value is (4.65 ± 0.89) × 10 −14 cm −2 s −1 above 25 TeV with χ 2 /d.o.f. = 1.9/3, consistent with the null hypothesis that no time variation of the gamma-ray flux takes place during our observation.  2021)). The flux of eHWC J1842-035 is calculated from the integral flux above 56 TeV assuming a differential spectral index of −2.7. The magenta dotted curve shows the best-fit power-law function with an exponential cutoff fitted to the combined spectra of HESS J1843-033, LHAASO J1843-0338, and TASG J1844-038.

DISCUSSION
Among the sources near TASG J1844-038, SNR G28.6-0.1, whose center is within the positional error circle of TASG J1844-038 as shown in Figure 4, is one of the prominent candidates for the counterpart. A radio complex G28.60-0.13 was first discovered by Helfand et al. (1989), some regions of which (regions C and F) were found to be emitting non-thermal radio emissions. Subsequently, an X-ray counterpart AX J1843.8-0352 was detected by ASCA and Chandra (Bamba et al. 2001;Ueno et al. 2003). The authors claimed that the observed non-thermal radio and X rays are synchrotron radiation of high-energy electrons accelerated by a shell-type SNR. The extension of AX J1843.8-0352 (mean diameter of 9 ′ (Ueno et al. 2003)) is smaller than that of TASG J1844-038 at the 2.3 σ level. The distribution of gamma rays of leptonic origin is expected to have approximately the same extension as that of the observed non-thermal X rays. The possible discrepancy thus suggests that, if SNR G28.6-0.1 is the counterpart of TASG J1844-038, radiation from hadronic cosmic rays accelerated by the SNR contributes to the observed gamma-ray emission. The maximum energy of accelerated cosmic-ray protons could be ≃ 500 TeV, roughly an order of magnitude larger than the cutoff energy of ≃ 50 TeV in the gamma-ray energy spectrum. In addition, the SNR is found to be interacting with molecular clouds in the velocity channel of around 86 km s −1 (Ranasinghe & Leahy 2018). Considering the peculiar motion of the Galactic objects (Reid et al. 2014), the distribution of molecular clouds in the integrated velocity channel from 75 km s −1 to 95 km s −1 is investigated with 12 CO (J = 1 − 0) line emission using the public FUGIN data (Umemoto et al. 2017). The result is shown in Figure 4 and some clouds are found to overlap with the TASG J1844-038 region. From the aforementioned reasons and the estimated age of SNR G28.6-0.1 (2.7 kyr by Bamba et al. (2001) or 19 kyr by Ranasinghe & Leahy (2018)), this SNR could be a similar system to a PeVatron candidate SNR G106.3+2.7 (Amenomori et al. 2021a) and thus could have been accelerating protons up to the PeV energy range in the past. Here the diffusion time scale of cosmic-ray protons in molecular clouds is estimated following Gabici et al. (2007): where R sys is the radius of the system, B the magnetic field strength, and χ the factor indicating the suppression of diffusion from the Galactic average. The extension of TASG J1844-038 corresponds to ≃ 60 pc at the distance of 9.6 kpc which is an estimated distance to SNR G28.6-0.1 (Ranasinghe & Leahy 2018), and the magnetic field can be assumed as 10 µG (Crutcher 1999). Cosmic-ray protons with E ≳ 250 TeV, which will generate π 0 -decay gamma rays with E ≳ 25 TeV in interactions with the clouds, can diffusively propagate up to the size of TASG J1844-038 within τ diff ∼ 2 kyr even if the case of χ ∼ 0.1 is considered. Under the same condition, the extension of HESS J1843-033 can be similarly explained by the gamma-ray emission from the accelerated protons with E ≃ 10 TeV diffusing for τ diff ∼ 5 kyr. These diffusion time scales do not conflict with the estimated age of SNR G28.6-0.1.
Another prominent candidate is the pulsar PSR J1844-0346 nearly centered at TASG J1844-038. This pulsar, also named 4FGL J1844.4-0345 (Abdollahi et al. 2020), was discovered by Fermi-LAT (Clerk et al. 2017), and its spindown luminosity, characteristic age, and spin period are 4.2 × 10 36 erg s −1 , 12 kyr, and 113 ms, respectively. Assuming the distance to the pulsar as 4.3 kpc (Devin et al. 2021), the gamma-ray energy flux of 1.1 × 10 −11 erg cm −2 s −1 in 1 TeV < E < 10 TeV (H.E.S.S. collaboration 2018a) is translated into the luminosity of 2.4 × 10 34 erg s −1 . The aforementioned properties of PSR J1844-0346 and the surrounding gamma-ray emission including its flat spectral index of 2.02 ± 0.06 in the TeV range are similar to those of "TeV pulsar wind nebulae (PWNe)" (H.E.S.S. collaboration 2018b) as proposed by Sudoh et al. (2021). The energies of electrons that generate gamma rays with E ≃ 50 TeV by Inverse Compton Scattering (ICS) off the cosmic microwave background (CMB) photons are ≃ 90 TeV (see the equation (10) of Hinton & Hofmann (2009)). Assuming the magnetic field as 3 µG and the diffusion coefficient of the 90 TeV electrons as 4.4 × 10 27 cm 2 s −1 following the observation of Geminga (Abeysekara et al. 2017a), the time scale of the electrons diffusing up to the size of TASG J1844-038 (≃ 26 pc at the assumed distance to PSR J1844-0346 of 4.3 kpc) is estimated at ∼ 8 kyr, which is within the characteristic age of PSR J1844-0346 and the cooling time of ∼ 11 kyr due to the synchrotron emission and ICS estimated from the equations (5) and (7) of Hinton & Hofmann (2009). Similarly, the extension of HESS J1843-033 in the TeV range can be explained by ICS off the CMB photons by ≃ 10 TeV electrons diffusing for ∼ 8 kyr. Devin et al. (2021) found no radio nor X-ray emission that indicates the existence of a PWN. Given the characteristic age of PSR J1844-0346, synchrotron emission from the PWN would not be bright enough to be observed due to the decay of the magnetic field (Tanaka & Takahara 2010). Future studies with wide field-of-view and high-sensitivity observations of X rays will be a key to constraining the TeV PWN scenario.
It should be noted that there are additional SNR candidates near TASG J1844-038 discovered by THOR (Anderson et al. 2017). Out of these candidates, G28.56+0.00, G28.64+0.20, and G28.78-0.44 overlap with TASG J1844-038, and future research on these SNRs is expected. Moreover, as pointed out by Devin et al. (2021), the star-forming region N49 (Dirienzo et al. 2012) is also located within TASG J1844-038 (see Figure 4). Several observations support the acceleration of cosmic rays in star-forming regions (for example, see Ackermann et al. (2011);Aharonian et al. (2019)), and detailed morphological studies of the gamma-ray emission of TASG J1844-038 by imaging atmospheric Cherenkov telescopes and its comparison with the stellar and gas distributions observed in other wavelength ranges will be needed to discuss a possible association with N49.
The origin of TASG J1844-038 is also discussed in detail for the first time assuming its associations with SNR G28.6-0.1 and PSR J1844-0346. If SNR G28.6-0.1 is assumed to be the counterpart, the nature of TASG J1844-038 can be explained by π 0 -decay gamma rays generated in hadronic interactions between adjacent molecular clouds and cosmic-ray protons with E ≲ 500 TeV that are accelerated by the SNR and diffusely propagate through the clouds. Given the similarities with SNR G106.3+2.7 in terms of the maximum energy of accelerated protons, the partial overlap of the gamma-ray distribution with molecular clouds, and the SNR's age, SNR G28.6-0.1 could have been a PeVatron and accelerating cosmic-ray protons up to the PeV energy range in the past. On the other hand, if associated with PSR J1844-0346, TASG J1844-038 can be explained by the gamma-ray emission from a TeV PWN generated by ICS off the CMB photons by high-energy electrons. Other nearby sources including SNR candidates and the star-forming region N49 should also be studied to investigate their contribution to the observed gamma-ray emission.  (Umemoto et al. 2017). The color scale shows the main beam brightness temperature. The solid and dashed white circles denote the extension of TASG J1844-038 with a radius of σext and the positional uncertainty at the 68% confidence level with a radius of R0.68, respectively (see Table 1). Positions and extensions of nearby sources are shown in the same way as in Figure 1 except that the orange hexagon denotes the star-forming region N49 (Dirienzo et al. 2012).
The collaborative experiment of the Tibet Air Shower Arrays has been conducted under the auspices of the Ministry of Science and Technology of China and the Ministry of Foreign Affairs of Japan. This work was supported in part by a Grant-in-Aid for Scientific Research on Priority Areas from the Ministry of Education, Culture, Sports, Science and Technology, and by Grants-in-Aid for Science Research from the Japan Society for the Promotion of Science in Japan. This work is supported by the National Key R&D Program of China (No. 2016YFE0125500), the Grants from the National Natural Science Foundation of China (No. 11533007, No. 11673041, No. 11873065, No. 11773019, No. 11773014, No. 11633007, No. 11803011, and No. 11851305), and the Key Laboratory of Particle Astrophysics, Institute of High Energy Physics, CAS. This work is also supported by the joint research program of the Institute for Cosmic Ray Research (ICRR), the University of Tokyo. S. Kato is supported by JST SPRING, Grant Number JPMJSP2108. This publication makes use of data from FUGIN, FOREST Unbiased Galactic plane Imaging survey with the Nobeyama 45-m telescope, a legacy project in the Nobeyama 45-m radio telescope. A total number of 1.1 × 10 8 gamma rays are generated in 300 GeV < E < 10 PeV using Corsika v7.4000 (Heck et al. 1998) following a simple power-law spectrum with an index of −2.0 along the path in the sky of a hypothetical source with δ = 0 • . These gamma rays are injected into the atmosphere, and the development of the extensive air showers is simulated until they reach the site of the experiment. Shower cores are randomly thrown over the circular region with a radius of 300 m centered at the Tibet AS array. The detector responses to the showers are simulated with Geant4 v10.0 (Agostinelli et al. 2003). The detector simulation includes the calculation of energy loss of shower particles in the plastic scintillation detectors, the soil, and the concrete layers of the MD array, and their Cherenkov-light emission in the water layer of the MD array. The output is fed into our analysis pipeline in the same way as experimental data. In the Monte Carlo data analysis, the gamma-ray events are appropriately weighted accounting for the spectral index of the source. The resultant angular and energy resolutions for 100 TeV gamma rays from HESS J1843-033 are estimated at ≃ 0. • 28 (50% containment) and ≃ 30%, respectively. Compared with the previous work (Amenomori et al. 2019), these resolutions worsen by ≃ 40% due to the large meridian zenith angle of HESS J1843-033 (≃ 33 • ) and the extension of the analyzed zenith-angle range from 40 • to 50 • (see Section 2).
The muon cut is optimized for this work using the Monte Carlo gamma-ray data and experimental background data. Following the equi-zenith angle method (Amenomori et al. 2003), the numbers of background events are averaged over 100 OFF regions opened around the ON-source region which contains the HESS J1843-033 region. Both the ON-source and OFF regions are circular regions with a radius of 0. • 7. The number of gamma-ray events is normalized to that of excess counts in the ON-source region. The gamma-ray and background events are then binned in Σρ which denotes the total number density of shower particles recorded with the plastic scintillation detectors of the AS array, and the optimum cut on the number of muons (ΣN µ ) is determined so that the figure of merit S/ √ S + B (S is the number of gamma-ray events and B that of background events) is maximized in each bin. The resultant muon cut requires events to satisfy ΣN µ < 1.8 × 10 −3 (Σρ/m −2 ) 1.1 or ΣN µ < 0.4. After applying the muon cut to events, the survival ratio of gamma-ray events and the rejection power of background events are estimated at ≃ 80% and ≳ 99.9%, respectively, at the gamma-ray equivalent energy of 100 TeV. Due to the tight muon cut to improve the figure of merit for the low gamma-ray flux, the gamma-ray survival ratio worsens by ≃ 10% from the previous work (Amenomori et al. 2019).