Revision of the GeV γ-Ray Emission in the Region of HESS J1813-178 with Fermi-LAT

HESS J1813-178 is one of the brightest and most compact TeV γ-ray sources, and whether its γ-ray emission is associated with supernova remnant (SNR), pulsar wind nebula (PWN), or young stellar cluster (YSC) is still under debate. By analyzing the GeV γ-ray data in the field of HESS J1813-178 using 14 yr of PASS 8 data recorded by the Fermi Large Area Telescope (Fermi-LAT), we report the discovery of three different sources with different spectra in this region. The hard source with a PL spectral index of 2.11 ± 0.08 has a small size extension, which is spatially and spectrally coincident with the TeV γ-ray emission from HESS J1813-178. CO observations display the dense molecular clouds surrounding HESS J1813-178 in the velocity range of 45–60 km s−1. The possible origins of the γ-ray emission from HESS J1813-178 are discussed, including SNR G12.820.02, the PWN driven by the energetic X-ray pulsar PSR J1813-1749 and YSC Cl 1813-178. However, none of them can be ruled out clearly. Note that the maximum energy of protons in the hadronic model should exceed a few hundred TeV, which makes HESS J1813-178 a promising PeVatron. A detailed LHAASO data analysis about the morphology and spectrum would be helpful to investigate the origin of the γ-ray emission in this region and test its PeVatron nature.


INTRODUCTION
Very high energy (VHE; >100 GeV) surveys of our Galaxy with the ground based γ-ray detectors, like H.E.S.S., MAGIC, HAWC and LHAASO, have uncovered a population of γ-ray sources in the TeV regime and revealed the very energetic cosmic accelerators in our Galaxy.More than two hundred of VHE γ-ray sources have been detected(H. E. S. S. Collaboration et al. 2018a;Albert et al. 2020;Cao et al. 2023), most of which have been identified to be PWNe, SNRs, binaries, TeV halos and so on.However, a large fraction of the sources still remain without firm association, and the multi-wavelength observations of these sources are crucial for revealing their nature and elucidate the origin of cosmic rays (CRs) around the knee energy range.
As one of the brightest and most compact objects detected by the H.E.S.S. Galactic Plane Survey (HGPS) observations, HESS J1813-178 was detected to be nearly point-like with an extension of σ = 2. ′ 2 (Aharonian et al. 2005(Aharonian et al. , 2006)).The TeV γ-ray emission from HESS J1813-178 shows a power-law spectrum with a rather hard photon index of 2.09±0.08,and the following MAGIC observations give the similar spectrum with an index of 2.1±0.2 in the energy range of 0.4-10 TeV (Albert et al. 2006).In H. E. S. S. Collaboration et al. (2018a), the extension of HESS J1813-178 was updated to be σ = 3. ′ 0 ± 0 ′ .24,and the spectrum above 1 TeV follows a power law with an exponential cutoff of ∼ 7 TeV.
In the First LHAASO Catalog of Gamma-Ray Sources, both the Water Cherenkov Detector Array (WCDA) and Kilometer Squared Array (KM2A) of LHAASO detected a γ-ray source in the region of HESS J1813-178, 1LHAASO J1814-1719u, in the energy range of 1-25 TeV and above 25 TeV, respectively (Cao et al. 2023).The spatial template of 1LHAASO J1814-1719u for WCDA is described by a two dimensional (2D) Gaussian with σ = 0.71 • ± 0.07 • , and its γ-ray spectrum in 1-25 TeV is modeled by a power law with an index of 2.83 ± 0.06.KM2A gives an upper limit of the extension of 1LHAASO J1814-1719u with σ = 0.27 • , together with the power law spectral index of 3.49 ± 0.31 above 25 TeV (Cao et al. 2023).HESS J1813-178 is positionally coincident with a shell-type SNR G12.82-0.02,which lies in the vicinity of a bright star forming region (SFR) W33 (Brogan et al. 2005).And non-thermal X-ray emission was also detected in this region (Brogan et al. 2005;Landi et al. 2006).Deep X-ray observations from XMM-Newton and Chandra revealed a PWN embedded near the center of SNR G12.82-0.02,which makes G12.82-0.02 to be a composite system (Funk et al. 2007;Helfand et al. 2007;Gotthelf & Halpern 2009).The PWN component is powered by an energetic pulsar, PSR J1813-1749, and its spin-down luminosity and characteristic age are Ė = 5.6 × 10 37 erg s −1 and τ c = 5600 yr, respectively (Halpern et al. 2012;Ho et al. 2020).A young stellar cluster at a kinematic distance of 4.8 kpc, Cl 1813-178, was discovered in the TeV γ-ray emission region of HESS J1813-178 (Messineo et al. 2008(Messineo et al. , 2011)).PSR J1813-1749/SNR G12.82-0.02 was first suggested to be associated with Cl 1813-178 at the same distance (Messineo et al. 2011).While the newest radio observation from Green Bank Telescope shows a large dispersion measure for PSR J1813-1749, which indicates a far distance of either 6.2 or 12 kpc according to the different models of the electron density distribution (Camilo et al. 2021).
The origin of the γ-ray emission from HESS J1813-178 is still not identified.The scenarios of PWN, SNR and stellar cluster are discussed in previous works.Using the NANTEN 12 CO(J=1-0) observations, Funk et al. (2007) detected a giant molecular cloud (MC) in the south of HESS J1813-178, which is suggest to be associated with SFR W33.With the X-ray and molecular observations, the scenarios in which the γ-rays of HESS J1813-178 are emitted from the shell of the SNR and one in which they are emitted from the central PWN are modeled.None of them can be clearly ruled out.Fang & Zhang (2010) also predicted that the TeV γ-ray emission is mainly originated from the PWN associated with PSR J1813-1749, although the contribution from the SNR shell could be enhanced with a denser medium.With the Fermi-LAT data, Araya (2018) found extended GeV γ-ray emission around HESS J1813-178 with a radius of ∼0.6 • .Such size is much larger than that of the TeV γ-ray emission of HESS J1813-178.The global GeV spectrum in the energy range of 0.5-500 GeV was fitted by a power-law with an index of 2.14±0.04,which is not consistent with the inverse Compton scattering (ICS) emission characteristic from leptons in a PWN.Araya (2018) argued that the extended GeV emission may be related to the SFRs around HESS J1813-178, like W33.
Considering the uncertain origin and the complexity of HESS J1813-178, the detailed analysis with more GeV observational data and the multi-wavelength observations would be helpful to investigate the origin of the γ-ray emission in this region, we re-analysed the Fermi-LAT PASS 8 data around HESS J1813-178 in this work.The data analysis procedure and results, including the spatial and spectral analysis, are shown in Section 2. In Section 3, the observation of the molecular clouds around HESS J1813-178 is introduced.We discuss the origin of the γ-ray emission from HESS J1813-178 based on the multi-wavelength observations in Section 4. The conclusion of this work is presented in Section 5.

