# Unveiling extended gamma-ray emission around HESS J1813-178

F Aharonian, F. Ait Benkhali, J Aschersleben, H Ashkar, M Backes, A Baktash, V. Barbosa Martins, J Barnard, R Batzofin, Y Becherini, et al.

► To cite this version:

F Aharonian, F. Ait Benkhali, J Aschersleben, H Ashkar, M Backes, et al. Unveiling extended gamma-ray emission around HESS J1813-178. Astronomy & Astrophysics - A&A, 2024, 686, pp.A149. 10.1051/0004-6361/202348374. [hal-04537734](https://hal.science/hal-04537734v1)

HAL Id: hal-04537734

https://hal.science/hal-04537734v1

Submitted on 27 Sep 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Unveiling extended gamma-ray emission around HESS J1813-178

H.E.S.S. Collaboration: F. Aharonian<sup>1,2,4</sup>, F. Ait Benkhali<sup>5</sup>, J. Aschersleben<sup>6</sup>, H. Ashkar<sup>7</sup>, M. Backes<sup>8,9</sup>, A. Baktash<sup>10</sup>, V. Barbosa Martins<sup>11</sup>, J. Barnard<sup>12</sup>, R. Batzofin<sup>13</sup>, Y. Becherini<sup>14,15</sup>, D. Berge<sup>11,16</sup>, K. Bernlöhr<sup>2</sup>, B. Bi<sup>17</sup>, M. Böttcher<sup>9</sup>, C. Boisson<sup>18</sup>, J. Bolmont<sup>19</sup>, M. de Bony de Lavergne<sup>20</sup>, J. Borowska<sup>16</sup>, M. Bouyahiaoui<sup>2</sup>, M. Breuhaus<sup>2</sup>, R. Brose<sup>1</sup>, F. Brun<sup>20</sup>, B. Bruno<sup>21</sup>, T. Bulik<sup>22</sup>, C. Burger-Scheidlin<sup>1</sup>, S. Caroff<sup>23</sup>, S. Casanova<sup>24</sup>, R. Cecil<sup>10</sup>, J. Celic<sup>21</sup>, M. Cerruti<sup>14</sup>, P. Chambery<sup>25,*</sup>, T. Chand<sup>9</sup>, A. Chen<sup>26</sup>, J. Chibueze<sup>9</sup>, O. Chibueze<sup>9</sup>, G. Cotter<sup>27</sup>, J. Damascene Mbarubucyeye<sup>11</sup>, A. Djannati-Ataï<sup>14</sup>, A. Dmytriiev<sup>9</sup>, V. Doroshenko<sup>17</sup>, S. Einecke<sup>28</sup>, J.-P. Ernenwein<sup>29</sup>, K. Feijen<sup>14</sup>, M. Filipovic<sup>40</sup>, G. Fontaine<sup>7</sup>, M. Füßling<sup>11</sup>, S. Funk<sup>21</sup>, S. Gabici<sup>14</sup>, Y. A. Gallant<sup>30</sup>, S. Ghafourizadeh<sup>5</sup>, G. Giavitto<sup>11</sup>, D. Glawion<sup>21</sup>, J. F. Glicenstein<sup>20</sup>, P. Goswami<sup>14</sup>, G. Grolleron<sup>19</sup>, M.-H. Grondin<sup>25</sup>, J. A. Hinton<sup>2</sup>, W. Hofmann<sup>2</sup>, T. L. Holch<sup>11</sup>, M. Holler<sup>31</sup>, M. Jamrozy<sup>32</sup>, F. Jankowsky<sup>5</sup>, V. Joshi<sup>21,*</sup>, I. Jung-Richardt<sup>21</sup>, K. Katarzyński<sup>33</sup>, R. Khatoon<sup>9</sup>, B. Khélifi<sup>14</sup>, S. Klepser<sup>11</sup>, W. Kluźniak<sup>34</sup>, Nu. Komin<sup>26</sup>, K. Kosack<sup>20</sup>, D. Kostunin<sup>11</sup>, A. Kundu<sup>9</sup>, R. G. Lang<sup>21</sup>, S. Le Stum<sup>29</sup>, F. Leitl<sup>21</sup>, A. Lemière<sup>14</sup>, M. Lemoine-Goumard<sup>25</sup>, J.-P. Lenain<sup>19</sup>, F. Leuschner<sup>17</sup>, J. Mackey<sup>1</sup>, D. Malyshev<sup>17</sup>, D. Malyshev<sup>21</sup>, V. Marandon<sup>20</sup>, P. Marinos<sup>28</sup>, G. Martí-Devesa<sup>31</sup>, R. Marx<sup>5</sup>, A. Mehta<sup>11</sup>, A. Mitchell<sup>21,*</sup>, R. Moderski<sup>34</sup>, L. Mohrmann<sup>2</sup>, A. Montanari<sup>5</sup>, E. Moulin<sup>20</sup>, T. Murach<sup>11</sup>, M. de Naurois<sup>7</sup>, J. Niemiec<sup>24</sup>, A. Priyana Noel<sup>32</sup>, P. O'Brien<sup>35</sup>, S. Ohm<sup>11</sup>, L. Olivera-Nieto<sup>2</sup>, E. de Ona Wilhelmi<sup>11</sup>, M. Ostrowski<sup>32</sup>, S. Panny<sup>31</sup>, M. Panter<sup>2</sup>, R. D. Parsons<sup>16</sup>, D. A. Prokhorov<sup>36</sup>, G. Pühlhofer<sup>17</sup>, M. Punch<sup>14</sup>, A. Quirrenbach<sup>5</sup>, M. Regeard<sup>14</sup>, P. Reichherzer<sup>20</sup>, A. Reimer<sup>31</sup>, O. Reimer<sup>31</sup>, H. Ren<sup>2</sup>, M. Renaud<sup>30</sup>, B. Reville<sup>2</sup>, F. Rieger<sup>2</sup>, G. Roellinghoff<sup>21</sup>, B. Rudak<sup>34</sup>, V. Sahakian<sup>3</sup>, H. Salzmann<sup>17</sup>, M. Sasaki<sup>21</sup>, F. Schüssler<sup>20</sup>, H. M. Schutte<sup>9</sup>, J. N. S. Shapopi<sup>8</sup>, A. Specovius<sup>21</sup>, S. Spencer<sup>21</sup>, R. Steenkamp<sup>8</sup>, S. Steinmassl<sup>2</sup>, C. Steppa<sup>13</sup>, I. Sushch<sup>9</sup>, H. Suzuki<sup>37</sup>, T. Takahashi<sup>38</sup>, T. Tanaka<sup>37</sup>, R. Terrier<sup>14</sup>, M. Tluczykont<sup>10</sup>, N. Tsuji<sup>39</sup>, T. Unbehaun<sup>21</sup>, C. van Eldik<sup>21</sup>, M. Vecchi<sup>6</sup>, J. Veh<sup>21</sup>, C. Venter<sup>9</sup>, J. Vink<sup>36</sup>, T. Wach<sup>21,*</sup>, S. J. Wagner<sup>5</sup>, A. Wierzcholska<sup>24</sup>, M. Zacharias<sup>5,9</sup>, D. Zargaryan<sup>1</sup>, A. A. Zdziarski<sup>34</sup>, S. Zouari<sup>14</sup>, and N. Żywucka<sup>9</sup>

*(Affiliations can be found after the references)*

Received 24 October 2023 / Accepted 21 March 2024

## ABSTRACT

**Context.** HESS J1813-178 is a very-high-energy γ-ray source spatially coincident with the young and energetic pulsar PSR J1813-1749 and thought to be associated with its pulsar wind nebula (PWN). Recently, evidence for extended high-energy emission in the vicinity of the pulsar has been revealed in the *Fermi* Large Area Telescope (LAT) data. This motivates revisiting the HESS J1813-178 region, taking advantage of improved analysis methods and an extended dataset.

**Aims.** Using data taken by the High Energy Stereoscopic System (H.E.S.S.) experiment and the *Fermi*-LAT, we aim to describe the γ-ray emission in the region with a consistent model, to provide insights into its origin.

**Methods.** We performed a likelihood-based analysis on 32 hours of H.E.S.S. data and 12 yr of *Fermi*-LAT data and we fitted a spectro-morphological model to the combined datasets. These results allowed us to develop a physical model for the origin of the observed γ-ray emission in the region.

**Results.** In addition to the compact very-high-energy γ-ray emission centred on the pulsar, we find a significant yet previously undetected component along the Galactic plane. With *Fermi*-LAT data, we confirm extended high-energy emission consistent with the position and elongation of the extended emission observed with H.E.S.S. These results establish a consistent description of the emission in the region from GeV energies to several tens of TeV.

**Conclusions.** This study suggests that HESS J1813-178 is associated with a γ-ray PWN powered by PSR J1813-1749. A possible origin of the extended emission component is inverse Compton emission from electrons and positrons that have escaped the confines of the pulsar and form a halo around the PWN.

**Key words.** pulsars: individual: PSR J1813-1749 – gamma rays: general – gamma rays: ISM

## 1. Introduction and multiwavelength context

Results of the Galactic plane survey, conducted by the Imaging Atmospheric Cherenkov Telescope (IACT) Array; High Energy

Stereoscopic System (H.E.S.S.); Aharonian et al. 2005, 2006b; H.E.S.S. Collaboration 2018), have significantly improved our knowledge of the Galactic very-high-energy (VHE) γ-ray sky. This survey detected many previously unknown sources, many of which have multiwavelength counterparts, shedding light on the acceleration mechanism of charged particles. However, due

<footnote>* Corresponding authors; contact.hess@hess-experiment.eu</footnote>
to a large number of possible acceleration sites, a firm identifi-
cation of VHE sources often remains elusive, with high-energy
$\gamma$ rays being mainly produced by either the interaction of leptons
with interstellar radiation fields or inelastic collisions of hadronic
cosmic rays with nuclei of the ISM (ISM) such as molecular
clouds. Therefore, $\gamma$ rays are a valuable probe for understanding
the acceleration of particles to extreme energies.

HESS J1813–178 is a bright source discovered in 2006, located directly on the Galactic plane (GLON, 12.82 ± 0.03°; and GLAT, −0.03 ± 0.02°) with a 1 σ radius of σ = (0.050 ± 0.004)° (H.E.S.S. Collaboration 2018). This discovery was confirmed by the IACT system MAGIC (Albert et al. 2006), and the detection of an X-ray counterpart was reported (Ubertini et al. 2005; Funk et al. 2007; Helfand et al. 2007). HESS J1813–178 is positionally coincident with the young shell-type radio supernova remnant (SNR) G012.8–00.0, with an estimated age of 300–2500 yr (Brogan et al. 2005a; Camilo et al. 2021) and the SNR G012.7–00.0, as well as the young, energetic pulsar PSR J1813–1749 with a characteristic age of 5600 yr and a spin-down luminosity of $\dot{E}$ = 5.6 × 10³⁷ ergs⁻¹ (Gotthelf & Halpern 2009; Camilo et al. 2021).

In addition to the compact emission, extended emission, labelled as HGPSC 063, was observed. However, it was believed to be a potential background artefact and it was therefore discarded (H.E.S.S. Collaboration 2018).

