© The Authors 2024

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 to a large number of possible acceleration sites, a firm identification of VHE sources often remains elusive, with high-energy γ 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, γ 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 (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 γ-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 π₀ particles which then decay into two γ-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 interesting 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 γ-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 γ-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 crosschecked 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: (1)

with σ_M being the major axis of the Gaussian and Θ the angular distance to the model centre. In the case of an elongated Gaussian, the major axis σ_M is defined by σ_eff, the effective containment radius: (2)

with e being the eccentricity and ∆ϕ the difference between the position angle of the Gaussian ϕ and the position angle of the evaluation point.

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

where N₀ is the flux normalisation factor, Γ the spectral index and E₀ the reference energy, which is set to 1 TeV for the whole study. Alternatively, a logarithmic parabola spectral model was used: (4)

where β is the curvature, as well as a power law with a generalised exponential cut-off (5)

with λ, 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.

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 (gll_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 ( 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σ. 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 ATS = 59.91, while the elongated Gaussian model is preferred by 3.4σ 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σ. The best-fit spectral and spatial 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.

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.

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.

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, 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⁻¹² cm⁻² s⁻¹ TeV⁻¹ for 4FGL J1813–1737e and 1.79 × 10⁻¹² cm⁻² s⁻¹ TeV⁻¹ 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 a 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⁻¹² cm⁻² s⁻¹ TeV⁻¹ for 4FGL J1813–1737e and 0.36 × 10⁻¹² cm⁻² s⁻¹ TeV⁻¹ for component C. Both errors were then added in quadrature.

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.

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).

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 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 σ 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 . 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 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 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.

Fig. 5 Best-fit angular offset and area of the ellipse of HESS J1813–178A and HESS J1813–178B with 1σ uncertainties, as well as 4FGL J1813–1737e data in independent energy bins (see Table 3) H1 – H3 and F1 – F6.

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.

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 cosmicray 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 C1J1813-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 detected as extended high-energy γ-ray emission. The second electron generation, the middle-aged electrons, were injected since t_e = t_true · t_PWN, with t_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_true · t_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 α, a break energy of E_b and a cutoff E_cut: (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 Ṗ, and spin-down power 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 α, the birth period P₀, and the conversion fraction of the spin-down luminosity to electrons ϑ, as well as the magnetic field present at the time of the observation of the γ-ray emission B_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: (7)

with τ₀ = P₀/(2 Ṗ₀), the spin-down time at the birth of the pulsar and the initial spin-down energy loss rate at birth of the pulsar. Similarly the pulsar period was assumed to evolve as , with P₀ the birth period of the pulsar and the magnetic field as (Gaensler & Slane 2006; Venter & de Jager 2007) with B₀ the magnetic field at birth of the pulsar. The spin down luminosity is defined as: (8)

with the inital spin down luminosity and .

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 optimisation, 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 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): (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 , with D the diffusion coefficient. The diffusion coefficient was derived using the electron energy E_e, the diffusion index δ and the diffusion coefficient D₀ at the reference energy of 1 TeV: (10)

The diffusion coefficient is estimated to D₀ = (6.98 ± 0.69) × 10²⁸ cm² s⁻¹ with a diffusion index δ = (0.30 ± 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 (δ = 0.5, Kraichnan 1965) and Kolmogorov turbulence (δ = 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 of the electrons was estimated using , or using the properties of the pulsar (11)

with the present-day spin-down power, τ_c the present-day characteristic age, R₁ the extension of the semi-major axis of the ellipse and R₂ 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_inj = 3.62 × 10⁴⁸ 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).

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.

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σ containment radius for different diffusion coefficients.

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 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: (12)

where N_p is the normalisation factor, E is the energy of the injected particles, E_min = m_p, the rest mass of the proton, E_max is assumed to be 1 PeV, and α_p is the spectral index of the injected protons. The injection spectrum of the electrons was defined as: (13)

with ∆α the difference between the injection index of the electrons and protons, k_ep the electron-to-proton ratio, and E_min = 100m_e, with the electron rest mass, as well as E_max = 1 TeV. This particle population is evolved in an environment with the ambient magnetic field B_now and the particle density d. We assume a particle density of d = 60 cm⁻³, 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 cm⁻³, corresponding to ISM level, however, results in the need for a very powerful SNR with the energy of the protons required to be log₁₀(E_p/erg) = [51.2–51.7], to describe the emission well. Since there is no evidence supporting this scenario, a 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 × 10³⁸ 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₁₀(E_p/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 γ-ray emission, a detailed examination is, however, beyond the scope ofthis paper.

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.

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.

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 γ-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 γ-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 γ-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 γ 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 γ-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 γ-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 γ-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 & Rodriguez 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 γ-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 (χ²) goodness-of-fit was estimated using the SED of component B derived in the joint-fit. This study finds for the hadronic model and for the leptonic model. The hadronic model is therefore preferred at a 3.47σ 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 γ-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 γ-ray flux.

For the data below 10 GeV, observed as component C, we cannot find evidence that convincingly 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 understanding 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 University 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 Education 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 Innovation, the Centre National de la Recherche Scientifique (CNRS/IN2P3 and CNRS/INSU), the Commissariat à l’Énergie atomique et aux Énergies alternatives (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 Science, 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 Austrian Science Fund (FWF), the Australian Research Council (ARC), the Japan Society for the Promotion of Science, the University of Amsterdam and the Science Committee of Armenia grant 21AG-1C085. We appreciate the excellent work of the technical support staff in Berlin, Zeuthen, Heidelberg, Palaiseau, Paris, Saclay, Tubingen 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.

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 Φ and spectral index δ 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 at Energy E, with a reference energy E₀ = 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.

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.

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 δ and amplitude Φ are indicated.

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.

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 ϱ_E to the energy axes of the IRFs.

To estimate the uncertainties introduced through mis-modeling of the hadronic background, a linear gradient with a direction angle α_bk𝑔 and a gradient amplitude Φ_bk𝑔 is applied to the background model. In addition to the gradient parameters, the amplitude Φ and index δ 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 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).

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.

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

Fig. C.1 Signficiance map of the ROI, computed with a correlation radius of 0.06º.

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.

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 convert 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.

Fig. D.1 Best-fit spectrum of a power-law model applied to the Fermi-LAT data in fermipy and gammapy.

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.

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.

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 C1 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 visualise in a significance map computed from the H.E.S.S. data (see figure E.3). The 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.

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 indicate by the pink and orange dashed circles.

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.

References

All Tables

All Figures

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.

In the text

	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.
In the text

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.

In the text

	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).
In the text

	Fig. 5 Best-fit angular offset and area of the ellipse of HESS J1813–178A and HESS J1813–178B with 1σ uncertainties, as well as 4FGL J1813–1737e data in independent energy bins (see Table 3) H1 – H3 and F1 – F6.
In the text

	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.
In the text

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.

In the text

	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σ containment radius for different diffusion coefficients.
In the text

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.

In the text

	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.
In the text

	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.
In the text

	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 δ and amplitude Φ are indicated.
In the text

	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.
In the text

	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.
In the text

	Fig. C.1 Signficiance map of the ROI, computed with a correlation radius of 0.06º.
In the text

	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.
In the text

	Fig. D.1 Best-fit spectrum of a power-law model applied to the Fermi-LAT data in fermipy and gammapy.
In the text

	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.
In the text

	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.
In the text

	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 indicate by the pink and orange dashed circles.
In the text

	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.
In the text