Data Reduction
In the following analysis, the latest Pass 8 data with "Source" event class are selected from 2008 August 4 (MET 239557418) to 2022 August 4 (MET 681264005).The energies of events are cut between 500 MeV and 1 TeV to avoid a too large point spread function (PSF) in the lower energy band.Meanwhile, the maximum zenith angle is adopted to be 90 • to minimize the contamination from the Earth Limb.The analysis is performed in a 14 • × 14 • square region of interest (ROI) with the standard LAT analysis software ScienceTools.The binned likelihood analysis method with gtlike is adopted, together with the instrumental response function of "P8R3 SOURCE V3".The Galactic and isotropic diffuse backgrounds are modeled by gll iem v07.fits and iso P8R3 SOURCE V3 v1.txt, respectively.All sources listed in the incremental version of the fourth Fermi-LAT source catalog (4FGL-DR4; Abdollahi et al. 2022;Ballet et al. 2023) within a radius of 20 • from the ROI center, together with the two diffuse backgrounds, are included in the model, which is generated with the user-contributed script make4FGLxml.py.

Spatial Analysis
In 4FGL-DR4 catalog, there is an extended GeV γ-ray source (4FGL J1813.1-1737e)described by an uniform disk with a radius of 0.6 • , which is regarded as the potential GeV counterpart of HESS J1813-178 (Araya 2018;Abdollahi et al. 2022).However, the size of 4FGL J1813.1-1737e is much larger than that of TeV γ-ray emission from HESS J1813-178 (σ = 3. ′ 0; H. E. S. S. Collaboration et al. 2018a).To investigate the GeV γ-ray emission around HESS J1813-178, we first refit the centroid and size of the uniform disk for 4FGL J1813.1-1737e with Fermipy, a PYTHON package that automates analyses with ScienceTools (Wood et al. 2017).The best-fit centroid position and 68% containment radius (R 68 ) are fitted to be R.A. = 273 • .419± 0 • .019,decl.= −17 • .625 ± 0 • .018,and 0 • .475+0 • .015−0 • .014, respectively.With the refit spatial model (Model 1), we get its GeV spectrum by dividing the data into 12 equal logarithmic energy bins from 500 MeV to 1 TeV.The likelihood fitting analysis is done for each energy bin and the 95% upper limits are calculated for the energy bins with Test Statistic (TS; Lande et al. 2012) values smaller than 5.0.The resulting spectral energy distribution (SED) of 4FGL J1813.1-1737e is shown in Figure 1, and an obvious spectral upturn is shown below an energy of about 10 GeV compared with the results from Araya (2018).
To test whether the upturn spectrum is intrinsic or is due to two overlapping components, we do the same likelihood fitting using the events with energies of 500 MeV -10 GeV (low energy) and 10 GeV -1 TeV (high energy), respectively.In the low energy range, we only select the "PSF3" events with the best quality of the reconstructed direction to reduce the uncertainties caused by large PSF.In the high energy range, all events are collected in order to have sufficient counts for the detailed analysis.For each analysis, we first create a TS map with all sources (except 4FGL J1813.1-1737e)included in the best-fit model (hence subtracted from the map), which are shown in Figure 2. The TS maps in the low and high energy bands show a clear difference of the centroids and sizes of the γ-ray emission.We thus expect two different γ-ray components in this region, which are labeled as SrcA for the low energy band and SrcB for the high energy band.To determine the morphology of the γ-ray emission from SrcA/SrcB, we used Fermipy to test the size of the uniform disk with the events of 500 MeV -10 GeV ("PSF3") and 10 GeV -1 TeV ("Source"), respectively.The spectrum of SrcA or SrcB are first assumed to be a power-law (PL; dN/dE ∝ E −α ) model.In the low energy band, the centroid position and R 68 of SrcA are fitted to be R.A. = 273 −0 • .017, respectively.The spatial extension of SrcB is more compact than that of SrcA, and the GeV γ-ray emission of SrcB concentrates in the TeV γ-ray emission region of HESS J1813-178 (Aharonian et al. 2006), while the fitted γ-ray spectrum of SrcB in the high energy band is hard with an index of 1.97 ± 0.17.In addition, the TS map in the high energy range shows an additional point source in the north region of SrcB, which is not included in 4FGL-DR4 catalog, but is in 4FGL-DR2 catalog (4FGL J1814.1-1710;Ballet et al. 2020).The point source (hereafter labeled as SrcC) is suggested to be associated with SNR G13.5+0.2 in 4FGL-DR2 catalog.Then we use the updated model including the two extended spatial templates of SrcA and SrcB (Model 2) and the additional SrcC (Model 3) to fit the data in the energy range of 500 MeV -1 TeV.Meanwhile, the models for 4FGL J1813.1-1737e and SrcA/B with two dimensional Gaussian template (Model 4-6) are also tested, and the overall maximum likelihood values of the different models are listed in Table 1.We adopt the Akaike information criterion (AIC; Akaike 1974)1 to compare the different models by calculating ∆AIC, the difference between Model 1 and Model 2-6.It is evident from Table 1 that the model including the two uniform disks for SrcA and SrcB and the additional point source SrcC provides the minimum ∆AIC value, which is adopted as the spatial model in the following spectral analysis.Note-In Model 3/6, the centered positions and R68 of SrcA and Src B are same to that in Model 2/5.