In contrast to the 0.036° extended source detected at TeV, an analysis of GeV Fermi Large Area Telescope (LAT) data revealed emission that is positionally coincident with HESS J1813–178, but with an extension of 0.6° ± 0.06° (Araya 2018). A recent study of the region with Fermi-LAT supports the detection of extended emission and estimates a spatial extent of the emission of 0.56° in an energy range of 1–20 GeV (Xin & Guo 2021). While previous studies of the TeV emission favour a leptonic origin connected to the pulsar (Funk et al. 2007), Araya (2018) concluded that the origin of the GeV $\gamma$-ray emission is more likely to be cosmic rays accelerated at the shock fronts of the SNR or the HII star-forming region W33, located at a distance of 10' from HESS J1813–178. In this case, cosmic rays collide with ambient gas, producing $\pi_0$ particles which then decay into two $\gamma$-ray photons.

Studies of the region have also been conducted by HAWC, revealing the extended emission 3HWC J1813-174, with an extension of less than 0.5° (Albert et al. 2020). LHAASO observations of this region also reveal extended emission. The source LHAASO J1814-1719u* observed with an extension of (0.71 ± 0.07)° by the water Cherenkov detector array (WCDA) and an extension of <0.27° in the one square kilometer array (KM2A) data, as well as the source LHAASO J1814-1636u only observed by the KM2A array with an extension of (0.68 ± 0.08)° (Cao et al. 2024). It is important to note, however, that there are intrinsic uncertainties in the association of sources detected by the WCDA and KM2A LHAASO facilities, with LHAASO J1814-1719u* indicated as an uncertain merger, such that the extended LHAASO J1814-1719u* observed by WCDA may actually correspond to the extended LHAASO J1814-1636u observed by KM2A.

The confirmation of the detection of large, extended emission
in the region around HESS J1813-178 by the water Cherenkov
detectors, and improvements in background rejection and event
reconstruction have recently led to the discovery of previously
undetected, large extended emission in IACT data (Abdalla et al.
2021; H. E. S. S. Collaboration 2023). This is a particularly inter-
esting development for the region around HESS J1813-178 and

warrants a reanalysis of the region to resolve the disagreement
between the different observed morphologies.

While pulsar halos have been detected in evolved systems, regions powered by the emission from young pulsars usually present a confined pulsar wind nebula (PWN; Giacinti et al. 2020). Therefore, a significant detection of extended $\gamma$-ray emission in the TeV energy range around PSR J1813-1749 will help to resolve the discrepancy in extension between GeV and TeV energies. Observing the escape of electrons from the PWN into the ISM at such a young age is highly unlikely therefore, the detection of extended emission resulting from highly relativistic electrons could hint at an untypical behaviour of the pulsar or an early interaction of the reverse shock of the SNR with the pulsar wind, disrupting the system and facilitating an early release of particles into the surrounding ISM.

We present an updated analysis of the H.E.S.S. data in the region around PSR J1813-1749. To obtain a consistent description of the region, we carried out an analysis of the Fermi-LAT data, with increased exposure compared to the study conducted in Araya (2018), and a joint likelihood minimisation of these datasets, allowing us to analyse the data from both experiments simultaneously and obtain a consistent description. We additionally present a leptonic and hadronic model in order to explain the emission observed around PSR J1813-1749.

## 2. Data analysis

2.1. H.E.S.S.

2.1.1. Data selection

The H.E.S.S. array consists of four imaging atmospheric Cherenkov telescopes with a mirror diameter of 12 m, constructed in 2003, located in the Khomas Highland in Namibia (Aharonian et al. 2006a). A fifth 28 m diameter telescope was added to the array in 2012. The telescopes are capable of detecting photons in an energy range from a few tens of GeV up to 100 TeV. The data used in this analysis was obtained between 2004 and 2010 as part of different observation campaigns in the sky region, therefore the fifth telescope has not been used in the following analysis.

H.E.S.S. observations consist of data collection intervals of ~28 min, called observation runs. The selection of suitable data is based on the observation characteristics for each run. For this data analysis, participation of at least three telescopes per run was required. An additional requirement was the pointing position of the telescopes, which was required to not be offset by more than 2.0° from RA = 273.40°, Dec = -17.84°, the position of HESS J1813-178 as derived in Aharonian et al. (2006b).

After applying standard selection cuts used in spectral analysis (Aharonian et al. 2006a), this selection resulted in a dataset with a total exposure of 32 h. The average angular distance of the position of HESS J1813-178 to the centre of the field of view (FoV) in all observations of the dataset is 1.28°, as most observations were taken around a neighbouring source, HESS J1809-193 (Fig. C.2, H.E.S.S. Collaboration 2023).

The data reduction was performed using the H.E.S.S. analysis package (see Aharonian et al. 2006a). The data was reconstructed using the ImPACT (Image pixel-wise fit for atmospheric Cherenkov Telescope) algorithm, which determines the best-fit parameters of the $\gamma$-ray-like events by using a template-based maximum-likelihood approach (Parsons & Hinton 2014).

An important source of background in the analysis of IACT
data are atmospheric air showers initiated by charged cosmic-ray
particles. Cosmic-ray showers outnumber the showers initiated by γ rays by several orders of magnitude. To describe this background across the whole field of view, a template model was used. This model is constructed from a large set of runs taken in regions without known γ-ray sources, following the scheme described in Mohrmann et al. (2019). The flux normalisation and spectral index of this field of view background template are used to account for the different observation conditions in each run (for more information see Appendix A).

For each individual run, a safe energy threshold was defined to exclude the threshold region where effective detection areas vary steeply with energy, and where the energy reconstruction is biased. The threshold was chosen as the energy at which the energy bias (the deviation between true and reconstructed event energy; Aharonian et al. 2006a) is below 10%. A second energy threshold was defined as the peak in the field of view background template spectrum (Mohrmann et al. 2019); the higher of the two thresholds was applied for the run.

The events passing the cut criteria were then stacked into a single dataset and the field of view background template, as discussed in Mohrmann et al. (2019), was applied. In addition to the run-wise threshold, a global energy threshold for the dataset of 0.4 TeV was set. All events below this energy were excluded from the binned-likelihood estimation.

## 2.1.2. Data processing

The open-source Python package gammapy, version 0.18.2 (Deil et al. 2020; Donath et al. 2023) was used for the data analysis, to characterise the spectral and morphological properties of the γ-ray emission. Gammapy enables a three-dimensional likelihood analysis that stores the data and models in data cubes, allowing a simultaneous parameter optimisation of the model in the two spatial dimensions and energy. All results were cross-checked using an independent calibration and reconstruction scheme (Khelifi et al. 2015).

The background model flux normalisation and spectral index were fitted to each run used in the analysis to account for differences in atmospheric conditions and changes in the optical efficiency of the telescopes over time. The list of events passing the selection criteria was then combined into a data cube of 4° × 4°, centred on the previously derived position of HESS J1813−178 in the equatorial coordinate system. This size ensures a sufficiently large off-source region for estimating the rate of remaining hadronic background events, without adding a large number of unrelated γ-ray sources. Spatial bins of 0.02° size, and a logarithmic energy binning of 25 bins between 0.1 TeV, and 100 TeV were chosen for this study.

A likelihood fit to determine the best description of the emission in the FoV was then conducted. For this purpose symmetric and elongated 2-dimensional Gaussian source models were used to describe the morphology. The symmetric Gaussian is defined as:

$$ \Phi(\Theta) = \frac{1}{2\pi\sigma_M^2} \exp\left(-\frac{1}{2}\frac{\Theta^2}{\sigma_M^2}\right), \qquad (1) $$

with $\sigma_M$ being the major axis of the Gaussian and $\Theta$ the angular distance to the model centre. In the case of an elongated Gaussian, the major axis $\sigma_M$ is defined by $\sigma_{eff}$, the effective containment radius:

$$ \sigma_{eff}(\Delta\phi) = \sqrt{(\sigma_M \sin(\Delta\phi))^2 + (\sigma_M \cdot \sqrt{1 - e^2}) \cos(\Delta\phi)^2}, \quad (2) $$

with $e$ being the eccentricity and $\Delta\phi$ the difference between the position angle of the Gaussian $\phi$ and the position angle of the evaluation point.

As a spectral model, a power law was used. This model is described by:

$$ \frac{dN}{dE} = N_0 \cdot \left(\frac{E}{E_0}\right)^{-\Gamma}, \qquad (3) $$

where $N_0$ is the flux normalisation factor, $\Gamma$ the spectral index and $E_0$ the reference energy, which is set to 1 TeV for the whole study. Alternatively, a logarithmic parabola spectral model was used:

$$ \frac{dN}{dE} = N_0 \cdot \left(\frac{E}{E_0}\right)^{-\Gamma - \beta \log\left(\frac{E}{E_0}\right)}, \qquad (4) $$

where $\beta$ is the curvature, as well as a power law with a generalised exponential cut-off

$$ \frac{dN}{dE} = N_0 \cdot \left(\frac{E}{E_0}\right)^{-\Gamma} \cdot \exp(-\lambda E), \qquad (5) $$

with $\lambda$, the inverse of the cut-off energy.

Significance maps of the region were computed using the method described in Li & Ma (1983). For all sky maps, a correlation radius of 0.4° was applied, since this has proven beneficial for detecting extended emission in previous studies (Abdalla et al. 2021). The choice of correlation radius has no influence on the analysis, it only improves the visibility of large extended emission.

## 2.1.3. Estimation of systematic uncertainties

The Instrument Response Functions (IRFs) describe the response of the system under various observation conditions. The IRFs for each run were obtained by interpolating between IRFs generated from Monte Carlo simulations covering a grid of observational conditions (Bernlöhr 2008). This method cannot describe all variations of observational conditions in detail, which leads to systematic uncertainties affecting the likelihood minimisation. In this analysis, three possible sources of uncertainties were considered in our results.

Firstly, a discrepancy between the true and reconstructed pointing positions will lead to uncertainties in source position estimation. As indicated in Acero et al. (2010), the uncertainty of the pointing is 20", estimated using the position of known γ-ray sources and stars.

Secondly, the reconstruction of the energy of a γ ray strongly depends on the optical efficiency of the telescopes, influenced by effects such as degradation of the photomultiplier tubes or mirrors over time. Additionally, the transparency of the atmosphere, affected by effects such as the presence of aerosols, strongly influences the accuracy of the energy reconstruction. A discrepancy in energy reconstruction particularly influences the estimated spectral parameters of the source.

Thirdly, an important aspect of the data analysis is the correct statistical description of the inferred background event rate as a function of energy. The background model used in this analysis was constructed using observation runs taken under not fully identical observation conditions from those runs used in this analysis. The statistical and systematic uncertainties in the inferred background model will lead to uncertainty in the prediction of background counts and affect all fit parameters.
<table><caption>**Table 1.** Best-fit parameters obtained for the likelihood minimisation of H.E.S.S. data, following the spatial and spectral description introduced in Sect. 2.1.2.</caption><thead><tr><th></th><th>HESS J1813-178A</th><th>HESS J1813-178B</th></tr></thead><tbody><tr><td>Γ</td><td>2.17 ± 0.05<sub>stat</sub> ± 0.03<sub>sys</sub></td><td>2.36 ± 0.09<sub>stat</sub> ± 0.05<sub>sys</sub></td></tr><tr><td>N<sub>0</sub> at 1 TeV [10<sup>-12</sup>cm<sup>-2</sup> s<sup>-1</sup> TeV<sup>-1</sup>]</td><td>3.16 ± 0.18<sub>stat</sub> ± 0.24<sub>sys</sub></td><td>9.89 ± 1.40<sub>stat</sub> ± 1.23<sub>sys</sub></td></tr><tr><td>RA [°]</td><td>273.400 ± 0.004<sub>stat</sub> ± 0.001<sub>sys</sub></td><td>273.61 ± 0.06<sub>stat</sub> ± 0.07<sub>sys</sub></td></tr><tr><td>Dec [°]</td><td>-17.831 ± 0.004<sub>stat</sub> ± 0.001<sub>sys</sub></td><td>-17.39 ± 0.07<sub>stat</sub> ± 0.08<sub>sys</sub></td></tr><tr><td>σ<sub>M</sub> [°]</td><td>0.056 ± 0.003<sub>stat</sub> ± 0.001<sub>sys</sub></td><td>0.72 ± 0.08<sub>stat</sub> ± 0.09<sub>sys</sub></td></tr><tr><td>e</td><td>—</td><td>0.80 ± 0.06<sub>stat</sub> ± 0.04<sub>sys</sub></td></tr><tr><td>Δφ [°]</td><td>—</td><td>40 ± 7<sub>stat</sub> ± 6<sub>sys</sub></td></tr></tbody></table>