Spectral Analysis
With the best-fit spatial model including SrcA/B/C in this region, we compared fits using power-law and log-parabola (LogPb; dN/dE ∝ E −(α+βlog(E/E b )) ) spectral functions to find the best model among the two for SrcA/B in the energy range of 500 MeV -1 TeV.SrcC is spatially coincident with SNR G13.5+0.2, not with HESS J1813-178, and its spectrum is fixed to be a power-law model.The spectral index of SrcC is fitted to be 2.65 ± 0.29, with the integrated  photon flux of (4.43 ± 1.21) × 10 −9 photon cm −2 s −1 .The fitting results with the different spectral models for SrcA/B are listed in Table 2.We compared the overall likelihood values with the PL models for SrcA/B, and the variation suggests an evidence of spectral curvature for the γ-ray spectrum of SrcA with TS curve =2(ln L LogPb−PL -ln L PL−PL ) = 9.0, which corresponds to a significance of ∼3.0σ with one additional degree of freedom.The fitting of the LogPb model for SrcA gives the spectral parameters of α = 2.52 ± 0.10 and β = 0.10 ± 0.05.The integrated photon flux of SrcA is calculated to be (3.71± 0.35) × 10 −8 photon cm −2 s −1 with the TS value of 1358.9.For SrcB, a power-law spectrum is enough to describe its γ-ray emission with the spectral index of 2.11 ± 0.08.The corresponding integrated photon flux is calculated to be (8.72 ± 1.48) × 10 −9 photon cm −2 s −1 with the TS value of 270.6.Note-TScurve is the difference between the TS values obtained for each spectral function.
To obtain the SED of SrcA/B/C, we divide the data from 500 MeV to 1 TeV into 12 bins in energy equally-spaced logarithmically, and perform the same likelihood fitting analysis for each energy bin.For the energy bin with TS value of SrcA/B/C lower than 5.0, an upper limit with 95% confidence level is calculated.The results of the SED are shown in Figure 4, and the global fitting with the PL or LogPb models are also overplotted.The SED of SrcA shows a soft GeV γ-ray spectrum, which is not consistent with the TeV γ-ray SED of HESS J1813-178.While the hard GeV γ-ray spectrum of SrcB smoothly connects with the TeV spectrum of HESS J1813-178.Also considering the spatial coincidence and compact γ-ray emission region between SrcB and HESS J1813-178 shown in Figure 3, we suggest that SrcB is the GeV counterpart of HESS J1813-178.(Albert et al. 2006) are marked by the red, blue and green dots, as shown in the legend.In the middle panel, the magenta and cyan filled butterflies present the spectra of 1LHAASO J1814-1719u observed by WCDA and KM2A of LHAASO, respectively (Cao et al. 2023).In the left and right panels, the magenta solid lines show the hadronic models by adopting a power-law spectrum of protons for SrcA and SrcC, respectively.And the dashed magenta line in the left panel represents the predicted γ-ray emission assuming that the CR spectra therein is the same as that measured locally by AMS-02 (Aguilar et al. 2015).

CO OBSERVATIONS
To search for molecular clouds, spatially and morphologically coincident with the γ-ray emission around HESS J1813-178, we analyzed the 12 CO(J=1-0) and 13 CO(J=1-0) line data observed by the FOREST Unbiased Galactic plane Imaging survey with the Nobeyama 45 m telescope (FUGIN; Umemoto et al. 2017).The Nobeyama telescope has a beam size of 14" at 115 GHz, which results in a postreduction 12 CO(J=1-0) angular resolution of 20" and a 13 CO(J=1-0) angular resolution of 21" for FUGIN survey data with the velocity resolution of 1.3 km s −1 .
The spectra of 12 CO(J=1-0) and 13 CO(J=1-0) emission are extracted toward SrcA within a radius of 0 • .557and 0 • .2, and HESS J1813-178 within a radius of 0 • .206,which are shown in the top panels of Figure 5.The spectra display three distinctive peaks in the velocity intervals of ∼ 10-25, 30-40, and 45-60 km s −1 , and the 12 CO(J = 1-0) intensity maps in the corresponding ranges are shown in the bottom panels of Figure 5.Note that the molecular cloud at 30-40 km s −1 has been studied by Funk et al. (2007) with the NANTEN 12 CO(J=1-0) observations, which is suggested to be related with the star-forming region W33.Although the molecular gas in the velocity ranges of 10-25 and 30-40 km s −1 are partly overlapping with the γ-ray emission, the wide distribution of molecular gas in the range of 45-60 km/s in the γ-ray emission region makes the correlation more likely for the scenario in which the γ-rays are produced by interactions of cosmic rays with the gas.
Adopting the standard Galactic rotation model in Reid et al. (2014), we calculated the possible kinematic distances associated with the velocity range of 45-60 km s −1 .It should be noted that the velocity of each MC could indicate two candidate kinematic distances, the near-side one and the far-side one.The near and far distances of the MC in 45-60 km s −1 are calculated to be 3.9-4.6kpc and 11.5-12.1 kpc, which are in accord with the values of 4.8 ± 0.3 kpc for the young stellar cluster, Cl 1813-178, and the higher value of 12 kpc for PSR J1813-1749/SNR G12.82-0.02.
Adopting the CO-to-H 2 conversion factor X CO = 2 × 10 20 cm −2 (K km s −1 ) −1 (Bolatto et al. 2013), we calculated the total mass of molecular gas within the 68% containment radius of extended GeV γ-ray emission for SrcA (r = 0 • .557)and HESS J1813-178 (r = 0 • .206)with the different distances.For the distance of 4.8 kpc, the total masses of gas within SrcA and HESS J1813-178 are estimated to be ∼ 1.82 × 10 6 d 2 4.8 M ⊙ and ∼ 3.86 × 10 5 d 2 4.8 M ⊙ , respectively.The corresponding average gas number densities are about n gas = 170d −1 4.8 cm −3 and 720d −1 4.8 cm −3 by assuming a spherical geometry of the gas distribution.For the distance of 12 kpc, the calculated total masses of gas within SrcA and HESS J1813-178 are ∼ 1.14 × 10 7 d 2 12 M ⊙ and ∼ 2.41 × 10 6 d 2 12 M ⊙ with the corresponding average gas number densities of 70d −1 12 cm −3 and 300d −1 12 cm −3 .The spatial and spectral data analyses with the different energy ranges reveal that the diffuse GeV γ-ray emission around HESS J1813-178 can be distinguished into three separate components.SrcA with an extension of 0 • .557has a soft GeV γ-ray spectrum, which has no identified associated counterpart, and a point source in the north of HESS J1813-178, SrcC, is spatially associated with the SNR G13.5+0.2.For SrcB, the GeV γ-ray emission region is more compact than SrcA, which is spatially coincident with the TeV γ-rays of HESS J183-178.The GeV γ-ray spectrum of SrcB is hard with an index of 2.11, which also smoothly connects with the TeV spectrum of HESS J1813-178.Both the spatial morphology of the γ-ray emission and the GeV-TeV γ-ray spectrum suggest that SrcB could be the GeV counterpart of HESS J1813-178.Moreover, the CO observations display that the intensity of the γ-ray emission from SrcA and HESS J1813-178 spatially correlate with the molecular gas distribution in the velocity range of 45-60 km s −1 .Such velocity range indicates two candidate kinematic distances, and the near one is in accord with the distance of the young stellar cluster, Cl 1813-178, of ∼ 4.8 kpc.The far one with ∼ 12 kpc is in accord with the distance of PSR J1813-1749/SNR G12.82-0.02system.