To assess the systematic errors, a large number of pseudo-datasets with a shift of the energy scale, as well as a shift of the amplitude and index of the background model were created, and a linear gradient to the hadronic background model was applied. Then, a likelihood minimisation for these pseudo-datasets was performed. Using this method, the spread of the resulting fit values can be interpreted as the systematic error of the respective parameter. In addition, the systematic error in the pointing position is added. A more detailed description of the method can be found in Appendix B, and the systematic errors estimated for the source parameters are presented in Table 1.

## 2.2. Fermi data

The LAT on board the Fermi satellite detects γ rays between 20 MeV and more than 300 GeV from the entire sky by pair-conversion within the detector (Atwood et al. 2009). In this study, the data from the beginning of the mission in August 2008 until October 2021 was analysed. The most recent IRFs from Pass 8 version 3 – P8R3_SOURCE_V2 (Ajello et al. 2021) were used and the catalogue models were taken from the 12-yr 4FGL source catalogue (Abdollahi et al. 2020). A 6° region of interest (ROI) around the source position of 4FGL J1813–1737e, as derived in Araya (2018), a bin size of 0.025° and 8 energy bins per decade with logarithmic spacing were used. Only events arriving with a maximum zenith angle of 90° were considered, in order to avoid including secondary γ rays from the Earth’s horizon.

The background was modelled using a isotropic background model (iso_P8R3_SOURCE_V3_v1.txt), and a Galactic diffuse background model (gl1_iem_v07.fits). In order to describe the region well, all known catalogue sources in a region of 10° around the source position were included.

For the analysis, the Python package fermipy, version 1.0.1 (Wood et al. 2017), optimised for performing a binned likelihood analysis of Fermi-LAT data using the Fermi Science Tools, was used. With fermipy a counts cube and IRF cubes in an energy range between 1 GeV and 1 TeV was created. For the analysis of the source, the data, IRFs, and background models were exported to gammapy. In contrast to the standard analysis chain of the Fermi-LAT data, this allows us to fit both morphological and spectral parameters simultaneously.

# 3. Results

## 3.1. H.E.S.S. results

A significance map of the region, computed with a correlation radius of 0.4°, is shown in Fig. 1. A significance map computed with a correlation radius of 0.06°, which is equivalent to

the point-spread-function (PSF) of the analysis, has also been computed and can be seen in Fig. C.1. Bright γ-ray emission around the pulsar PSR J1813–1749 is visible. Additionally, one can observe diffuse, extended emission around the pulsar, as well as γ-ray emission from the source HESS J1809–193 an extended source potentially powered by the pulsar PSR J1809–1917. The γ-ray emission from this source was accounted for by using the best-fit description derived in H.E.S.S. Collaboration (2023). The $1\sigma$ extends of the two Gaussian models used in the analysis are indicated by the grey dashed lines in the right panel of Fig. 1.

The compact emission can be accounted for using a Gaussian model, centred at a position of RA = (273.396 ± 0.004)°, Dec = (-17.831 ± 0.004)° and σ = (0.056 ± 0.003)°, hereafter referred to as HESS J1813–178A. The emission is well described with a power-law spectral model and detected with a significance of $38\sigma$. The best-fit position and extension derived in this analysis are compatible with the results reported in Aharonian et al. (2006b).

In addition to the compact emission centred at the location of PSR J1813–1749, a second extended emission component was detected during the Galactic plane survey (H.E.S.S. Collaboration 2018), but discarded because of large uncertainties in the background estimation. The emission was reported to be located at RA = (273.46 ± 0.05)°, Dec = (-17.77 ± 0.06)° with a Gaussian width of σ = (0.31 ± 0.06)° (H.E.S.S. Collaboration 2018). In addition to the improved image reconstruction technique ImPACT (Parsons & Hinton 2014), this study used a background template created from a large number of archival observations in order to estimate the background events, as opposed to the former approach, which only used a ring-shaped area around the source from the observation run itself and is therefore heavily biased by the initial assumption of an extension of the emission (Aharonian et al. 2006b). With these improvements, a firm detection of significant extended emission in the source region is possible. The morphological model which best describes the emission is evaluated, testing both a disc-like model and a Gaussian model, each with and without the possibility of an elongation. This study finds that an elongated Gaussian model centred at RA = (273.61 ± 0.06)°, Dec = (-17.39 ± 0.07)° with a semi-major axis of σ = (0.72 ± 0.08)° describes the emission best. The gaussian model is preferred over the disc-like model with $\Delta T_S = 59.91$, while the elongated Gaussian model is preferred by $3.4\sigma$ compared to the symmetric Gaussian model. The increased size compared to the first report of possible extended emission can be attributed to the difference in background estimation. A power-law model was used to describe the emission. This additional extended emission, hereafter referred to as HESS J1813–178B, was detected with a significance of $13\sigma$. The best-fit spectral and spatial
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Significance maps of the data taken by H.E.S.S. in the energy range between 0.4 and 100 TeV. For both maps, a correlation radius of 0.4° was used. Left: significance map of the region around PSR J1813−1749. The contours correspond to the 3 and 5σ regions. The dashed line indicates the position of the Galactic plane. Right: significance map of the region after subtracting the emission using four components. The blue dots indicate the positions of both pulsars. The 1σ Gaussian extend of the models used to describe the emission are indicated by the grey and black lines.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Spectra and spectral energy distribution (SED) of the H.E.S.S. γ-ray emission from HESS J1813−178A and the extended emission from HESS J1813−178B. The SED derived in Aharonian et al. (2006b) are shown as filled triangles.</figcaption></figure>

parameters of both models are shown in Table 1. The best-fit position and extension along the major axis of components A and B are shown in Fig. 1 by the black lines. The morphology of HESS J1813−178B is in good spatial agreement with the extended emission 3HWC J1813−174 observed by HAWC, as well as 1LHAASO J1814−1719u*, extended emission observed by the WCDA of LHAASO (see Fig. E.3). The spectra of both components are depicted in Fig. 2. Additionally, the morphological structure is tested for energy-dependence, the results of this study are presented in Sect. 3.3.

## 3.2. Fermi-LAT results

A study using nine years of data revealed extended emission, later referred to as 4FGL J1813−1737e, in the region around PSR J1813−1749 (Araya 2018). They reported significant, symmetric, γ-ray emission with an extension of (0.60 ± 0.06)° and a power-law spectral index of 2.14 ± 0.04. This study reanalyses the region using an increased exposure of 12 yr.

In Fig. 3, a significance map of the emission in the region around 4FGL J1813−1737e is shown. After the verification of the data export from fermipy to gammapy (Appendix D), the morphological and spectral properties of the emission are estimated. The emission in the ROI is best described by an elongated Gaussian model centred at RA = (273.33 ± 0.05)° and Dec = (−17.57 ± 0.05)°, with σ = (0.56 ± 0.07)°. This position and size are similar to the source position of HESS J1813−178B observed with H.E.S.S. (see Fig. E.1), indicating that the emission from 4FGL J1813−1737e and HESS J1813−178B could be connected. The best-fit parameters are shown in Table 2. The significance of the fitted model is 35σ. After accounting for the emission from this component, residual emission spatially overlapping with the pulsar is observed. Therefore a second component with a symmetric Gaussian model was added to the data, to account for this residual emission. Because of the spectral shape below 10 GeV, a simple power-law spectral model cannot reflect the observed emission well. Therefore a power law with exponential cut-off was chosen as a spectral model. Minimisation yields a compact component, positionally coincident with PSR J1813−1749. The residual significance map after the application of both models is shown on the right-hand side of Fig. 3. Because the shape of the SED does not suggest a possible connection to HESS J1813−178A, this compact component will be referred to as component C. The best-fit parameters for 4FGL J1813−1737e and component C are given in Table 2. The cut-off of this model was estimated to be 1.81 GeV. Component C improved the likelihood of the model with a significance of 3.8σ. To detect this second component with certainty and probe the nature of this emission, events with energies below 1 GeV would need to be included. The SEDs are given in Fig. 4. Differences between the SED derived in this analysis and the analysis performed in Araya (2018) can be attributed to the changes from the 3FGL to the 4FGL catalogue, the usage of newly computed IRFs from Pass 8 version 3, and four years of additional data used in this analysis.

In order to estimate the systematic uncertainties on the best-fit flux normalisation factor, two separate factors are considered. In order to estimate the error introduced by the instrument itself,
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Significance maps of the Fermi-LAT data around the position of 4FGL J1813-1737e in the energy range from 1 GeV to 1 TeV, with a correlation radius of 0.1°. Left: the position of the Galactic plane is indicated by the dashed line, the pulsar position is indicated in blue, and 3σ and 5σ contours are depicted. Right: significance map of the region after subtracting the emission using an elongated Gaussian model and a symmetric Gaussian model. The 1σ Gaussian extent of the models is depicted by the black dashed lines.</figcaption></figure>

<table><caption>**Table 2.** Best-fit source parameters obtained for the Fermi-LAT data for a model consisting of an elongated Gaussian model with a power law and a symmetric Gaussian model with an exponential cut-off power law as a spectral model.</caption><thead><tr><th></th><th>Component C</th><th>4FGL J1813-1737e</th></tr></thead><tbody><tr><td>Γ</td><td>0.95 ± 0.22<sub>stat</sub></td><td>2.17 ± 0.09<sub>stat</sub></td></tr><tr><td>λ [TeV<sup>-1</sup>]</td><td>551.5 ± 84.4<sub>stat</sub></td><td>—</td></tr><tr><td>N<sub>0</sub> at 2 GeV [10<sup>-12</sup> cm<sup>-2</sup> s<sup>-1</sup> TeV<sup>-1</sup>]</td><td>8.4 ± 2.5<sub>stat</sub> ± 1.83<sub>sys</sub></td><td>2.7 ± 0.6<sub>stat</sub> ± 0.40<sub>sys</sub></td></tr><tr><td>RA [°]</td><td>273.452 ± 0.004<sub>stat</sub></td><td>273.33 ± 0.05<sub>stat</sub></td></tr><tr><td>Dec [°]</td><td>-17.73 ± 0.03<sub>stat</sub></td><td>-17.57 ± 0.05<sub>stat</sub></td></tr><tr><td>σ<sub>M</sub> [°]</td><td>0.21 ± 0.03<sub>stat</sub></td><td>0.56 ± 0.07<sub>stat</sub></td></tr><tr><td>e</td><td>—</td><td>0.78 ± 0.07<sub>stat</sub></td></tr><tr><td>Δφ [°]</td><td>—</td><td>33 ± 8<sub>stat</sub></td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** SED and best-fit spectral models of the two-component description derived from the analysis of the data obtained by the Fermi-LAT satellite, and SED derived in Araya (2018).</figcaption></figure>