The possible γ-ray origin of HESS J1813-178
Considering the association between HESS J1813-178 and PSR J1813-1749/SNR G12.82-0.02 or Cl 1813-178 at the different distances, here we discuss these two possibilities as the origin of HESS J1813-178.The γ-ray emission can be produced by the decay of neutral pions due to the inelastic nuclei-nuclei collisions, which is the hadronic process.In addition, the leptonic process, that is the inverse Compton scattering or the bremsstrahlung (Brem) process of relativistic electrons, can also produce the intense γ-ray emission.For the model fitting, the TeV observations of HESS J1813-178 from H.E.S.S. (Aharonian et al. 2006;H. E. S. S. Collaboration et al. 2018a) and MAGIC (Albert et al. 2006) are adopted.The observations by LHAASO-WCDA show a larger extension and a higher flux of the γ-ray emission around HESS J1813-178 in the energy range of 1-25 TeV (Cao et al. 2023), which are not consistent with the results of H.E.S.S. and MAGIC.Such results could be due to the contributions from the unresolved γ-ray sources in the more extended γ-ray region of LHAASO-WCDA.LHAASO-KM2A gives an upper limit of the extension of the γ-ray emission above 25 TeV, and its spectrum also could connect with that of H.E.S.S. and MAGIC.Therefore, only LHAASO-KM2A data are used in the following model fitting.Funk et al. (2007).The gray and magenta butterflies mark the X-ray spectra of PWN observed by Funk et al. (2007) and Ho et al. (2020), together with the INTEGRAL observations shown as the purple dots (Brogan et al. 2005).For the left and middle panels, the green and magenta dashed lines represent the synchrotron and bremsstrahlung components, respectively.The contributions from ICS process with the different radiation fields are shown as the dotted lines, and the black solid line represents the sum of the contributions from different radiation components.For the right panel, the radiation components of zone 1 and zone 2 are marked by the dotted and dot-dashed lines, as shown in the legend.The contributions from ICS-OPT and bremsstrahlung components of zone 1 are too low to show.