the effective area has been scaled up and down by 3%, and the analysis has been repeated. This procedure is explained in detail in Fermi-LAT Collaboration (2019). The errors estimated from this method are 0.36 × 10<sup>-12</sup> cm<sup>-2</sup> s<sup>-1</sup> TeV<sup>-1</sup> for

4FGL J1813-1737e and 1.79 × 10<sup>-12</sup> cm<sup>-2</sup> s<sup>-1</sup> TeV<sup>-1</sup> for component C.

In order to estimate the error introduced by the background estimation, the residual flux inside a region corresponding to the source region, at an empty position in the ROI, has been calculated. The systematic uncertainty on the number of measured counts can then be derived by dividing the total amount of residual counts by the total amount of background in the shifted region and then multiplying by the total amount of background counts in the original source region. The number of counts per bin is then scaled up and down for the whole ROI and the likelihood minimisation is repeated. This process has been repeated for ten different regions. This estimation yields an systematic uncertainty of 0.18 × 10<sup>-12</sup> cm<sup>-2</sup> s<sup>-1</sup> TeV<sup>-1</sup> for 4FGL J1813-1737e and 0.36 × 10<sup>-12</sup> cm<sup>-2</sup> s<sup>-1</sup> TeV<sup>-1</sup> for component C. Both errors were then added in quadrature.

### 3.3. Energy-dependent morphology

The spatial association of the γ-ray emission with PSR J1813-1749, and an X-ray PWN detected by the INTEGRAL satellite, XMM-Newton and Chandra (Ubertini et al. 2005; Funk et al. 2007; Helfand et al. 2007; Fig. E.2), suggest that a leptonic origin of the emission is plausible. Previous analyses of evolved PWNe suggest an energy-dependent
<table><caption>**Table 3.** Best-fit parameters for the morphology of HESS J1813–178B in different energy bands, for both *Fermi*-LAT (F1 – F6) and H.E.S.S. datasets (H1 large – H3 large).</caption><thead><tr><th></th><th>Energy [GeV]</th><th>RA [°]</th><th>Dec [°]</th><th>σ<sub>M</sub> [°]</th><th>e</th></tr></thead><tbody><tr><td>F1</td><td>1.0–2.0</td><td>273.46 ± 0.02</td><td>-17.67 ± 0.02</td><td>0.39 ± 0.04</td><td>0.81 ± 0.05</td></tr><tr><td>F2</td><td>2.0–4.0</td><td>273.48 ± 0.03</td><td>-17.69 ± 0.03</td><td>0.42 ± 0.04</td><td>0.85 ± 0.04</td></tr><tr><td>F3</td><td>4.0–7.5</td><td>273.38 ± 0.04</td><td>-17.66 ± 0.06</td><td>0.39 ± 0.08</td><td>0.87 ± 0.07</td></tr><tr><td>F4</td><td>7.5–18</td><td>273.26 ± 0.09</td><td>-17.52 ± 0.12</td><td>0.53 ± 0.11</td><td>0.75 ± 0.19</td></tr><tr><td>F5</td><td>18–58</td><td>273.17 ± 0.10</td><td>-17.77 ± 0.08</td><td>0.57 ± 0.16</td><td>0.77 ± 0.21</td></tr><tr><td>F6</td><td>58–1.0 × 10<sup>3</sup></td><td>273.22 ± 0.20</td><td>-17.75 ± 0.13</td><td>0.71 ± 0.19</td><td>0.92 ± 0.08</td></tr><tr><td>H1 small</td><td>(0.4–1.3) × 10<sup>3</sup></td><td>273.393 ± 0.005</td><td>-17.832 ± 0.005</td><td>0.054 ± 0.004</td><td>—</td></tr><tr><td>H2 small</td><td>(1.3–5.7) × 10<sup>3</sup></td><td>273.397 ± 0.006</td><td>-17.834 ± 0.007</td><td>0.064 ± 0.005</td><td>—</td></tr><tr><td>H3 small</td><td>(5.7–100) × 10<sup>3</sup></td><td>273.408 ± 0.009</td><td>-17.817 ± 0.009</td><td>0.035 ± 0.006</td><td>—</td></tr><tr><td>H1 large</td><td>(0.4–1.3) × 10<sup>3</sup></td><td>273.61 ± 0.08</td><td>-17.48 ± 0.09</td><td>0.87 ± 0.12</td><td>0.88 ± 0.04</td></tr><tr><td>H2 large</td><td>(1.3–5.7) × 10<sup>3</sup></td><td>273.41 ± 0.09</td><td>-17.42 ± 0.11</td><td>0.48 ± 0.09</td><td>0.48 ± 0.38</td></tr><tr><td>H3 large</td><td>(5.7–100) × 10<sup>3</sup></td><td>273.65 ± 0.11</td><td>-17.18 ± 0.10</td><td>0.41 ± 0.09</td><td>0.00 ± 0.02</td></tr></tbody></table>

**Notes.** An elongated Gaussian model was assumed in all cases and the $1\sigma$ statistical uncertainties are given. Additionally, the best-fit values for the symmetrical Gaussian were used to account for the emission from HESS J1813–178A.

morphology caused by electron diffusion and cooling (H.E.S.S. Collaboration 2019; MAGIC Collaboration 2020; Principe et al. 2020). To test for energy dependence of the extended emission, the data observed with H.E.S.S. was divided into three energy ranges (H1 – H3, specified in Table 3) and the source models for HESS J1813–178A and HESS J1813–178B were fitted in each range. The *Fermi*-LAT data was divided into six energy ranges (F1 - F6, also specified in Table 3) and the source model for 4FGL J1813–1737e was refitted in each range. The spacing of these energy ranges was chosen to account for the small amount of statistics in the H.E.S.S. data and the high-energy range in the *Fermi*-LAT data. While F1 - F6 span 3 energy bins each, F6 spans 11 energy bins. The H.E.S.S. data was divided such that H1 and H2 contain five energy bins of the dataset, while H3 spans 10 energy bins.

Table 3 shows the extension along the semi-major axis, as well as the eccentricity of the ellipse describing the emission in each energy bin and the distance of the centre of the spatial model to the pulsar. In Fig. 5, the dependence of the $1\sigma$ containment area of the ellipse and the offset to the centre of the fitted model from the pulsar position as a function of the energy is depicted. The area of the ellipse is defined as $A = \pi\sigma_M^2\sqrt{1 - e^2}$. The binning was chosen based on the available statistics for each dataset.

The enclosed area is compatible within errors in the overlapping energy ranges of the *Fermi*-LAT and H.E.S.S. data. There is no significant indication for energy-dependence of the best-fit containment area of HESS J1813–178A and HESS J1813–178B. This study does however reveal a dependence of the distance between the pulsar and the centre of the emission, which increases with increasing energy. As the morphology of HESS J1813–178A remains fixed throughout all energy bands, the increased offset of the best-fit position towards higher energies may indicate that the particle transport occurs preferentially towards a single direction before spreading out more isotropically as the particles lose energy and cool. Alternatively, this could indicate that a second faint source is present in the region.

### 3.4. Joint model fit

A consistent description of the ROI through the GeV-TeV range can be reached by fitting the *Fermi*-LAT and H.E.S.S. data

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Best-fit angular offset and area of the ellipse of HESS J1813–178A and HESS J1813–178B with $1\sigma$ uncertainties, as well as 4FGL J1813–1737e data in independent energy bins (see Table 3) H1 – H3 and F1 – F6.</figcaption></figure>

jointly. Following the results derived in the analysis of the respective datasets, three source models were used to describe the data across five decades of energy. First a symmetric Gaussian model, as well as an elongated Gaussian model was added to both datasets. These models will be referred to as component A and B respectively. Additionally, a symmetric Gaussian model with an exponential cut-off power-law spectral model was added only to the *Fermi*-LAT data. This model component, referred to as component C in the analysis of the *Fermi*-LAT data. The likelihood minimisation was then performed on both datasets at the same time. The best-fit parameters are shown in Table 4.

Figure 6 depicts the spectra and SED of the best-fit models, together with the sensitivity of the LAT for 12 yr of exposure in blue. The curve shows the broadband sensitivity for sources located in the Galactic plane (Ajello et al. 2021). While the best-fit parameters of the models where estimated using the data from both datasets at the same time, the flux points were computed in the respective datasets. This results in a energy range where both the H.E.S.S. and *Fermi*-LAT data show significant flux or upper limits. If the spectrum of component
<table><caption>**Table 4.** Best-fit parameters derived in the joint analysis of the H.E.S.S. and Fermi-LAT data.</caption><thead><tr><th></th><th>Comp. A</th><th>Comp. B</th><th>Comp. C</th></tr></thead><tbody><tr><td>Γ</td><td>2.05 ± 0.03</td><td>2.17 ± 0.03</td><td>-0.12 ± 0.02</td></tr><tr><td>β</td><td>0.0620 ± 0.003</td><td>0.0437 ± 0.008</td><td>—</td></tr><tr><td>α</td><td>—</td><td>—</td><td>0.66 ± 0.01</td></tr><tr><td>λ [10<sup>3</sup> · TeV<sup>-1</sup>]</td><td>—</td><td>—</td><td>3.93 ± 0.20</td></tr><tr><td>N<sub>0</sub> at 1 TeV [10<sup>-12</sup> cm<sup>-2</sup> s<sup>-1</sup> TeV<sup>-1</sup>]</td><td>3.16 ± 0.14</td><td>6.47 ± 0.43</td><td>9.89 ± 1.64</td></tr><tr><td>RA [°]</td><td>273.400 ± 0.003</td><td>273.39 ± 0.03</td><td>273.48 ± 0.02</td></tr><tr><td>Dec [°]</td><td>-17.832 ± 0.003</td><td>-17.50 ± 0.04</td><td>-17.67 ± 0.02</td></tr><tr><td>σ<sub>M</sub> [°]</td><td>0.056 ± 0.003</td><td>0.54 ± 0.03</td><td>0.29 ± 0.02</td></tr><tr><td>e</td><td>—</td><td>0.73 ± 0.04</td><td>—</td></tr><tr><td>Δφ [°]</td><td>—</td><td>33 ± 7</td><td>—</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** SED from a joint-analysis of the combined H.E.S.S. and Fermi-LAT data. The 12 yr sensitivity of the LAT is indicated by the blue line.</figcaption></figure>

A does not change drastically at the lower energies, measuring it in the LAT data would not be possible as components A and C are positionally coincident and the detector sensitivity is insufficient. Together these three different source models describe the emission of the ROI well. The advantage of this joint model becomes evident in the GeV energy regime. The spatial separation between 4FGL J1813–1737e and the compact emission component detected in the Fermi-LAT data is complicated, but using the assumption that 4FGL J1813–1737e and HESS J1813–178B have the same origin, the influence of the second component on the model parameters can be reduced. This becomes evident when comparing the flux around 1 GeV, which has been reduced by a factor two in the joint model.

## 4. Spectral modelling using GAMERA

Two scenarios were investigated for the origin of the extended γ-ray emission in the region of PSR J1813–1749. In the first scenario, a magnetised wind of charged particles, primarily electrons and positrons, originating from PSR J1813–1749, forms a PWN. This leptonic emission scenario has been favoured in previous analyses of the X-ray and TeV emission (Ubertini et al. 2005; Funk et al. 2007). In systems older than 10 kyr, the electrons can escape the confines of the PWNe and diffuse through the ISM, although the age at which this occurs may vary for individual systems. The diffuse γ-ray emission produced by these electrons is typically observed around middle-aged pulsars and

one can distinguish between a majority of the particles still confined in the PWN and one where most particles have escaped into the ISM (Giacinti et al. 2020). This evolutionary stage of a leptonic system could potentially correspond to the extended emission observed as HESS J1813–178B.

Another possible origin of the extended emission are cosmic-ray nuclei, accelerated at shock fronts and which have since escaped from their accelerator. A possible candidate for such an acceleration site is the SNR G12.82–0.02, which is positionally coincident with the pulsar PSR J1813–1749 and believed to be the pulsar’s host SNR (see Fig. E.2). After acceleration and propagation, the nuclei may interact within dense molecular clouds in the region, producing high-energy γ rays.

This emission scenario has been favoured in an analysis of nine years of Fermi-LAT data, performed by Araya (2018). He suggested the stellar cluster Cl J1813-178 as a possible origin, but also states that the extension of the γ-ray emission exceeds the extension of the star-forming region, and therefore an unknown component of star formation would be necessary.

### 4.1. Leptonic scenario