PWN or SNR scenario at the distance of 12 kpc
The PWN component associated with PSR J1813-1749 was first resolved by Helfand et al. (2007) using the observations of Chandra.The spectrum of the PWN follows a power law with an index of ∼1.3, and an absorbed flux of 5.6 × 10 −12 erg cm −2 s −1 is obtained in the energy range of 2-10 keV.Following XMM-Newton observations also resolved the PWN with an spectral index of ∼1.8 and an absorbed flux of 7 × 10 −12 erg cm −2 s −1 from 2 to 10 keV (Funk et al. 2007).In the high energy band of 20-100 keV, INTEGRAL discovered a soft γ-ray source with an index of ∼1.8, which is spatially consistent with the PWN component (Brogan et al. 2005).The newest X-ray data analysis combining Chandra and NICER data by Ho et al. (2020) modeled the PWN with an spectral index of 1.5 ± 0.1, which is slight harder than the results of previous works.And an absorbed 2-10 keV flux of 7.6 × 10 −12 erg cm −2 s −1 is calculated.
Adopting the updated γ-ray spectrum and the new X-ray observations for PWN, we first used a one-zone leptonic model to fit the multi-wavelength data of HESS J1813-178.The radio upper limits toward the PWN region by VLA was also considered (Funk et al. 2007).For the ICS process in the leptonic model, the radiation field includes the cosmic microwave background (CMB), the infrared (IR) component from interstellar dust and gas with the typical values of T = 30K & u = 1 eV cm −3 , and the optical (OPT) component from stellar sources with T = 3000K & u = 1 eV cm −3 (Porter et al. 2006(Porter et al. , 2008)).For the leptonic model, the spectra of electrons was assumed to be a broken power-law model with an exponential cutoff (BPL) for typical PWNe (Bucciantini et al. 2011), which follows: Here, α e,1 and α e,2 are the spectral indices below and above the break energy of E e,break .The cutoff energy of electrons is denoted by E e,cut .The distance of HESS J1813-178 as a PWN is adopted to be 12.0 kpc (Camilo et al. 2021).The CO observations in Section 3 give an average gas density of 300d −1 12 cm −3 toward HESS J1813-178.Such density value would produce a quite high flux of bremsstrahlung component, which is much higher than that of ICS component considering the same electronic distribution.Therefore, to reduce the contribution from bremsstrahlung process, we also adopted a low density value of 1.0 cm −3 by assuming that the similarity of the pulsar distance to the molecular distance from the observer is a coincidence.The naima package was used to fit the multi-wavelength data with the Markov Chain Monte Carlo (MCMC) algorithm (Zabalza 2015).
The best-fit parameters of one-zone leptonic models for PWN are shown in Table 3, and the corresponding modeled SEDs are given in the left panel of Figure 6.For the ICS-dominated model with the gas density of 1.0 cm −3 , the break energy of electrons is fitted to be ∼140 GeV.The spectral indices below and above it are 1.12 and 2.78, respectively.The cutoff energy of ∼700 TeV is obtained, with the total energy of electrons above 1 GeV of ∼10 49 erg.A magnetic field strength of about 7.3 µG is needed to explain the flux in the X-ray band.Nonetheless, the hard X-ray spectrum given by Ho et al. (2020) suggests a harder electronic spectrum above the break energy, which can not be in accord with the flat γ-ray spectrum.
With the gas density of 300 cm −3 , the γ-ray spectrum of HESS J1813-178 is dominated by the bremsstrahlung process, as shown in the middle panel of Figure 6.And the spectral indices of electrons below and above the break energy of ∼125 GeV are fitted to be 1.98 and 2.33, together with the cutoff energy of ∼400 TeV.The magnetic field strength is similar to that of the one-zone ICS-dominated model, while the total energy of electrons is about one magnitude lower than that.Although the synchrotron component of the Brem-dominated model could better accord with the radio upper limits and the hard X-ray spectrum without a low-energy cutoff of electrons, the curved γ-ray spectrum produced by the bremsstrahlung process can not match the γ-ray observations well, particularly at the lowest energies.
Given the hard X-ray spectrum and the flat GeV γ-ray spectrum, we also considered a two-zone leptonic model: one zone (zone 1) contains the newly accelerated electrons with the highest energies, which produces the X-ray emission.Another zone (zone 2) contains the old electrons accumulated in the region of GeV γ-ray emission.Such two-zone leptonic model could also explain the smaller extension of X-ray emission compared with the GeV or TeV γ-ray emission.For the electronic distribution of each zone, a single power-law with an exponential cutoff (PL) model is adopted in the form of where α e,i and E e,cut,i (i = 1 or 2) are the spectral index and the cutoff energy of electrons for each zone.The low density value of 1.0 cm −3 is adopted for this model.The best-fit parameters of two-zone leptonic models for PWN are shown in Table 3, with the modeled SED shown in the right panel of Figure 6.For zone 1, the spectral index and cutoff energy of electrons are fitted to be 1.64 and ∼200 TeV to explain the hard X-ray and TeV γ-ray spectra.A high magnetic field strength of ∼8.5 µG is needed.And the total energy of electrons for zone 1 is about 7.85×10 46 erg.The GeV γ-ray emission could be explained by the electrons in zone 2 with an index of 2.77.The fitted cutoff energy of electrons of ∼90 TeV is lower than that of zone 1, while total energy of electrons of 2.97×10 50 erg in zone 2 is much higher than that of zone 1.The magnetic field strength in zone 2 is much lower with ∼1.5 µG.And the much lower upper limits in the radio band need a low-energy cutoff of ∼50 GeV for the electronic spectrum in zone 2. The soft GeV spectrum of SrcA suggests that its γ-ray emission could be related to the escaped electrons from the PWN associated with PSR J1813-1749, like Vela X (Hinton et al. 2011;Grondin et al. 2013).The steep GeV spectrum of SrcA indicates an absence of ≳100 GeV energy particles from the extended γ-ray emission region.Considering a typical magnitude of magnetic field with ∼ µG, the synchrotron cooling time-scale for particles emitting in the GeV range is up to ∼ 10 8 yr, which makes the synchrotron cooling effect to be implausible for the steep GeV spectrum of SrcA.And the possible scenario is that the bulk of the ≳100 GeV particles have escaped from the GeV emission region into the interstellar medium (Hinton et al. 2011).The required diffusion coefficient is calculated to be D = r 2 /4t ≃ 2×10 29 cm 2 s −1 , where the radius of GeV γ-ray emission region of SrcA is estimated to be r≃120pc, and the diffusion time is adopted to be the characteristic age of PSR J1813-1749 with t ∼ 5600 yr.This value is similar to the Galactic diffusion coefficient, but much faster than Bohm diffusion for the magnetic field of ∼ µG in the γ-ray emission region of SrcA.Now, the scenario that the γ-ray emission of HESS J1813-178 is associated with SNR G12.82-0.02 is considered.The γ-ray SNRs could be roughly divided into two classes with different spectra (Funk 2015;Zeng et al. 2019): the young-aged SNRs, including RX J1713.7-3946(H.E. S. S. Collaboration et al. 2018b) and RX J0852.0-4622(H.E. S. S. Collaboration et al. 2018c), have relatively harder GeV γ-ray spectra, which are suggested to be from the leptonic process.The old-aged SNRs interacting with molecular clouds, like IC 443 and W44 (Ackermann et al. 2013), have softer γ-ray spectra with a break at ∼ GeV, whose γ-ray emission are suggested to be from the hadronic process.Considering the spatial association between the γ-ray emission from HESS J1813-178 and the molecular gas distribution, the hadronic model is adopted.The spectrum of protons is assumed to be a single power-law with an exponential cutoff in the form of where α p and E p,cut are the spectral index and the cutoff energy of protons, respectively.β describes the sharpness of the cutoff, and the typical values of 0.5, 1.0, and 2.0 are adopted to constrain the parameters in the model.The best-fit parameters of hadronic models for HESS J1813-178 are shown in Table 4, and the corresponding modeled γ-ray spectrum is given in Figure 7.The hadronic models with the three different values of β can reproduce the γ-ray spectrum with little differences.The spectral indices of protons are fitted to be about 2.1, while the cutoff energy varies from 261 TeV for β=0.5 to 737 TeV for β=1.0.The total energy of protons above 1 GeV, W p , is estimated to be about (1.26-1.51)× 10 49 (n gas /300 cm −3 ) −1 (d/12.0kpc) 2 erg, corresponding to ∼1% particle acceleration efficiency for the explosion energy of a typical supernova (SN; E SN ∼ 10 51 erg).Brogan et al. (2005) limited the age of SNR G12.82-0.02from 285 to 2520 yrs with the assumptions of the distance of 4 kpc and the ambient density of 1.0 cm −3 .Considering the far distance of 12 kpc and the high density associated with the molecular cloud around the remnant, the age of G12.82-0.02should be much older.Such old-aged SNR would not expect to accelerate the needed protons with energy of hundreds of TeV (Reynolds 2008).However, considering the much larger extension of the GeV γ-ray emission and the clouds compared with the radio size of SNR G12.82-0.02,cosmic rays that escaped from the SNR in earlier evolutional epochs would be the possible explanation.Such scenario is similar to a proposed origin of the γ-ray emission around SNR DA 530 (Xin & Guo 2023) and the SNR associated with PSR J0837-2454 (Zhang & Xin 2023).Given the typical value of Galactic diffusion coefficient, D(E) = 3 × 10 28 (E/10GeV) δ (δ = 1/3 or 1/2 for Kolmogorve or Kraichnan turbulence; Blasi 2013), and the radius of the GeV γ-ray emission region of r = 43 pc, the diffusion time for protons with energy of 100 GeV is calculated to be t diff = r 2 /4D(E) ≃ 500-2000 yr for different values of δ.Such time is still lower than the age of SNR G12.82-0.02,which is reasonable for such explanation.And the molecular clouds can accumulate those cosmic rays in the past since the diffusion coefficient is expected to be low in dense environments (Aharonian & Atoyan 1996).

Young Stellar Cluster scenario at the distance of 4.8 kpc
Over the last decade, more and more observations support the young stellar clusters associated with star forming regions to be an important class of factories of Galactic cosmic rays (Aharonian et al. 2019).Cosmic rays could be efficiently accelerated by strong fast winds of young massive stars and shocks caused by the core-collapse SNs in clusters.The acceleration efficiency and maximum energy of particles are expected to be enhanced with respect to the standard values derived for single SNR shock, considering the fact that the SN blast wave could interact with fast stellar winds (Bykov et al. 2020).Meanwhile, YSCs typically host the dense molecular gas to drive the strong star formation, which makes the hadronic interaction of accelerated CRs with the surrounding dense gas a natural explanation for their γ-ray emission.So far, several such systems, e.g.Cygnus cocoon associated with the compact cluster Cygnus OB2 (Ackermann et al. 2011;Aharonian et al. 2019), Westerlund 1 (Abramowski et al. 2012), Westerlund 2 (Yang et al. 2018), NGC 3603 (Yang & Aharonian 2017;Saha et al. 2020), and 30 Dor C (H. E. S. S. Collaboration et al. 2015), have been detected with the extended γ-ray structures from GeV to TeV bands.Especially, LHAASO recently detected the ultra high energy (UHE; >100TeV) γ-ray emission toward Cygnus cocoon in the energy range extending to few PeV (LHAASO Collaboration 2023).The γ-ray spectra of YSCs are relatively hard with the typical index of 2.1-2.3 (Aharonian et al. 2019).
The association between the γ-ray emission from HESS J1813-178 and YSC Cl 1813-178 has been discussed in the previous works (Messineo et al. 2011;Araya 2018).With the updated GeV γ-ray spectrum, the YSC scenario with the hadronic model is also considered here.The fitting results shown in Figure 7 is same to that of the SNR scenario.However, for the YSC scenario, the total energy of protons above 1 GeV, W p , is calculated to be (0.84-1.01) × 10 48 (n gas /720 cm −3 ) −1 (d/4.8 kpc) 2 erg.Here, the gas density of 720 cm −3 calculated with the distance of 4.8 kpc for Cl 1813-178 is adopted.The spectral index of proton distribution is similar to that of the typical YSCs, while the total energy for HESS J1813-178/Cl 1813-178 is about one order of magnitude lower than that of other YSCs (Aharonian et al. 2019).
The near-infrared spectroscopic survey in the direction of Cl 1813-178 cluster confirmed 25 massive stars, including 2 Wolf-Rayet (WR) stars, a candidate luminous blue variable (cLBV) and 21 OB stars (Messineo et al. 2011).Different from the other known LBVs in clusters, which are typically among the most luminous members, the cLBV in Cl 1813-178 appears to be rather faint and has a mass smaller than that of the two detected WR stars.Given the typical wind powers expected for a 35 M ⊙ star on and off the main sequence and in the WR stage (Parizot et al. 2004) and ignoring the contribution from the cLBV, we get the total kinetic wind power of Cl 1813-178 with ∼ 10 38 erg s −1 .The kinetic power is similar to that of Cygnus Cocoon with 2×10 38 erg s −1 (Ackermann et al. 2011), but much lower than that of Westerlund 1 with 1×10 39 erg s −1 (Aharonian et al. 2019) and Westerlund 2 with 2×10 40 erg s −1 (Rauw et al. 2007).With the distance of 4.8 kpc, the GeV γ-ray luminosity of HESS J1813-178 is calculated to be about 9.6×10 34 (d/4.8kpc) 2 erg s −1 , which represents ∼0.1% of the stellar wind power in Cl 1813-178.The conversion efficiency is slightly higher than that of Cygnus Cocoon with ∼0.03% (Ackermann et al. 2011), Westerlund 1 with ∼0.002-0.02%(Härer et al. 2023) and Westerlund 2 with ∼0.005% (Yang et al. 2018).However, the GeV extension of HESS J1813-178 with r = 17 pc for d = 4.8 kpc is much smaller that that of other YSCs (e.g.r = 50 pc for Cygnus Cocoon, r = 60 pc for Westerlund 1, and r = 210 pc for Westerlund 2; Aharonian et al. 2019;Yang et al. 2018).If the GeV extension of r = 17 pc reflects the propagation depth of CRs, the diffusion coefficient is calculated to be D = r 2 /4T ≃ 5.5×10 24 cm 2 s −1 , where the age of Cl 1813-178 cluster is adopted to be T = 4×10 6 yr (Messineo et al. 2011).While for the interstellar magnetic field B ∼ 10 µG, the Bohm diffusion coefficient is D B = R L c/3 ≃ 3 × 10 24 (E/1TeV)(B/10µG) −1 cm 2 s −1 .And for the 1 TeV protons responsible for 100 GeV γ-rays, the diffusion coefficient in HESS J1813-178 is close to Bohm limit, which seems unrealistic.Such result implies that the radius of the CR halo around Cl 1813-178 is most likely significantly larger than that of the γ-ray emission region, like Westerlund 1 (Aharonian et al. 2019).
Considering the molecular gas around SrcA, we also used the hadronic model to fit the γ-ray spectrum of SrcA, by assuming a power-law distribution for protons.The modeled spectrum is shown in Figure 1, together with the predicted γ-ray emission assuming that the CR spectra therein are the same as those measured locally by AMS-02 (Aguilar et al. 2015).The best-fit parameters of hadronic models are listed in Table 4.For SrcA, the spectral of protons with an index of ∼ 2.71 is need to explain its soft γ-ray spectrum detected.With the different values of gas density derived from the different distances for SrcA, the total energy of protons above 1 GeV is estimated to be W p ∼ 3.76 × 10 50 (n gas /70 cm −3 ) −2 (d/12.0kpc) 2 erg or ∼ 2.4 × 10 49 (n gas /170 cm −3 ) −2 (d/4.8 kpc) 2 erg.The extended γ-ray morphology and the association with the molecular gas for SrcA suggest that the high energy protons could be accelerated and escaped from YSC Cl 1813-178 associated with HESS J1813-178.Because of the energy-dependence of the diffusion coefficient, the initial particle spectrum of dN/dE ∝ E −α would be modified to be dN/dE ∝ E −(α+δ) for the escaped particle spectrum with the continuous injection (Aharonian et al. 2019).With the values of α ∼ 2.1 and δ = 1/2, the escaped proton spectrum follows dN/dE ∝ E −2.6 , which is close to the fitted proton spectrum above.The radius of the γ-ray emission region for SrcA is about 47 pc with the distance of 4.8 kpc.If the escaping scenario between SrcA and YSC Cl 1813-178 is real, the diffusion length of 47 pc constrains that the diffusion coefficient should be as low as 10 25 cm 2 s −1 by equating the diffusion time to be the age of cluster, which is also close to Bohm limit.
SrcC is spatially coincident with SNR G13.5+0.2.And the hadronic model is also used to fit its γ-ray spectrum.By adopting the distance of 13.0 kpc (Lee et al. 2020) and assuming the gas density of 10 cm −3 , the total energy of protons for SNR G13.5+0.2 is calculated to be ∼ 2.61 × 10 50 (n gas /10 cm −3 ) −2 (d/13.0kpc) 2 erg, which means that ∼26% of the SN kinetic energy of ∼10 51 erg is transferred to the energy of particles.