To describe the γ-ray emission in the region using a time-dependent model of electrons originating from PSR J1813–1749, the GAMERA package (Hahn 2015) was used. Additionally to the SED derived in the joint analysis, measurements of the region in the 20–100 keV band, showing an unresolved hard X-ray source observed by the INTEGRAL satellite (Ubertini et al. 2005), XMM-Newton (Funk et al. 2007) (with an extraction region of 75"), and Chandra (Helfand et al. 2007, with an elliptical extraction region of 6" × 8") were used, as well as observation data from the Very Large Array in the 20 and 90 cm band, revealing a shell-like non-thermal radio source identified as SNR G12.82–0.02 (Brogan et al. 2005a). The radio data is only shown for comparison and was not included in the modelling.

The SED used for this model were not derived using the same spatial extent. While one possible solution for this discrepancy would be to scale the observed flux for the different observation regions, this would require the assumption that a linear scaling can be applied. Additionally, X-ray emission generally traces higher energy electrons, and hence younger electrons, than the TeV emission. To account for these differences, three electron ‘generations’ were defined, following the procedure in H.E.S.S. Collaboration (2023). The first electron generation, the relic electrons, were injected since the birth of the pulsar and are now
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** SED derived in this analysis, as well as X-ray SED by XMM-Newton (Funk et al. 2007), Chandra (Helfand et al. 2007), and INTEGRAL (Ubertini et al. 2005), and the SED derived in the analysis of VLA data (Brogan et al. 2005a), compared to the SED curves expected from the leptonic model obtained with GAMERA. Left: full energy range. Right: comparison between the γ ray flux estimated from the leptonic model from 0–1.6 kyrs and the last 2.5 kyrs and the emission observed from HESS J1813–178B and HESS J1813–178A respectively. The shaded grey error bands indicate the possible parameter space, given in Table 6.</figcaption></figure>

<table><caption>**Table 5.** Assumed properties of PSR J1813–1749 following newest estimates from Camilo et al. (2021).</caption><thead><tr><td>Parameter</td><td>Assumed value</td></tr></thead><tbody><tr><td>Distance, d</td><td>6.2 kpc</td></tr><tr><td>Spin-down power, <i>E</i>&#775;</td><td>5.6 × 10<sup>37</sup> erg s<sup>-1</sup></td></tr><tr><td>Spin period, <i>P</i></td><td>44.7 × 10<sup>-3</sup> s</td></tr><tr><td>Change in spin period, <i>P&#775;</i></td><td>1.27 × 10<sup>-13</sup> s s<sup>-1</sup></td></tr><tr><td>Characteristic age <i>τ</i><sub>c</sub></td><td>5.585 kyr</td></tr><tr><td>Braking index, <i>n</i></td><td>3.0</td></tr><tr><td>Break energy, <i>E</i><sub>b</sub></td><td>100 GeV</td></tr></tbody></table>

**Notes.** For the braking index and break energy canonical values are assumed.

detected as extended high-energy γ-ray emission. The second electron generation, the middle-aged electrons, were injected since $t_e = t_{\text{true}} \cdot t_{\text{PWN}}$, with $t_{\text{true}}$ the age of the system, and are now detected as PWN around PSR J1813–1749. The third electron generation, the young electrons injected since $t_e = t_{\text{true}} \cdot t_{\text{X-ray}}$ are responsible for the X-ray synchrotron emission. For this model, an electron injection spectrum following a power law with a spectral index of $\alpha$, a break energy of $E_b$ and a cutoff $E_{\text{cut}}$:

$$N(E_e) = \left(1 + \frac{E_e}{E_b}\right)^\alpha \cdot \exp\left(\frac{-E_e}{E_{\text{cut}}}\right), \qquad (6)$$

where $E_e$ the electron energy, was used.

The distance to the pulsar, as well as pulsar-specific parameters, such as spin-down period $P$, change in spin period $\dot{P}$, and spin-down power $\dot{E}$ were taken from Camilo et al. (2021). All assumed properties of the pulsar used for this model are given in Table 5.

The injection index $\alpha$, the birth period $P_0$, and the conversion fraction of the spin-down luminosity to electrons $\vartheta$, as well as the magnetic field present at the time of the observation of the γ-ray emission $B_{\text{now}}$, and the age of the electron populations were free parameters of the model. The spin-down power of the pulsar was assumed to evolve in time with:

$$\dot{E} = \dot{E}_0 \left(1 + \frac{t}{\tau_0}\right)^{-2}, \qquad (7)$$

with $\tau_0 = P_0/(2 \dot{P}_0)$, the spin-down time at the birth of the pulsar and $\dot{E}_0$ the initial spin-down energy loss rate at birth of the pulsar. Similarly the pulsar period was assumed to evolve as $P(t) = P_0 \left(1 + \frac{t}{\tau_0}\right)^{0.5}$, with $P_0$ the birth period of the pulsar and the magnetic field as $B(t) = B_0 \left(1 + \left(\frac{t}{\tau_0}\right)^{0.5}\right)^{-1}$ (Gaensler & Slane 2006; Venter & de Jager 2007) with $B_0$ the magnetic field at birth of the pulsar. The spin down luminosity is defined as:

$$L = L_0 \cdot \vartheta \cdot \left(1 + \frac{t}{\tau_0}\right)^{-2} \qquad (8)$$

with the initial spin down luminosity $L_0 = \dot{E} \cdot (1 + (t_{\text{true}}/t_0))^2$ and $t_0 = \frac{P_0^2}{2P\dot{P}}$.

Since the time at which the escape of the particles into the ISM first occurs is unknown, this modeling approach assumes that the particles have been able to escape instantly and therefore the relic electrons have been injected since the birth of the pulsar. An additional caveat of this model is, that while this model is time-dependent and cooling losses of electrons were taken into account, no spatial dependence or spatial evolution was added. The fit of the different electron generations was performed only on the SED.

The free parameters of the model were fitted to the observation data from the X-ray satellites, as well as the SEDs derived in the joint-analysis from components A and B. For optimization, the package EMCEE (Foreman-Mackey et al. 2013), which applies a Markov chain Monte Carlo (MCMC) method, was used.

The optimised model is shown in Fig. 7, while the parameter space that yields the highest numerical probability (16–84% range of the probability distribution) for the MCMC method is listed in Table 6. This morphology independent model can account well for the observed emission.

The observed X-ray data in the vicinity of PSR J1813–1749 can be well explained by young electrons that were produced only in the last 0.4 kyr, the electrons of the third generation. The electrons produced in the last 2.7 kyr, the second generation, were observed by H.E.S.S. as compact emission around the pulsar. Taking into account the low statistics at high energies, the population of relic electrons from the first generation, released from the birth of the pulsar up to an age of 1.6 kyrs, can describe
<table><caption>**Table 6.** Validity range of the parameters used in the evolution of the leptonic model.</caption><thead><tr><th>Parameter</th><th>Validity range</th></tr></thead><tbody><tr><td>B<sub>now</sub></td><td>[11.36–13.18] µG</td></tr><tr><td>P<sub>0</sub></td><td>[20.19 – 22.97] ms</td></tr><tr><td>ϑ</td><td>[0.15–0.17]</td></tr><tr><td>α</td><td>[2.35–2.40]</td></tr><tr><td>Δt<sub>X-ray</sub></td><td>[0.88–0.91]</td></tr><tr><td>Δt<sub>PWN</sub></td><td>[0.35–0.39]</td></tr></tbody></table>

**Notes.** The parameter combination presented here should be interpreted as a range of possible combinations, not best-fit values.

the extended emission from HESS J1813–178B reasonably well but fails to describe the emission from component C. To describe this emission in a leptonic scenario, a second electron population and an unreasonably high particle density, only experienced by the second electron generation, would need to be introduced.

Using the parameters of the model, the true age of the system was estimated, following Gaensler & Slane (2006):

$$t_{\text{true}} = \frac{P}{(n-1)\dot{P}} \left(1 - \left(\frac{P_0}{P}\right)^{n-1}\right) \qquad (9)$$

with *n* the braking index of the pulsar. For the validity range of parameters in this analysis, the true age of the system was estimated to be ~4.7 kyrs. With this information, together with the measurement of the extension of 4FGL J1813–1737e and HESS J1813–178B in the energy bands specified in Table 3, an estimation of the diffusion coefficient and diffusion index can be made by comparing the observed angular size to that expected for the relic electron population at different γ-ray energies. In order to estimate a diffusion coefficient radial symmetric diffusion was assumed. To get a estimation of the radial symmetric extend of the emission from HESS J1813–178B, a symmetric Gaussian model was fitted to the emission in the energy bands defined in Sect. 3.3. This derived extension, as well as the age of the system, derived using Eq. (9), were then used to estimate the distance which the relic electrons have diffused since they were injected, with the expected angular size calculated as $r_{\text{diff}} = \sqrt{2t_{\text{age}} \cdot D}$, with *D* the diffusion coefficient. The diffusion coefficient was derived using the electron energy $E_e$, the diffusion index $\delta$ and the diffusion coefficient $D_0$ at the reference energy of 1 TeV:

$$D = D_0 \left( \frac{E_e}{1 \text{ GeV}} \right)^{\delta}. \qquad (10)$$

The diffusion coefficient is estimated to $D_0 = (6.98 \pm 0.69) \times 10^{28} \text{ cm}^2 \text{s}^{-1}$ with a diffusion index $\delta = (0.30 \pm 0.06)$. The estimated radius as a function of γ-ray energy, as well as the measurement, can be seen in Fig. 8. The estimated diffusion index lies in the range between Kraichnan turbulence ($\delta = 0.5$, Kraichnan 1965) and Kolmogorov turbulence ($\delta = 0.33$, Kolmogorov 1991), which is the canonical diffusion index typically assumed for PWNe. The best-fit diffusion coefficient is comparable with the diffusion coefficient estimated for the ISM (Strong & Moskalenko 1998).

To compare the results to the population of known pulsars, the energy density of both HESS J1813–178A and HESS J1813–178B in the TeV energy range was estimated, following the scheme in Giacinti et al. (2020). The energy density

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Measured extension of HESS J1813–178B in continuous energy bands, specified in Table 3, is shown in red. The figure also shows the expected $1\sigma$ containment radius for different diffusion coefficients.</figcaption></figure>

<table><caption>**Table 7.** Estimated electron density following Eq. (11).</caption><thead><tr><th>ε<sub>e</sub> [eVcm<sup>-3</sup>]</th><th>HESS J1813–178A</th><th>HESS J1813–178B</th></tr></thead><tbody><tr><td>Estimated from E&#x0307;</td><td>65.29 ± 4.49</td><td>0.11 ± 0.02</td></tr><tr><td>Estimated from L<sub>γ</sub></td><td>23.95 ± 1.64</td><td>0.04 ± 0.01</td></tr></tbody></table>

of the electrons was estimated using $\epsilon_e = \frac{E_{\text{inj}}}{V}$, or using the properties of the pulsar

$$\epsilon_e = \frac{\dot{E} \tau_c}{4\pi(R_1^2 \cdot R_2)/3}, \qquad (11)$$

with $\dot{E}$ the present-day spin-down power, $\tau_c$ the present-day characteristic age, $R_1$ the extension of the semi-major axis of the ellipse and $R_2$ the size the semi-minor axis. The energy density using the total TeV γ-ray luminosity was estimated by integrating over the electron population inferred by our best-fit model. This study finds $E_{\text{inj}} = 3.62 \times 10^{48} \text{ erg}$. The derived electron densities are given in Table 7. The energy density calculated from the γ-ray luminosity is lower than that estimated from the properties of the pulsar. This effect was already observed by Giacinti et al. (2020), a possible reason could be that only the estimation using the γ-ray luminosity takes the evolution into account, while the estimation using the current pulsar properties cannot sufficiently describe this evolution. The energy density calculated for HESS J1813–178A is compatible with the energy densities estimated for other PWNe, while the energy density of HESS J1813–178B is close to 0.1 eV cm⁻³ and therefore comparable to the values estimated for pulsar halos (Giacinti et al. 2020).

## 4.2. Hadronic origin

Due to the positional coincidence of the emission and the SNRs G012.8–00.0 and G012.7–00.0 or the star-forming region W33, there is also the possibility that the observed emission is a superposition of γ-ray emission produced by different sources that overlap along the line of sight. While there is no information available for SNR G012.7–00.0, the distance to PSR J1813–1749 and the SNR G012.8–00.0 is estimated to be
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Significance map computed from the H.E.S.S data with an correlation radius of 0.06°. The estimated position and density of molecular clouds in the region are indicated by the dashed lines. Additionally, the morphology of the best-fit model derived in the analysis of the Fermi-LAT and H.E.S.S. data is indicated. Only clouds with a distance between 4 kpc and 12 kpc and a large positional overlap with the extended emission observed in Fermi-LAT and H.E.S.S. are depicted in the counts map.</figcaption></figure>

6.2–12 kpc (Camilo et al. 2021), while the stellar cluster is expected to be located at a much smaller distance of 4.8 kpc (Messineo et al. 2011). It is therefore possible that the emission observed as HESS J1813–178A is a PWN powered by PSR J1813–1749, while the extended emission observed in the Fermi-LAT and H.E.S.S. data is caused by protons accelerated at the SNR shock front or the stellar cluster.

In order to produce high energy γ rays via hadronic interactions, target material is necessary. Using the molecular cloud catalogues provided in Rice et al. (2016) and Miville-Deschênes et al. (2017), three molecular clouds, positionally coincident with the best-fit position of the emission from HESS J1813–178B were identified. Their location compared to the emission observed in the analysis of the Fermi-LAT data can be seen in Fig. 9. In order to investigate the feasibility of this scenario as the origin of the emission observed from 4FGL J1813–1737e and HESS J1813–178, a proton and electron population, accelerated at the shock front of the SNR were evolved over time. Following Diesing & Caprioli (2019), the protons were injected with:

$$f(p) = N_p \cdot E^{-\alpha_p} \cdot \exp\left(\frac{-E}{E_{\text{max}}}\right) \cdot \exp\left(-\frac{E_{\text{min}}}{E}\right) \quad (12)$$

where $N_p$ is the normalisation factor, $E$ is the energy of the injected particles, $E_{\text{min}} = m_p$, the rest mass of the proton, $E_{\text{max}}$ is assumed to be 1 PeV, and $\alpha_p$ is the spectral index of the injected protons. The injection spectrum of the electrons was defined as:

$$f_e(p) = N_p \cdot k_{ep} \cdot E^{-(\alpha_p + \Delta\alpha)} \cdot \exp\left(\frac{E}{E_{\text{max}}}\right) \cdot \exp\left(-\frac{E_{\text{min}}}{E}\right) \quad (13)$$

with $\Delta\alpha$ the difference between the injection index of the electrons and protons, $k_{ep}$ the electron-to-proton ratio, and $E_{\text{min}} = 100m_e$, with the electron rest mass, as well as $E_{\text{max}} = 1$ TeV. This particle population is evolved in an environment with the ambient magnetic field $B_{\text{now}}$ and the particle density $d$. We assume

<table><caption>**Table 8.** Validity range of the free parameters of the hadronic model.</caption><tr><td>α<sub>p</sub></td><td>[2.01–2.08]</td></tr><tr><td>log<sub>10</sub>(E<sub>p</sub>)</td><td>[48.49–49.65] log (erg)</td></tr><tr><td>Δα</td><td>−[0.15–0.27]</td></tr><tr><td>k<sub>ep</sub></td><td>−[1.72–6.05]</td></tr><tr><td>B<sub>now</sub></td><td>[73–271] μG</td></tr></table>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Estimated γ-ray flux for the hadronic model compared to the γ-ray flux expected for the leptonic model computed in Sect. 4.1. The deviation of the SED to the respective model curves are shown in the bottom panel.</figcaption></figure>

a particle density of $d = 60 \text{ cm}^{-3}$, which corresponds to the density of one of the molecular clouds in the region (see Fig. 9).

We evolve this particle population with time following the approach in Sect. 4.1, using a MCMC-Chain. The fit parameters of the model, as well as the adjusted parameter range can be seen in Table 8. A comparison of the leptonic model derived in Sect. 4.1 and the hadronic model can be seen in Fig. 10, while Fig. E.4 shows the comparison of both models with the SED derived in this analysis, as well as the spectra derived from the emission observed in HAWC (Albert et al. 2020) and LHAASO (Cao et al. 2024).

Both models exhibit good agreement with the SED from component B as derived in Sect. 3.4, while only the leptonic model shows good agreement with the spectrum of 1LHAASO J1814–1719u* measured by KM2A. The spectrum of 3HWC J1813–174 does not show good agreement with either model.

Additionally to the discrepancy with the spectra derived by HAWC and LHAASO, the hadronic model has further caveats. The estimated distance between the molecular cloud at 4.6 kpc and the SNR at 6.2–12 kpc are not compatible, and the location of the molecular cloud can account for the emission observed in the GeV energy range, but not the extension of the TeV emission. Assuming a particle density of $d = 1 \text{ cm}^{-3}$, corresponding to ISM level, however, results in the need for a very powerful SNR with the energy of the protons required to be $\log_{10}(E_p/\text{erg}) = [51.2–51.7]$, to describe the emission well. Since there is no evidence supporting this scenario, an origin of the observed emission only from acceleration of protons on the SNR seems unlikely.
While such a high energy output seems unlikely in the case of only one object, collective effects of a population of young stars in the stellar cluster CL J1813-178 (located at a distance of 4.8 kpc; Messineo et al. 2011) could account for the required energy. For example, for a kinetic luminosity of $3 \times 10^{38}$ erg s$^{-1}$ and an age of 3–10 Myr, as may be typical for a young, massive stellar cluster, the available total energy budget is $\log_{10}(E_p/\text{erg}) = [52.5-53]$ (Morlino et al. 2021). Acceleration of protons in the stellar cluster remains therefore a valid possible origin for the observed extended $\gamma$-ray emission, a detailed examination is, however, beyond the scope of this paper.

## 5. Discussion

In previous analyses of the region around PSR J1813-1749 the analysis of H.E.S.S. data revealed compact emission with an extension of 0.06° around the pulsar (Aharonian et al. 2006b), while the analysis of Fermi-LAT data, as well as data acquired by HAWC and LHAASO showed largely extended emission. With these results, no connection between the Fermi-LAT source 4FGL J1813-1737e and the H.E.S.S. source HESS J1813-178 could be established. We reanalysed the region using an improved reconstruction for the TeV data and increased exposure in the GeV energy range.

In the TeV energy range, we confirm the detection of HESS J1813-178, a bright, TeV $\gamma$-ray source, centred at the position of the PSR J1813-1749, with an extension of 0.06°, observed as HESS J1813-178A in this analysis. The results derived for this emission component are consistent with the results derived in H.E.S.S. Collaboration (2018). We also detect a fainter $\gamma$-ray structure, with an extension of 0.7° enclosing the pulsar and HESS J1813-178A, which we refer to as HESS J1813-178B.

In the GeV energy range observed with Fermi-LAT, we find $\gamma$-ray emission with an extension of 0.4°, positionally coincident with the extended emission observed in the H.E.S.S. dataset. Additionally, we observe compact emission positionally coincident with the pulsar. This emission has already been reported by Araya (2018), but was not significant. While we also cannot claim a significant detection, adding component C improves the region's description.

By combining the datasets in a joint analysis we find that the extended emission detected in the H.E.S.S. and Fermi-LAT data can be connected and described by a single source model, while the emission from HESS J1813-178A can only be observed in the H.E.S.S. data, since the emission drops below the sensitivity of the detector in the Fermi-LAT energy range. We can therefore establish a consistent description of the region over five decades of energy and conclude that the emission observed by Fermi-LAT and H.E.S.S. have the same origin.

Araya (2018) concluded that the emission was most likely of hadronic origin due to the spectral shape. While Araya (2018) only took $\gamma$ rays in the GeV energy range into account, in this work, a combined model of GeV and TeV photons was used, as well as X-ray data included. We computed a physical model that can describe the observed $\gamma$-ray emission by evolving an electron population, or alternatively an proton population, over time. We find that the leptonic emission scenario describes the combined data well, while a acceleration at the SNR shock front seems less justified since the present target material cannot explain the observed extension of the TeV emission. The energetics necessary to produce the observed $\gamma$-ray emission could also be

explained by acceleration of protons inside the stellar cluster, but has not been further investigated in this study.

We find, that the data from HESS J1813-178A, as well as HESS J1813-178B can be described well by $\gamma$-ray emission from synchrotron and IC emission of electrons originating from the pulsar. This indicates that HESS J1813-178A, is likely a PWN, while the extended emission is possibly caused by electrons and positrons escaping the confines of the PWN and diffusing into the ISM. The effect of constructing a model that considers only time dependence, but not spatial dependence, leads to an over-prediction in the keV energy range. Since the measurement of the flux observed by INTEGRAL and XMM-Newton only includes a small area (with an extension of less than 80′′), while the measurement of the flux from HESS J1813-178B was conducted in an area of ~0.5°, we can assume that a low surface brightness, as well as the limitations of the detectors cause the difference between predicted flux and measured flux in the X-ray energy range.

Diffuse leptonic emission around a pulsar is usually only observed in systems older than ~10 kyr (Giacinti et al. 2020). The true age estimated in this study, as well as previous observations of PSR J1813-1749 (Dzib & Rodríguez 2021), indicate that the system is younger than 5 kyrs, but it is not implausible that environmental factors, such as the local ISM density distribution or the pulsar proper motion, might lead to particles diffusing outwards faster than in previously observed systems.

The estimated diffusion coefficient for HESS J1813-178B is within the expected range of diffusion in the ISM (Strong & Moskalenko 1998) and comparable to the diffusion coefficient estimated for the $\gamma$-ray emission observed around PSR B1823-13, the pulsar powering the highly extended PWN HESS J1825-137 (H.E.S.S. Collaboration 2019). Assuming the same diffusion coefficient for HESS J1813-178A, would lead to a extension far bigger than the observed emission, implying that the diffusion coefficient is non uniform, which is expected for a pulsar environment. Additionally, we find that the energy density estimated for the extended emission is comparable to the densities estimated for the halos observed around the evolved systems around PSR J0633+1746 and PSR B0656+14 (Giacinti et al. 2020). While we would expect the diffusion coefficient in a TeV halo to be below the diffusion observed in the ISM, the low energy density and lack of other accelerators in the area suggest that the system around PSR J1813-1749 shows characteristics of an older, evolved system despite its young age. Although it would be unusual to note electron escape from young PWNe, this may occur if the system becomes disrupted at an early stage, such as due to an early return of the supernova reverse shock preferentially from one direction.

This study also investigates the possibility of a hadronic origin of the emission, by evolving a proton population originating from SNR G012.8-00.0. In order to estimate the statistical description of both models, the absolute chi-squared ($\chi^2$) goodness-of-fit was estimated using the SED of component B derived in the joint-fit. This study finds $\chi^2_{\text{had}} = 18.56$ for the hadronic model and $\chi^2_{\text{lep}} = 30.62$ for the leptonic model. The hadronic model is therefore preferred at a $3.47\sigma$ level. Whilst the hadronic model yields a marginally better statistical description for the emission, the target material identified in the region is not sufficient to account for the extent of the detected $\gamma$-ray emission and target material with a density of the level of the ISM would require an unreasonably high proton energy to explain the detected $\gamma$-ray flux.
For the data below 10 GeV, observed as component C, we
cannot find evidence that convincing points towards either
a leptonic or hadronic origin. A hadronic scenario cannot be
supported since no association to a nearby interstellar cloud or
target material can be made. A leptonic scenario would imply
the existence of a second electron population, or an unreasonably
high particle density, which is only experienced by this second
electron population. Another possible leptonic scenario might be
a pulsar origin, which may alternatively account for the emission
component, but this study cannot find convincing evidence for
such an explanation. For a firm identification of the soft γ-ray
emission as originating from the pulsar, a dedicated study of the
emission at lower energies needs to be performed.

While we are not able to fully resolve the origin of the γ-
ray emission in the region around PSR J1813−1749, we present
a detailed morphological and spectral analysis of the region, as
well as possible emission scenarios. To further add to the under-
standing about this source, more observations at high-energy γ
rays are necessary, especially the addition of further data taken
by the Water Cherenkov detector arrays HAWC or LHAASO into
the joint-fit and the addition of further data taken by Fermi-LAT
will be of great value in constraining the possible origin of the
different components.

Acknowledgements. The support of the Namibian authorities and of the Uni-
versity of Namibia in facilitating the construction and operation of H.E.S.S. is
gratefully acknowledged, as is the support by the German Ministry for Edu-
cation and Research (BMBF), the Max Planck Society, the German Research
Foundation (DFG), the Helmholtz Association, the Alexander von Humboldt
Foundation, the French Ministry of Higher Education, Research and Inno-
vation, the Centre National de la Recherche Scientifique (CNRS/IN2P3 and
CNRS/INSU), the Commissariat à l'Énergie atomique et aux Énergies alterna-
tives (CEA), the U.K. Science and Technology Facilities Council (STFC), the
Irish Research Council (IRC) and the Science Foundation Ireland (SFI), the Knut
and Alice Wallenberg Foundation, the Polish Ministry of Education and Sci-
ence, agreement no. 2021/WK/06, the South African Department of Science and
Technology and National Research Foundation, the University of Namibia, the
National Commission on Research, Science & Technology of Namibia (NCRST),
the Austrian Federal Ministry of Education, Science and Research and the Aus-
trian Science Fund (FWF), the Australian Research Council (ARC), the Japan
Society for the Promotion of Science, the University of Amsterdam and the Sci-
ence Committee of Armenia grant 21AG-1C085. We appreciate the excellent
work of the technical support staff in Berlin, Zeuthen, Heidelberg, Palaiseau,
Paris, Saclay, Tübingen and in Namibia in the construction and operation of
the equipment. This work benefited from services provided by the H.E.S.S.
Virtual Organisation, supported by the national resource providers of the EGI
Federation.

References

Abdalla, H., Aharonian, F., Ait Benkhali, F., et al. 2021, ApJ, 917, 6

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33

Acero, F., Aharonian, F., Akhperjanian, A. G., et al. 2010, MNRAS, 402, 1877

Aharonian, F., Akhperjanian, A. G., Aye, K. M., et al. 2005, Science, 307, 1938

Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006a, A&A, 457, 899

Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006b, ApJ, 636, 777

Ajello, M., Atwood, W. B., Axelsson, M., et al. 2021, ApJS, 256, 12

Albert, J., Aliu, E., Anderhub, H., et al. 2006, ApJ, 637, L41

Albert, A., Alfaro, R., Alvarez, C., et al. 2020, ApJ, 905, 76

Araya, M. 2018, ApJ, 859, 69

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071

Bernlöhr, K. 2008, Astropart. Phys., 30, 149

Brogan, C. L., Gaensler, B. M., Gelfand, J. D., et al. 2005a, ApJ, 629, L105

Camilo, F., Ransom, S. M., Halpern, J. P., & Roshi, D. A. 2021, ApJ, 917, 67

Cao, Z., Aharonian, F., An, Q., et al. 2024, ApJS, 271, 25

Deil, C., Donath, A., Terrier, R., et al. 2020, https://doi.org/10.5281/
zenodo.4701500

Diesing, R., & Caprioli, D. 2019, Phys. Rev. Lett., 123, 071101

Donath, A., Terrier, R., Remy, Q., et al. 2023, A&A, 678, A157

Dzib, S. A., & Rodríguez, L. F. 2021, ApJ, 923, 228

Fermi-LAT Collaboration 2019, https://fermi.gsfc.nasa.gov/ssc/data/
analysis/scitools/Aeff_Systematics.html

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Funk, S., Hinton, J. A., Moriguchi, Y., et al. 2007, A&A, 470, 249

Gaensler, B. M., & Slane, P. O. 2006, ARA&A, 44, 17

Giacinti, G., Mitchell, A. M. W., López-Coto, R., et al. 2020, A&A, 636, A113

Gotthelf, E. V., & Halpern, J. P. 2009, ApJ, 700, L158

Hahn, J. 2015, Int. Cosmic Ray Conf., 34, 917

Helfand, D. J., Gotthelf, E. V., Halpern, J. P., et al. 2007, ApJ, 665, 1297

H.E.S.S. Collaboration 2018, A&A, 612, A1

H.E.S.S. Collaboration 2019, A&A, 621, A116

H.E.S.S. Collaboration 2023, A&A, 672, A103

H. E. S. S. Collaboration (Aharonian, F., et al.) 2023, A&A, 672, A103

Khelifi, B., Djannati-Ataï, A., Jouvin, L., et al. 2015, Int. Cosmic Ray Conf., 34, 837

Kolmogorov, A. N. 1991, Proc. R. Soc. London Ser. A, 434, 9

Kraichnan, R. H. 1965, Phys. Fluids, 8, 1385

Li, T. P., & Ma, Y. Q. 1983, ApJ, 272, 317

MAGIC Collaboration 2020, MNRAS, 497, 3734

Messineo, M., Davies, B., Figer, D. F., et al. 2011, ApJ, 733, 41

Miville-Deschênes, M.-A., Murray, N., & Lee, E. J. 2017, ApJ, 834, 57

Mohrmann, L., Specovius, A., Tiziani, D., et al. 2019, A&A, 632, A72

Morlino, G., Blasi, P., Peretti, E., & Cristofari, P. 2021, MNRAS, 504, 6096

Parsons, R. D., & Hinton, J. A. 2014, Astropart. Phys., 56, 26

Principe, G., Mitchell, A. M. W., Caroff, S., et al. 2020, A&A, 640, A76

Rice, T. S., Goodman, A. A., Bergin, E. A., Beaumont, C., & Dame, T. M. 2016, ApJ, 822, 52

Strong, A. W., & Moskalenko, I. V. 1998, ApJ, 509, 212

Ubertini, P., Bassani, L., Malizia, A., et al. 2005, ApJ, 629, L109

Venter, C., & de Jager, O. C. 2007, in WE-Heraeus Seminar on Neutron Stars and Pulsars 40 years after the Discovery, eds. W. Becker, & H. H. Huang, 40

Wood, M., Caputo, R., Charles, E., et al. 2017, Int. Cosmic Ray Conf., 301, 824

Xin, Y., & Guo, X. 2021, Proc. Int. Cosmic Ray Conf., 395, 625

<footnote n="1">¹ Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2, Ireland</footnote>

<footnote n="2">² Max-Planck-Institut für Kernphysik, PO Box 103980, D 69029 Heidelberg, Germany</footnote>

<footnote n="3">³ Yerevan Physics Institute, 2 Alikhanian Brothers St., Yerevan 0036, Armenia</footnote>

<footnote n="4">⁴ Yerevan State University, 1 Alek Manukyan St, Yerevan 0025, Armenia</footnote>

<footnote n="5">5 Landessternwarte, Universität Heidelberg, Königstuhl, 69117 Heidelberg, Germany</footnote>

<footnote n="6">⁶ Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands</footnote>

<footnote n="7">7 Laboratoire Leprince-Ringuet, École Polytechnique, CNRS, Institut Polytechnique de Paris, 91128 Palaiseau, France</footnote>

<footnote n="8">8 University of Namibia, Department of Physics, Private Bag 13301, Windhoek 10005, Namibia</footnote>

<footnote n="9">9 Centre for Space Research, North-West University, Potchefstroom 2520, South Africa</footnote>

<footnote n="10">10 Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149, 22761 Hamburg, Germany</footnote>

<footnote n="11">11 Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen, Germany</footnote>

<footnote n="12">12 Department of Physics, University of the Free State, PO Box 339, Bloemfontein 9300, South Africa</footnote>

<footnote n="13">13 Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Strasse 24/25, 14476 Potsdam, Germany</footnote>

<footnote n="14">14 Université Paris Cité, CNRS, Astroparticule et Cosmologie, 75013 Paris, France</footnote>

<footnote n="15">15 Department of Physics and Electrical Engineering, Linnaeus University, 351 95 Växjö, Sweden</footnote>

<footnote n="16">16 Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany</footnote>

<footnote n="17">17 Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, 72076 Tübingen, Germany</footnote>
18 Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université Paris Cité, 5 Pl. Jules Janssen, 92190 Meudon, France

19 Sorbonne Université, Université Paris Diderot, Sorbonne Paris Cité, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE, 4 Place Jussieu, 75252 Paris, France

20 IRFU, CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette, France

21 Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Nikolaus-Fiebiger-Str. 2, 91058 Erlangen, Germany

22 Astronomical Observatory, The University of Warsaw, Al. Ujazdowskie 4, 00-478 Warsaw, Poland

23 Université Savoie Mont Blanc, CNRS, Laboratoire d'Annecy de Physique des Particules - IN2P3, 74000 Annecy, France

24 Instytut Fizyki Jądrowej PAN, ul. Radzikowskiego 152, 31-342 Kraków, Poland

25 Université Bordeaux, CNRS, LP2I Bordeaux, UMR 5797, 33170 Gradignan, France

26 School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue, Braamfontein, Johannesburg, 2050 South Africa

27 University of Oxford, Department of Physics, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK

28 School of Physical Sciences, University of Adelaide, Adelaide 5005, Australia

29 Aix Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

30 Laboratoire Univers et Particules de Montpellier, Université Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, 34095 Montpellier Cedex 5, France

31 Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstraße 25, 6020 Innsbruck, Austria

32 Obserwatorium Astronomiczne, Uniwersytet Jagielloński, ul. Orla 171, 30-244 Kraków, Poland

33 Institute of Astronomy, Faculty of Physics, Astronomy and Informationics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland

34 Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18, 00-716 Warsaw, Poland

35 Department of Physics and Astronomy, The University of Leicester, University Road, Leicester, LE1 7RH, UK

36 GRAPPA, Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

37 Department of Physics, Konan University, 8-9-1 Okamoto, Higashinada, Kobe, Hyogo 658-8501, Japan

38 Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS), The University of Tokyo, 5-1-5 Kashiwa-no-Ha, Kashiwa, Chiba, 277-8583, Japan

39 RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan

40 School of Science, Western Sydney University, Locked Bag 1797, Penrith South DC, NSW 2751, Australia
# Appendix A: Estimation of the hadronic background

The hadronic background was estimated using a template model constructed from archival H.E.S.S. data (see Mohrmann et al. (2019)). To account for the variation in cosmic-ray flux for the respective observation conditions, a 3D likelihood fit of this template background model for every observation was performed, using the background amplitude $\Phi$ and spectral index $\delta$ as fit parameters. Using the spectral index small inacurracies in the spectral shape of the background model can be corrected by modifying the predicted background rate $R_{BG}^* = R_{BG} \cdot (E/E_0)^{-\delta_{BG}}$ at Energy $E$, with a reference energy $E_0 = 1$ TeV.

To get an accurate estimation of the background counts the regions in which significant emission is detected needs to be excluded from the background fit. These excluded regions are shown in Fig. A.1. Figure A.2 shows the distribution of the background amplitude and spectral index, as well as a Gaussian fit indicated in red, while Fig. A.3 shows the significance distribution outside of the exclusion regions. This distribution follows a Gaussian distribution with a mean close two zero and a width that is close to one, indicating that only statistical fluctuations are present outside of the exclusion regions.

<figure><img src="image_13.png" /><figcaption>Fig. A.3. Significance distribution of the data outside of the defined exclusion regions. The orange curve shows the expected behaviour for a normal distribution, indicating that only statistical fluctuations are present, while the red curve shows a fit of a Gaussian distribution to the data.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. A.1. Significance map of the region around HESS J1813–178, as seen after the fit of the background model. The exclusion regions used for the fitting are indicated by the dashed lines. A correlation radius of 0.4° is used for the computation of this map.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. A.2. Distribution of the fit parameters of the template background model and a Gaussian fit to the respective distribution. The mean and standard deviation of the background spectral index $\delta$ and amplitude $\Phi$ are indicated.</figcaption></figure>