CONCLUSIONS
In this work, we analyzed the GeV γ-ray emission in the field of HESS J1813-178, using 14 years of Fermi-LAT data, and found that the GeV emission around HESS J1813-178 is composed of three γ-ray components: SrcA, SrcB and SrcC.SrcC is spatially coincident with SNR G13.5+0.2.SrcB has a small extension, which is spatially coincident with HESS J1813-178.The GeV γ-ray spectrum follows a power-law model with an index of 2.11 ± 0.08, which can connect with the TeV SED of HESS J1813-178 smoothly.The spatial and spectral associations indicate that SrcB could be the GeV counterpart of HESS J1813-178.SrcA with a large size extension has a soft GeV γ-ray spectrum described by a log-parabola model.No clear identified counterparts in other wavelengths are found in the region of SrcA.While the CO observations display that the intensity of the γ-ray emission from SrcA and HESS J1813-178 spatially correlate with the molecular gas distribution in the velocity range of 45-60 km s −1 .The velocity range indicates two candidate kinematic distances: the near one is in accord with the distance of the young stellar cluster, Cl 1813-178, of ∼4.8 kpc, and the far one is in accord with the distance to the PSR J1813-1749/SNR G12.82-0.02system, of ∼12 kpc.
We then discussed the possible γ-ray origin of HESS J1813-178 with three scenarios: the PWN associated with PSR J1813-1749, SNR G12.82-0.02 or YSC Cl 1813-178.For the one-zone leptonic scenario with PWN, the electron distribution is adopted to be a broken power-law model with an exponential cutoff.The one-zone ICS-dominated and bremsstrahlung-dominated models with the different values of gas density are considered.However, none of them can well describe the multi-wavelength data from radio, X-ray to γ-ray bands.The one-zone ICS-dominated model expects a flat X-ray spectrum, which can not match the observed one with an spectral index of 1.5 ± 0.1 in the energy range of 2-10 keV.With the high density around HESS J1813-178, the one-zone bremsstrahlung-dominated model could better accord with the hard X-ray spectrum.While it predicts a curved γ-ray spectrum dominated by the bremsstrahlung process, which can not well match the γ-ray observations.A two-zone leptonic model for PWN could explain the multi-wavelength data.The electron distribution of each zone is adopted to be a single power-law with an exponential cutoff.In one zone, an electron distribution with an index of ∼1.64 could explain the hard X-ray and TeV γ-ray spectra, with a high magnetic field strength of ∼8.5µG.In another zone, an electron distribution with an index of ∼2.77 and a low-energy cutoff of ∼50 GeV is use to explain the radio upper limits and the GeV γ-ray emission, with a low magnetic field strength of ∼1.5µG.
For the SNR scenario, only the hadronic model is discussed, considering the association between the γ-ray emission of HESS J1813-178 and the surrounding molecular gas.The spectra of protons was assumed to be a single power-law with an exponential cutoff, and the different values of the sharpness of the cutoff, β, are also adopted.The hadronic models with the different values of β can reproduce the γ-ray spectra with little differences, and the proton distribution needs to be hard with an index of ∼2.1.The cutoff energy varies between 261 TeV and 737 TeV for the different values of β.The large extension of the GeV γ-ray emission and the clouds compared with the size of SNR G12.82-0.02suggests that these protons contributed to the γ-ray emission could escape from the SNR in earlier evolutional epochs.The calculated diffusion time is lower than the age of SNR G12.82-0.02,which is reasonable for such explanation.
Although the flat γ-ray spectrum of HESS J1813-178 is similar to that of YSCs with the typical spectral index of 2.1-2.3, the GeV extension of HESS J1813-178 is much smaller than that of other YSCs.Considering the GeV γ-ray size as the propagation depth of CRs, the calculated diffusion coefficient needs to be close to Bohm limit, which is not very reasonable.
The soft γ-ray spectrum of SrcA can be explained by the fast diffusion of electrons from the PWN associated with PSR J1813-1749, like Vela X, or the hadronic process due to the inelastic collisions between the molecular cloud and the high energy protons, which can be accelerated in and escaped from YSC Cl 1813-178.
It should be noted that the hadronic models show that the highest energy of protons should exceed a few hundred TeV, which makes HESS J1813-178 to be a promising PeV CR accelerator (PeVatron).LHAASO has detected the UHE γ-ray around HESS J1813-178, and more detailed analysis based on more LHAASO data, especially the morphological analysis, would be helpful to investigate the origin of the γ-ray emission in this region and test its PeVatron nature.Moreover, the sensitivity of LHAASO above 10 TeV could help to study in detail the shape of the particle spectrum cutoff, that encodes information about the particle acceleration and transport mechanisms (Zirakashvili & Aharonian 2007;Romoli et al. 2017).
We would like to thank the anonymous referee for very helpful comments, which help to improve the paper.This work is based on observations made with NASA's Fermi Gamma-Ray Space Telescope, and it also makes use of molecular line data from FOREST Unbiased Galactic plane Imaging survey with the Nobeyama 45 m telescope (FUGIN).This work is supported by the National Natural Science Foundation of China under the grants 12103040.