# Appendix B: Calculation of systematic uncertainties

In order to derive the systematic uncertainties introduced by a mismodeling of the energy axis of the IRFs, a randomly generated value is drawn from a Gaussian distribution with a mean of 1 and a standard deviation of 10%. We apply this factor $\rho_E$ to the energy axes of the IRFs.

To estimate the uncertainties introduced through mismodeling of the hadronic background, a linear gradient with a direction angle $\alpha_{bkg}$ and a gradient amplitude $\Phi_{bkg}$ is applied to the background model. In addition to the gradient parameters, the amplitude $\Phi$ and index $\delta$ of the model is varied by randomly drawing a value from Gaussian distributions. The mean and deviation of these distributions are indicated in table B.1. The standard deviation of these distributions was chosen such that they represent the deviation between the mean value estimated in the pre-fit of the background model (described in Appendix A) and the values estimated in the likelihood minimisation on the whole dataset.

Then, the shift, corresponding to the randomly drawn factor, was applied to the IRFs of the dataset and the background gradient introduced. The best-fit model, described in section 3.1, was applied, and the likelihood fit executed. This process was repeated 500 times.

An example of the resulting distribution of fit parameters can be seen in Fig. B.1, where the statistical error is indicated by the orange band and the orange dashed line represents the best-fit value derived in the analysis of the original dataset. The standard

<table><caption>Table B.1. Parameter distribution used for the computation of the systematic errors.</caption><thead><tr><th>parameter</th><th>distribution</th><th>variation</th></tr></thead><tbody><tr><td>ρ<sub>E</sub></td><td>Gaussian</td><td>μ = 1, σ = 0.1</td></tr><tr><td>Φ</td><td>Gaussian</td><td>μ = 1, σ = 0.05</td></tr><tr><td>δ</td><td>Gaussian</td><td>μ = 0, σ = 0.05</td></tr><tr><td>α<sub>bkg</sub></td><td>uniform</td><td>0° – 360°</td></tr><tr><td>Φ<sub>bkg</sub></td><td>Gaussian</td><td>μ = 1, σ = 0.01</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>**Fig. B.1.** Distribution of the best-fit values of the spectral index for HESS J1813−178B. A Gaussian fit to the data is indicated in red. The orange dashed line corresponds to the fit value derived in the analysis of the original dataset, the orange band indicates the statistical error.</figcaption></figure>

deviation of this Gaussian distribution represents the systematic error derived by this method. The small systematic uncertainties estimated for the spectral index and morphological parameters indicate that this method does not take into account all possible sources of systematic uncertainties. A more general estimation of uncertainties on the spectral index can be found in Aharonian et al. (2006a).

## Appendix C: Coverage of the region using H.E.S.S.

To visualize the coverage of the region using H.E.S.S., an exposure map (figure C.2) was computed. Additionally, a significance map computed with a small correlation radius of 0.06°, which is comparable to the PSF of the analysis, is depicted in figure C.1.

<figure><img src="image_17.png" /><figcaption>**Fig. C.1.** Significance map of the ROI, computed with a correlation radius of 0.06°.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. C.2.** Exposure map of the dataset used for the analysis of the H.E.S.S. data. Most of the exposure comes from observations centred on the neighbouring source HESS J1809-193.</figcaption></figure>

## Appendix D: Fermi-LAT data: Fermipy-Gammapy cross check

Whilst *gammapy* cannot be used to directly analyse *Fermi*-LAT data, the *fermipy* analysis package is not able to perform a three-dimensional fit to the data. To combine both methods, the count's cube, the IRFs and the background model cubes were computed using the *fermipy* package. Then, the cubes and background source models were extracted and converted into a format that is supported by *gammapy*. The robustness of exporting the data and models to *gammapy* needs to be validated before the analysis results are used for further investigation. For this purpose, a simple analysis of the source region in *fermipy* and *gammapy* is carried out. The morphology of the source given in the 4FGL-catalogue, a disc model with a source position of R.A. = 273.405°, Dec. = −17.653° and an extension of 0.6°, was used for this validation. Then, the ROI was optimised with *fermipy* by only taking into account sources with a test statistics (TS) greater than four and predicted counts of more than 50. These cuts yielded a total of 24 source models present in the analysis region. The spectral parameters of these sources in a 5° radius around the catalogue position of 4FGL J1813−1737e were refitted.

<figure><img src="image_19.png" /><figcaption>**Fig. D.1.** Best-fit spectrum of a power-law model applied to the *Fermi*-LAT data in *fermipy* and *gammapy*.</figcaption></figure>
For the analysis in **gammapy**, the morphology of the source was fixed to the position derived in **fermipy** and the best-fit spectral parameters were evaluated. The spectral parameters and SED derived in **gammapy** and **fermipy** agree within errors. In addition to the spectral fit parameters, a comparison of the counts spectra of the background models and other source models in the ROI, can be seen in Appendix, Fig. D.1. We find that the results from both likelihood minimisations are in agreement with each other.

## Appendix E: Multiwavelength context

To better understand the nature of the emission in the region around PSR J1813-1749, a comparison of the observed emission in different energy ranges was made. Comparing the best-fit model derived in the analysis of the *Fermi*-LAT data and the H.E.S.S. data yields good spatial agreement (Fig. E.1). The elongation is in both cases observed along the galactic plane and the pulsar is contained.

<figure><img src="image_20.png" /></figure>

**Fig. E.1.** Excess map of the H.E.S.S. data after the addition of the Gaussian and elongated Gaussian model. The position of PSR J1813-1749 is indicated by the blue dot. The best-fit model describing the emission observed by H.E.S.S. is indicated by the black dashed lines. The *Fermi*-LAT best-fit model is depicted by the green ellipse.

<figure><img src="image_21.png" /></figure>

**Fig. E.2.** Zoom-in on a significance map with a correlation radius of 0.06° of the region around PSR J1813-1749 as seen with H.E.S.S.. The best-fit morphology of HESS J1813-178A and HESS J1813-178B are indicated by the white lines, the black contours depict SNR G012.8-00.0 and CL J1813-178 as observed by VLA. The position of PSR J1813-1749 is indicated in blue.

Taking into account data taken by the VLA (red contours in Fig. E.2), a positional coincidence between HESS J1813-178A, the pulsar and SNR can be observed. Additionally shown in the figure is a part of the stellar cluster CL J1813-178 located in close vicinity to the pulsar. The region has also been observed by HAWC (Albert et al. 2020) and LHAASO (Cao et al. 2024), both detecting extended emission around the position of HESS J1813-178A. In order to obtain a better understanding of the region, the position and extension of these sources, as well as the best-fit model obtained in this work, are visualised in a significance map computed from the H.E.S.S. data (see figure E.3). The

<figure><img src="image_22.png" /></figure>

**Fig. E.3.** Best-fit morphology derived from the analysis of the H.E.S.S. data is indicated in black. The best-fit morphology of the source 3HWC J1813-174 is indicated in blue. The description of the region obtained with LHAASO is indicated by the pink and orange dashed circles.

<figure><img src="image_23.png" /></figure>

**Fig. E.4.** Estimated γ-ray flux for the hadronic and leptonic model together with the spectra derived in the analysis of the Fermi-LAT, H.E.S.S., HAWC and LHAASO data.

morphology of the different components show a good spatial agreement. Additionally to the morphology, the spectra, derived in Albert et al. (2020) and Cao et al. (2024) are visualised in figure E.4, together with the leptonic and hadronic model derived in this work. The spectra of 1LHAASO J1814-1719u*, observed by KM2A and 3HWC J1813-174 show a good agreement with the results derived in this study, while the spectrum of 1LHAASO J1814-1719u*, observed by WCDA as well as the spectrum of 1LHAASO J1814-1636u do not agree. A possible reason for this could be the differences in the separation between source emission and diffuse emission along the galactic plane.