Figure 1 .
Figure 1.The SED of 4FGL J1813.1-1737e with black dots.The arrows indicate the 95% upper limits, and the gray histogram denotes the TS value for each energy bin.The cyan dots are the results from Araya (2018).

Figure 2
Figure 2. 2.5 • × 2.5 • TS maps in the energy range of 500 MeV -10 GeV with "PSF3" events (left) and 10 GeV -1 TeV with "Source" events (right).The white and cyan circles represent the spatial size (R68) of SrcA and SrcB in the two energy bands.The positions of 4FGL-DR4 sources are shown as the red crosses.In the left panel, the magenta dashed circle shows R68 of 4FGL J1813.1-1737e in 4FGL-DR4 catalog.In the right panel, the white cross shows the best-fit position of SrcC, and the associated SNR G13.5+0.2 is marked by the red circle.The yellow circle represents R68 of 1LHAASO J1814-1719u observed by WCDA of LHAASO, and the upper limit of extension given by KM2A above 25 TeV is shown as the green circle(Cao et al. 2023).

Figure 3 .
Figure 3. Zoom-in of the TS map for a region of 1.0 • × 1.0 • with the energy of 10 GeV -1 TeV.The cyan circle shows R68 of SrcB.The upper limit size of 1LHAASO J1814-1719 observed by KM2A is marked by the green circle.The black solid and dashed circles show the radio sizes of SNR G12.82-0.02 and SNR G12.72-0.00,respectively.The region of the stellar cluster Cl 1813-178 is marked by the magenta circle (Messineo et al. 2011).The position of PSR J1813-1749 is shown as the blue cross (Gotthelf & Halpern 2009).The centroid position of the TeV γ-ray emission from HESS J1813-178 detected by MAGIC (Albert et al. 2006) is marked by the green diamond, and the white contours represent the TeV γ-ray emission detected by H.E.S.S. (H.E. S. S. Collaboration et al. 2018a).

Figure 4 .
Figure 4.The SEDs of SrcA (left), SrcB (middle) and SrcC (right).The black dots depict the results of Fermi-LAT data in the energy range of 500 MeV -1 TeV, together with the global best-fit spectra with 1σ propagated statistical errors shown as dashed gray lines.The black arrows indicate the 95% upper limits and the gray histogram shows the TS value for each energy bin.The TeV γ-ray data of HESS J1813-178 observed by H.E.S.S. (Aharonian et al. 2006; H. E. S. S. Collaboration et al. 2018a) and MAGIC(Albert et al. 2006) are marked by the red, blue and green dots, as shown in the legend.In the middle panel, the magenta and cyan filled butterflies present the spectra of 1LHAASO J1814-1719u observed by WCDA and KM2A of LHAASO, respectively(Cao et al. 2023).In the left and right panels, the magenta solid lines show the hadronic models by adopting a power-law spectrum of protons for SrcA and SrcC, respectively.And the dashed magenta line in the left panel represents the predicted γ-ray emission assuming that the CR spectra therein is the same as that measured locally by AMS-02(Aguilar et al. 2015).

Figure 6 .
Figure6.The leptonic model fitting as a PWN for HESS J1813-178 with one-zone ICS-dominated scenario (left), Bremdominated scenario (middle) and two-zone ICS-dominated scenario (right).The radio upper limits shown as the cyan arrows are fromFunk et al. (2007).The gray and magenta butterflies mark the X-ray spectra of PWN observed byFunk et al. (2007) andHo et al. (2020), together with the INTEGRAL observations shown as the purple dots(Brogan et al. 2005).For the left and middle panels, the green and magenta dashed lines represent the synchrotron and bremsstrahlung components, respectively.The contributions from ICS process with the different radiation fields are shown as the dotted lines, and the black solid line represents the sum of the contributions from different radiation components.For the right panel, the radiation components of zone 1 and zone 2 are marked by the dotted and dot-dashed lines, as shown in the legend.The contributions from ICS-OPT and bremsstrahlung components of zone 1 are too low to show.

Figure 7 .
Figure 7.The hadronic model fitting for HESS J1813-178 with the different values of the sharpness of the cutoff, β.The gray square data show the γ-ray spectrum of Cygnus Cocoon detected by Fermi-LAT and LHAASO (Yang & Aharonian 2017; LHAASO Collaboration 2023).

Table 3 .
Parameters of the leptonic model as a PWN for HESS J1813-178 Note-αe,1 and αe,2 represent the indices below and above E e,break for BPL spectrum of electrons in one-zone model, or the indices for two PL spectra of electrons in two-zone model.And the total energy of electrons, We, is calculated for Ee > 1 GeV.