# HESS J1809-193: a halo of escaped electrons around a pulsar wind nebula?

F Aharonian, F. Ait Benkhali, J Aschersleben, H Ashkar, M Backes, V

Barbosa Martins, R Batzofin, Y Becherini, D Berge, M Böttcher, et al.

► To cite this version:

F Aharonian, F. Ait Benkhali, J Aschersleben, H Ashkar, M Backes, et al., HESS J1809-193: a halo of escaped electrons around a pulsar wind nebula?. Astronomy and Astrophysics - A&A, 2023, 672, pp.A103. 10.1051/0004-6361/202245459. hal-04019936

HAL Id: hal-04019936

https://hal.science/hal-04019936

Submitted on 10 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
# HESS J1809-193: A halo of escaped electrons around a pulsar wind nebula?

H.E.S.S. Collaboration: F. Aharonian<sup>1,2</sup>, F. Ait Benkhali<sup>3</sup>, J. Aschersleben<sup>4</sup>, H. Ashkar<sup>5</sup>, M. Backes<sup>6,7</sup>, V. Barbosa Martins<sup>8</sup>, R. Batzofin<sup>9</sup>,

Y. Becherini<sup>10,11</sup>, D. Berge<sup>8,12</sup>, M. Böttcher<sup>7</sup>, C. Boisson<sup>13</sup>, J. Bolmont<sup>14</sup>, J. Borowska<sup>12</sup>, M. Bouyahiaoui<sup>2</sup>, F. Bradascio<sup>15</sup>, M. Breuhaus<sup>2</sup>,
R. Brose<sup>1</sup>, F. Brun<sup>15</sup>, B. Bruno<sup>16</sup>, T. Bulik<sup>17</sup>, C. Burger-Scheidlin<sup>1</sup>, T. Bylund<sup>11</sup>, S. Caroff<sup>18</sup>, S. Casanova<sup>19</sup>, J. Celic<sup>16</sup>, M. Cerruti<sup>10</sup>,

P. Chambery<sup>20</sup>, T. Chand<sup>7</sup>, A. Chen<sup>9</sup>, J. Chibueze<sup>7</sup>, O. Chibueze<sup>7</sup>, J. Damascene Mbarubucyeye<sup>8</sup>, A. Djannati-Ataï<sup>10</sup>, A. Dmytriiev<sup>7</sup>,
S. Einecke<sup>21</sup>, J.-P. Ernenwein<sup>22</sup>, K. Feijen<sup>21</sup>, M. Filipovic<sup>23</sup>, G. Fontaine<sup>5</sup>, M. Füßling<sup>8</sup>, S. Funk<sup>16</sup>, S. Gabici<sup>10</sup>, Y. A. Gallant<sup>24</sup>,

S. Ghafourizadeh<sup>3</sup>, G. Giavitto<sup>8</sup>, L. Giunti<sup>10,15</sup>, D. Glawion<sup>16</sup>, P. Goswami<sup>7</sup>, G. Grolleron<sup>14</sup>, M.-H. Grondin<sup>20</sup>, L. Haerer<sup>2</sup>, J. A. Hinton<sup>2</sup>,

W. Hofmann<sup>2</sup>, T. L. Holch<sup>8</sup>, M. Holler<sup>25</sup>, D. Horns<sup>26</sup>, Zhiqiu Huang<sup>2</sup>, M. Jamrozy<sup>27</sup>, F. Jankowsky<sup>3</sup>, V. Joshi<sup>16,*</sup> (©), I. Jung-Richardt<sup>16</sup>,

E. Kasai<sup>6</sup>, K. Katarzyński<sup>28</sup>, B. Khélifi<sup>10</sup>, W. Kluźniak<sup>29</sup>, Nu. Komin<sup>9</sup>, K. Kosack<sup>15</sup>, D. Kostunin<sup>8</sup>, R. G. Lang<sup>16</sup>, S. Le Stum<sup>22</sup>, F. Leitl<sup>16</sup>,

A. Lemière<sup>10</sup>, M. Lemoine-Goumand<sup>20</sup>, J.-P. Lenain<sup>14</sup>, F. Leuschner<sup>30</sup>, T. Lohse<sup>12</sup>, A. Luashvili<sup>13</sup>, I. Lypova<sup>3</sup>, J. Mackey<sup>1</sup>, D. Malyshev<sup>30</sup>,

D. Malyshev<sup>16</sup>, V. Marandon<sup>2</sup>, P. Marchegiani<sup>9</sup>, A. Marcowith<sup>24</sup>, P. Marinos<sup>21</sup>, G. Martí-Devesa<sup>25</sup>, R. Marx<sup>3</sup>, A. Mitchell<sup>16</sup>, R. Moderski<sup>29</sup>,

L. Mohrmann<sup>2,*</sup> (©), A. Montanari<sup>15</sup>, E. Moulin<sup>15</sup>, J. Muller<sup>5</sup>, K. Nakashima<sup>16</sup>, M. de Neurois<sup>5</sup>, J. Niemiec<sup>19</sup>, A. Priyana Noel<sup>27</sup>, S. Ohm<sup>8</sup>,

L. Olivera-Nieto<sup>2</sup>, E. de Ona Wilhelmi<sup>8</sup>, M. Ostrowski<sup>27</sup>, S. Panny<sup>25</sup>, M. Panter<sup>2</sup>, R. D. Parsons<sup>12</sup>, D. A. Prokhorov<sup>31</sup>, G. Pühlhofer<sup>30</sup>,

M. Punch<sup>10</sup>, A. Quirrenbach<sup>3</sup>, P. Reichherzer<sup>15</sup>, A. Reimer<sup>25</sup>, O. Reimer<sup>25</sup>, M. Renaud<sup>24</sup>, B. Reville<sup>2</sup>, F. Rieger<sup>2</sup>, G. Rowell<sup>21</sup>, B. Rudak<sup>29</sup>,

V. Sahakian<sup>32</sup>, A. Santangelo<sup>30</sup>, M. Sasaki<sup>16</sup>, H. M. Schutte<sup>7</sup>, U. Schwanke<sup>12</sup>, J. N. S. Shapopi<sup>6</sup>, H. Sol<sup>13</sup>, A. Specovius<sup>16</sup>, S. Spencer<sup>16</sup>,

Ł. Stawarz<sup>27</sup>, R. Steenkamp<sup>6</sup>, S. Steinmassl<sup>2</sup>, I. Sushch<sup>7</sup>, H. Suzuki<sup>33</sup>, T. Takahashi<sup>34</sup>, T. Tanaka<sup>33</sup>, R. Terrier<sup>10</sup>, C. Thorpe-Morgan<sup>30</sup>,

M. Tsirou<sup>2</sup>, N. Tsuji<sup>35</sup>, Y. Uchiyama<sup>36</sup>, C. van Eldik<sup>16</sup>, M. Vecchi<sup>4</sup>, J. Veh<sup>16</sup>, C. Venter<sup>7</sup>, J. Vink<sup>31</sup>, T. Wach<sup>16</sup>, S. J. Wagner<sup>3</sup>, R. White<sup>2</sup>,

A. Wierzcholska<sup>19</sup>, Yu Wun Wong<sup>16</sup>, M. Zacharias<sup>3,7</sup>, D. Zargaryan<sup>1</sup>, A. A. Zdziarski<sup>29</sup>, A. Zech<sup>13</sup>, S. Zouari<sup>10</sup>, and N. Żywucka<sup>7</sup>

*(Affiliations can be found after the references)*

Received 14 November 2022 / Accepted 24 February 2023

## ABSTRACT

**Context.** HESS J1809-193 is an unassociated very-high-energy γ-ray source located on the Galactic plane. While it has been connected to the nebula of the energetic pulsar PSR J1809-1917, supernova remnants and molecular clouds present in the vicinity also constitute possible associations. Recently, the detection of γ-ray emission up to energies of ~100 TeV with the HAWC observatory has led to renewed interest in HESS J1809-193.

**Aims.** We aim to understand the origin of the γ-ray emission of HESS J1809-193.

**Methods.** We analysed 93.2 h of data taken on HESS J1809-193 above 0.27 TeV with the High Energy Stereoscopic System (H.E.S.S.), using a multi-component, three-dimensional likelihood analysis. In addition, we provide a new analysis of 12.5 yr of Fermi-LAT data above 1 GeV within the region of HESS J1809-193. The obtained results are interpreted in a time-dependent modelling framework.

**Results.** For the first time, we were able to resolve the emission detected with H.E.S.S. into two components: an extended component (modelled as an elongated Gaussian with a 1-σ semi-major and semi-minor axis of ~0.62° and ~0.35°, respectively) that exhibits a spectral cutoff at ~13 TeV, and a compact component (modelled as a symmetric Gaussian with a 1-σ radius of ~0.1°) that is located close to PSR J1809-1917 and shows no clear spectral cutoff. The Fermi-LAT analysis also revealed extended γ-ray emission, on scales similar to that of the extended H.E.S.S. component.

**Conclusions.** Our modelling indicates that based on its spectrum and spatial extent, the extended H.E.S.S. component is likely caused by inverse Compton emission from old electrons that form a halo around the pulsar wind nebula. The compact component could be connected to either the pulsar wind nebula or the supernova remnant and molecular clouds. Due to its comparatively steep spectrum, modelling the Fermi-LAT emission together with the H.E.S.S. components is not straightforward.

**Key words.** acceleration of particles – radiation mechanisms: non-thermal – pulsars: individual: PSR J1809-1917 – gamma rays: general

## 1. Introduction

The High Energy Stereoscopic System (H.E.S.S.) Galactic Plane Survey (HGPS; H.E.S.S. Collaboration 2018a) has revealed a large number of Galactic γ-ray sources in the very-high-energy (VHE; $E > 100$ GeV) domain. While a number of these sources could be firmly associated with multi-wavelength counterparts, a large fraction of the sources remain without firm association. Among the firmly identified sources, the vast majority are either pulsar wind nebulae (PWNe) or supernova remnants (SNRs), or composite systems. In this paper, we study the unassociated

source HESS J1809-193, which has been discovered using early HGPS observations (Aharonian et al. 2007).

Identifying the physical counterpart of HESS J1809-193 is particularly challenging due to the presence of several plausible associations in its vicinity. For instance, the region harbours two energetic pulsars: PSR J1811-1925 ($E = 6.4 \times 10^{36}$ erg s$^{-1}$, $d \sim 5$ kpc) and PSR J1809-1917 ($E = 1.8 \times 10^{36}$ erg s$^{-1}$, $d \sim 3.3$ kpc; Manchester et al. 2005), which is particularly noteworthy as it powers an X-ray PWN with an extension of ~3' (Kargaltsev & Pavlov 2007; Anada et al. 2010; Klingler et al. 2018, 2020). Located nearby is also a transient X-ray magnetar, XTE J1810-197 (Alford & Halpern 2016). On the other hand, HESS J1809-193 is also spatially coincident with several SNRs, most notably G011.1+00.1 and G011.0-00.0

<footnote>* Corresponding authors;
e-mail: contact.hess@hess-experiment.eu</footnote>
(Green 2019), as well as with molecular clouds (Castelletti et al. 2016; Voisin et al. 2019). This leaves open the possibility to interpret the γ-ray emission as originating from high-energy electrons<fnref n="1" /> – most likely provided by one of the pulsars – that up-scatter photons from ambient radiation fields to γ-ray energies via the inverse Compton (IC) process ('leptonic scenario'), or as being due to interactions of high-energy cosmic-ray nuclei – accelerated, for example, at SNR shock fronts – within nearby molecular clouds ('hadronic scenario'). In the HESS J1809-193 discovery paper (Aharonian et al. 2007), as well as in two follow-up studies presented shortly afterwards (Komin et al. 2007; Renaud et al. 2008), the authors showed that the PWN surrounding PSR J1809-1917 can naturally explain the observed γ-ray emission in a leptonic scenario, and thus it represents the most likely association. However, Castelletti et al. (2016) and Araya (2018) have subsequently put forward an interpretation in a hadronic scenario involving the SNRs and molecular clouds found within the region.

Recently, the detection of γ-ray emission from HESS J1809-193 with the High Altitude Water Cherenkov Observatory (HAWC) above 56 TeV (and quite possibly above 100 TeV; Abeysekara et al. 2020; see also Goodman 2022) has added to the motivation to identify its origin. In a hadronic scenario, this detection would make HESS J1809-193 a good 'PeVatron' candidate – that is, a source capable of accelerating cosmic-ray nuclei to PeV energies. The identification of such PeVatrons is regarded as decisive in the quest for unveiling the origin of Galactic cosmic rays (e.g. Berezinskii et al. 1990; Aharonian et al. 2013; Cristofari 2021). On the other hand, should the leptonic interpretation hold, the detection by HAWC would demonstrate that HESS J1809-193 is a fascinating laboratory for the study of high-energy electrons and their propagation – and render it another extremely high-energy γ-ray source associated to a pulsar (see e.g. Sudoh et al. 2021; Albert et al. 2021a; de Oña Wilhelmi et al. 2022). In this context, we also note the recent discovery of extended halos around several energetic pulsars (e.g. Abeysekara et al. 2017a). While the term 'halo' has frequently been adopted in the literature for many γ-ray sources associated with pulsars (e.g. Linden et al. 2017), we follow here the stricter definition by Giacinti et al. (2020), who have defined it as a region where the pulsar no longer dominates the dynamics of the interstellar medium (ISM), yet where an over-density of relativistic electrons is present. The escape of the electrons from the PWN to the extended halo could, for example, be caused by an interaction of the reverse SNR shock with the pulsar wind (Blondin et al. 2001; Hinton et al. 2011).

We present here an updated study of HESS J1809-193 with H.E.S.S., based on a larger data set compared to previous publications, and employing improved analysis methods. To be able to interpret the results in a consistent manner, we complement this with a new analysis of data above 1 GeV from the Fermi-LAT space telescope for the same region. In doing so, we are able to gain new insights into the nature of HESS J1809-193.

In Sect. 2, we introduce the H.E.S.S. and Fermi-LAT data sets and analyses. The results of the analyses are presented in Sect. 3, followed by an interpretation in the framework of leptonic and hadronic models in Sect. 4. Finally, we conclude the paper in Sect. 5.

# 2. Data analysis

## 2.1. H.E.S.S. data analysis

H.E.S.S. is an array of five imaging atmospheric Cherenkov telescopes (IACTs), which detect the Cherenkov light produced in atmospheric air showers that are initiated by primary γ rays. It is situated on the Southern Hemisphere, in the Khomas highlands of Namibia (23°16′18″S, 16°30′00″E), at an altitude of 1800 m above sea level. The original array, referred to as 'HESS-I', was installed in 2000-2003 and comprised four telescopes with 12 m-diameter mirrors (CT1-4), arranged in a square layout with 120 m side length (Aharonian et al. 2006). The array was completed in 2012 by a fifth telescope, CT5, featuring a 28 m-diameter mirror, and placed in the centre of the array (Holler et al. 2015). With the full array, H.E.S.S. is sensitive to γ rays in the energy range between ~0.1 TeV and ~100 TeV.

### 2.1.1. Data set and low-level data analysis

H.E.S.S. observations on HESS J1809-193 have been carried out between 2004 and 2010, that is, exclusively during the HESS-I phase. The analysis presented here is therefore restricted to the CT1-4 telescopes. The observations are divided into 'runs' of typically 28 min duration. Selecting runs that encompass HESS J1809-193 within ~2.2° of the pointing position of the telescopes, and applying standard selection criteria for spectral studies (Aharonian et al. 2006), we obtained a data set comprising 201 runs, amounting to a total observation time of 93.2 h. This represents a significant increase with respect to the previous dedicated publications on HESS J1809-193, which used 25 h (Aharonian et al. 2007), 32 h (Komin et al. 2007), and 41 h (Renaud et al. 2008) of data.

In the data analysis, we have selected γ ray-like events using a machine learning-based method (Ohm et al. 2009), and have reconstructed their energy and arrival direction employing a maximum-likelihood fit, in which the recorded telescope images are compared to a library of simulated image templates (Parsons & Hinton 2014). We have repeated the entire analysis with an independent second analysis chain (Becherini et al. 2011), which employs different algorithms for the image calibration, event selection, and event reconstruction, obtaining compatible results. For the subsequent high-level analysis, we converted our data to the open 'GADF' format (Deil et al. 2018), and used the open-source analysis package GAMMAPY (Deil et al. 2017, 2020; v0.17).

Atmospheric air showers initiated by charged cosmic rays outnumber those resulting from γ rays by several orders of magnitude and they cannot be rejected completely in the event selection without a too severe loss of γ-ray efficiency. The modelling of the residual background due to these events (referred to as 'hadronic background' hereafter) represents one of the major challenges in any analysis of IACT data. Most established techniques rely on an estimation of the background from source-free regions in the observed field of view of the run itself (see Berge et al. 2007, for a review). We have chosen here an alternative approach, in which the residual hadronic background is provided by a background model, which we have constructed from archival H.E.S.S. observations, as detailed in Mohrmann et al. (2019). Together with the usage of GAMMAPY, this enabled us to carry out a 3-dimensional likelihood analysis of the data, that is, to model simultaneously the energy spectrum and spatial morphology of HESS J1809-193.

<footnote n="1">¹ In this paper, we use the term 'electrons' to refer to both electrons and positrons.</footnote>
The application of this analysis method to H.E.S.S. data has been validated by Mohrmann et al. (2019).

Owing to varying observation conditions (in particular the pointing zenith angle and the atmospheric transparency), a dedicated energy threshold needs to be computed for each observation run. We determined the thresholds requiring that the average bias of the energy reconstruction does not exceed 10%, and that the background model is not used below the energy at which the predicted background rate peaks (for more details, see Mohrmann et al. 2019). The resulting energy thresholds are below 0.9 TeV for all observations, while the lowest threshold, obtained for ~10% of the observations, is 0.27 TeV. Because the performance of the system degrades at large offset angles, we furthermore imposed a maximum angle between the direction of reconstructed events and the telescope pointing position of 2.2°. The value has been chosen such that the emission region is fully enclosed for all selected observations runs, many of which have been taken as part of the HGPS, implying that a considerable fraction exhibit relatively large (i.e. >1°) offset angles with respect to the centre of the emission region.

### 2.1.2. Likelihood analysis

In the likelihood analysis, the best-fit models have been obtained by minimising the quantity $-2 \log(\mathcal{L})$, where $\mathcal{L} = \prod_i P(n_i | \mu_i)$, and $P(n_i | \mu_i)$ is the Poisson probability of observing $n_i$ events in bin $i$, given a predicted number of events $\mu_i$ from the background model, and $\gamma$-ray source models if present (Mattox et al. 1996). To compute the number of events predicted by source models, we folded the source spatial model and energy spectrum with the instrument response functions (IRFs; effective area, point spread function, and energy dispersion matrix), which we derived for every observation run from extensive Monte Carlo simulations<fnref n="2" /> (Bernlöhr 2008). As spatial source models we have used two-dimensional Gaussians that can be either symmetric or elongated, represented by the GaussianSpatialModel class in GAMMAPY. As spectral models, we have used a power law (PL) of the form

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

with normalisation $N_0$, spectral index $\Gamma$, and reference energy $E_0$, as well as a power law with exponential cutoff (ECPL),

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

where $E_c$ is the cutoff energy.

We carried out the likelihood fit in a region of interest (RoI) of 6° × 6°, centred on HESS J1809–193 (see Fig. A.2). For the binning of our data, we used spatial pixels of 0.02° × 0.02° size, and an energy binning of 16 bins per decade of energy. Besides HESS J1809–193, the RoI also contains the known $\gamma$-ray sources HESS J1804–216, HESS J1813–178 (H.E.S.S. Collaboration 2018a), and HESS J1808–204 (H.E.S.S. Collaboration 2018b), which we have masked in the fit using circular exclusion regions (cf. Fig. A.2).

In the first step of the analysis, we have adjusted the background model for each observation. This background model fit is described in detail in Appendix A, where we also lay out the procedure for computing significance maps. The fit result indicates that we have achieved a very good description of the hadronic background after the adjustment.

For the further analysis, we have combined the observations into six ‘stacked’ data sets, where observations with the same energy threshold have been grouped together. This procedure effectively combines observations with similar observing conditions; further divisions of the data would lead to too many separate data sets. The six data sets are fitted jointly in the likelihood analysis. Then, we have modelled the $\gamma$-ray emission of HESS J1809–193 by adding source components to the model prediction. For nested<fnref n="3" /> models, the preference of one model over another one can be computed from the ‘test statistic’, $TS = -2 \log(\mathcal{L}_0 / \mathcal{L}_1)$, which – in the limit of sufficient statistics, and far enough from parameter boundaries – follows a $\chi^2$ distribution with $k$ degrees of freedom, where $k$ is the difference in the number of model parameters between the two tested models (Wilks 1938).

After the model fit, it is possible to extract flux points for each fitted source component. To do so, we re-ran the fit in narrow energy ranges, keeping all source model parameters except for the flux normalisation $\phi_0$ fixed to their best-fit values. The best-fit normalisation found in each energy range can then be taken as the measured flux in that range, and is quoted at its centre energy (in logarithmic space).

### 2.1.3. Estimation of systematic uncertainties

Despite the fact that we have computed customised IRFs for each observation run, due to necessary simplifying assumptions in their generation, these IRFs do not always describe the instrument and data-taking conditions perfectly. Discrepancies between the assumed IRFs and the true conditions can then lead to a systematic bias in the likelihood analysis. To assess the potential impact of mis-modelled IRFs on our fit results, we have estimated systematic uncertainties for all fit parameters. Specifically, we have considered two effects that together dominate the systematic uncertainty on our results: a shift of the global energy scale, and uncertainties of the hadronic background model. A shift of the energy scale may, for example, arise from a mis-modelling of the optical efficiency of the telescopes, or from variations in the transparency for Cherenkov radiation of the atmosphere. On the other hand, the background model has been constructed from observation runs that were taken under similar, but not identical conditions, and may therefore – despite its adjustment to the analysed observations (cf. Appendix A) – not predict the background rate perfectly.

We estimated the systematic uncertainties adopting a Monte Carlo-based approach, in which we randomly varied the IRFs according to the two systematic effects mentioned above, generated random pseudo data sets based on these IRFs and the best-fit source models, and re-fitted these pseudo data sets with the original, unmodified IRFs. The obtained spread in the fitted source model parameters then reflects their combined statistical

<footnote n="2">² Technically, custom IRFs for every observation run are obtained by interpolating between IRFs generated from Monte Carlo simulations that have been carried out for a grid of observational parameters.</footnote>

<footnote n="3">³ Two models $M_0$ and $M_1$ are considered nested if the parameters of $M_0$ are a subset of those of $M_1$, and $M_1$ can be reduced to $M_0$ for a particular choice of values for its additional parameters.</footnote>
<figure><img src="image_2.png" /></figure>

**Fig. 1.** Significance maps with best-fit models. We used an oversampling radius of 0.07° for smoothing. Significance values were obtained using the ‘Cash’ statistic (Cash 1979; see also Appendix A). Circle markers and coloured dashed lines display the best-fit position and 1-σ extent of the Gaussian model components, respectively. (a) Pre-modelling significance map, with best-fit model components of the 1-component model and the 2-component model indicated. (b) Residual significance map for the 1-component model. (c) Residual significance map for the 2-component model. In all panels, the grey dashed line marks the Galactic plane, white dashed circles show regions excluded from the analysis, and the black triangle marker denotes the position of PSR J1809−1917.

and systematic uncertainty. The procedure is described in detail
in Appendix B, and the resulting systematic uncertainties are
presented along with the analysis results in Sect. 3.1.

We note that the two considered effects potentially do not
encompass all possible sources of systematic errors in the anal-
ysis. For a more general estimate of the systematic uncertainties
of H.E.S.S., we refer to Aharonian et al. (2006), where system-
atic uncertainties of 20% on the flux normalisation and 0.1 on
the source spectral index have been derived.

## 2.2. Fermi-LAT data analysis

Fermi-LAT is a pair conversion telescope onboard the Fermi Gamma-Ray Space Telescope and is sensitive to γ rays in the energy range from ~20 MeV to several hundred GeV (Atwood et al. 2009). Here, we analysed 12 yr and 5 months of data, taken between August 4, 2008 and December 31, 2020. We used the ‘Pass 8’ IRFs (Atwood et al. 2013) and selected events passing the P8R3_SOURCE event selection (event class 128, event type 3). Because the angular resolution of Fermi-LAT substantially worsens below 1 GeV, we restricted the analysis to events above this energy. To suppress γ rays originating from the Earth’s limb, we furthermore excluded events with zenith angles above 90°. The data were analysed using FERMITOOLS<fnref n="4" /> version 2.2.0 and FERMIPY<fnref n="5" /> version 1.1.5 (Wood et al. 2017).

The analysis was carried out in a region of interest (ROI) of 10° × 10°, centred on the nominal position of 4FGL J1810.3–1925e provided in the 4FGL-DR2 catalogue (H.E.S.S. Collaboration 2020; Ballet et al. 2020). The events and exposure maps were binned using spatial bins

of 0.1° × 0.1° size and five bins per decade in energy.
γ ray source models were then fitted to the data with a
likelihood fit as described in Sect. 2.1.2, where we used
for the sources in the ROI the models provided in the
4FGL-DR2 catalogue by default. In addition, we included
standard templates for isotropic and Galactic diffuse γ-ray
emission<fnref n="6" />.

In the likelihood analysis, we have fixed the parameters of all source models with a TS value smaller than 25 or with fewer than 700 predicted events. We left free the normalisation parameter for all sources within 1° of 4FGL J1810.3–1925e, all parameters of the isotropic and Galactic diffuse model, as well as all parameters of the source models for 4FGL J1810.3–1925e and 4FGL J1811.5–1925, which immediately overlap with HESS J1809–193. Systematic uncertainties on the best-fit flux normalisations have been obtained by scaling up and down the effective area by 3% and repeating the analysis<fnref n="7" />.

3. Results

3.1. H.E.S.S. results

We show in Fig. 1a the residual significance map after the fit of
the hadronic background model, while the map in Fig. 2a dis-
plays the deduced flux of γ rays above the threshold energy
of 0.27 TeV. Extended γ-ray emission around the position of
PSR J1809−1917 (black triangle marker) is visible, we refer to

<footnote n="6">⁶ We used the file iso_P8R3_SOURCE_V2_v1.txt for the isotropic and the file gll_iem_v07.fits for the Galactic diffuse emission, respectively; see also http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html</footnote>

<footnote n="7">⁷ We have followed the procedure outlined at https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/Aeff_Systematics.html</footnote>

<footnote n="4">⁴ https://fermi.gsfc.nasa.gov/ssc/data/analysis/software</footnote>

<footnote n="5">⁵ https://fermipy.readthedocs.io</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Flux maps for HESS J1809−193, above a γ-ray energy of 0.27 TeV. (a) View of the entire emission of HESS J1809−193. (b) Zoom-in on core region. In both panels, light blue circles show the positions of known SNRs, the black triangle marker denotes the position of PSR J1809−1917, and the grey dashed line marks the Galactic plane. The position and extent of component A and component B of the two-component model are displayed in green and purple, respectively. The multi-wavelength data in panel (b) are from the JCMT (Castelletti et al. 2016), the Fugin survey (Umemoto et al. 2017), the Mopra telescope (Voisin et al. 2019), and the Suzaku X-ray telescope (Anada et al. 2010). The velocity interval for the Fugin data has been adopted from Castelletti et al. (2016). We computed the flux maps assuming a power law-type spectrum with index −2.2, and have employed a Gaussian kernel with 0.07° radius for smoothing.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Significance distributions after model fits. The red, unfilled histogram shows the distribution for all spatial pixels (outside the exclusion regions) of the significance map in Fig. 1b, while the blue, filled histogram corresponds to the map in Fig. 1c. The black, dashed line displays a Gaussian distribution with a mean $\mu = 0$ and width $\sigma = 1$.</figcaption></figure>

this source as HESS J1809−193. It is striking that besides the larger-scale emission with an extent of about 1°, a core of bright emission close to – but not fully coinciding with – the pulsar position is present.

First, we have attempted to model the emission using a 1-component model, which comprises a single source component, described by an elongated Gaussian spatial model and a

PL spectral model. The best-fit position and $1-\sigma$ extent of this component are shown in blue in Fig. 1a, and in Fig. 1b we show the residual significance map after subtracting the emission predicted by this model. The residual map shows significant remaining features close to the best-fit position, indicating that the 1-component model does not provide an acceptable description of the observed emission. This finding is confirmed by the corresponding distribution of pixel significance values, shown by the red line in Fig. 3, which clearly deviates from the expected distribution for a good description of the data (black, dashed line). This is because the larger-scale emission and the bright core cannot be simultaneously modelled by a single component that is described by a Gaussian, or any other reasonably basic spatial model. This finding holds even if we allow the extent of the Gaussian model to vary with energy, as detailed in Appendix C.1.

We have therefore adopted a 2-component model, which features in addition a second component that is described by a symmetric Gaussian spatial model and a PL spectral model. Figure 1a also shows the best-fit position and extent of both components of the 2-component model, while the residual significance map after the fit of this model is displayed in Fig. 1c, and the corresponding distribution of significance values is shown in Fig. 3 (blue histogram). As is clearly visible, the 2-component model provides a much better fit to the observed data than the 1-component model (statistically, it is preferred by $13.3\sigma$), and in fact no residual deviations except for those expected from statistical fluctuations can be made out. In the following, we
<table><caption>**Table 1.** Extent of H.E.S.S. component A in energy bands.</caption><thead><tr><th>E<sub>min</sub> - E<sub>max</sub><br>(TeV)</th><th>E<sub>mean</sub><br>(TeV)</th><th>σ<sub>major</sub><br>(deg)</th><th>σ<sub>minor</sub><br>(deg)</th></tr></thead><tbody><tr><td>0.27 - 0.75</td><td>0.43</td><td>0.69 ± 0.10</td><td>0.41 ± 0.10</td></tr><tr><td>0.75 - 2.1</td><td>1.2</td><td>0.62 ± 0.05</td><td>0.31 ± 0.04</td></tr><tr><td>2.1 - 5.6</td><td>3.2</td><td>0.62 ± 0.06</td><td>0.35 ± 0.05</td></tr><tr><td>>5.6</td><td>9.6</td><td>0.61 ± 0.11</td><td>0.46 ± 0.13</td></tr></tbody></table>

**Notes.** $E_{min}$ and $E_{max}$ denote the lower and upper boundary of the energy band, respectively, $E_{mean}$ the weighted mean energy. $\sigma_{major}$ and $\sigma_{minor}$ denote the 1-sigma extent of the semi-major and semi-minor axis of the elongated Gaussian spatial model for H.E.S.S. component A, respectively. See Appendix C.3 for further details.

refer to the extended component as ‘component A’, and to the compact component as ‘component B’. In Appendix C.2, we study the agreement between the 2-component model and the observed data as a function of energy, finding that the model provides a good fit at all energies. Finally, we have explored in Appendix C.3 how the parameters of component A vary when the model is fitted in separate energy bands. We find that the fitted parameters do not change significantly, and provide in Table 1 the fitted extent of component A in the four employed energy bands.

In Fig. 2b, we provide a detailed view of the inner region of HESS J1809–193, with multi-wavelength data overlaid. The peak of the emission (and, by that, the position of component B), is offset by ~7' from the position of PSR J1809–1917 and its surrounding X-ray nebula, indicated by the brown contour lines (Anada et al. 2010). On the other hand, component B lies with its centre point directly on the western edge of the SNR G011.0+00.0, and is furthermore spatially coincident with dense molecular clouds as indicated by $^{12}$CO ($J = 3-2$) observations by the James Clerk Maxwell Telescope (JCMT; Castelletti et al. 2016) and CS observations by the Mopra telescope (Voisin et al. 2019). Moreover, the contour lines from the FUGIN $^{12}$CO ($J = 1-0$) survey (Umemoto et al. 2017) illustrate that molecular gas is present throughout the region. In Appendix D, we provide a map of the FUGIN $^{12}$CO data with the two components of HESS J1809–193 overlaid.

Finally, we show in Fig. 4 the energy spectra and flux points obtained for the two components of HESS J1809–193, and compare these to previously obtained results from the literature. We note that, when comparing the fitted PL models for components A and B (solid lines), the spectrum of component B appears somewhat harder than that of components A ($\Gamma = 1.98 \pm 0.05$ vs. $\Gamma = 2.24 \pm 0.03$). However, the flux upper limits obtained above energies of ~20 TeV for component A seem to indicate the presence of a cutoff to the spectrum. We have therefore repeated the analysis adopting an ECPL spectral model for component A, which led to no changes in the best-fit parameters values for component B or for the spatial model of component A, but yielded a flatter spectrum ($\Gamma = 1.90 \pm 0.05$) also for component A at energies below ~5 TeV (dashed line in Fig. 4). With respect to the PL spectral model, the ECPL model is preferred with a statistical significance of ~8$\sigma$. For component B, on the other hand, we found no significant preference for a cutoff to the spectrum.

The spectrum of component A, which dominates the total emission except at the highest energies, is well compatible with

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

**Fig. 4.** H.E.S.S. energy spectrum results. We show in green and purple the flux points for component A and component B of the H.E.S.S. analysis, respectively. Upper limits are at 95% confidence level. The solid lines with shaded bands display the best-fit PL model and statistical uncertainty for each of the components. For component A, the dashed green line shows the best-fit ECPL model in addition. The energy spectra are compared to published results, taken from Aharonian et al. (2007), H.E.S.S. Collaboration (2018a), Zarić et al. (2021), Abeysekara et al. (2020), and Goodman (2022).

that published in H.E.S.S. Collaboration (2018a), and may also be reconcilable with the high-energy emission measured with HAWC (Abeysekara et al. 2020; Goodman 2022), although a more gradual decrease of the flux than predicted by the ECPL spectral model would be required in this case. That the flux points from Aharonian et al. (2007) indicate a lower flux compared to the one found here can be understood when considering that the flux was extracted only from within a circular region of 0.5° radius, and thus part of the larger-scale emission has been missed.

We summarise the best-fit parameter values found for the 2-component model in Table 2, along with the corresponding statistical and systematic uncertainties (the latter having been derived as described in Appendix B). For component A, we provide the results both for the PL spectral model and the ECPL spectral model.

## 3.2. Fermi-LAT results

The Fermi-LAT 4FGL-DR2 catalogue (H.E.S.S. Collaboration 2020; Ballet et al. 2020) lists two sources that are located in the immediate vicinity of HESS J1809–193$^8$: (i) 4FGL J1810.3–1925e is modelled as an extended source (using a two-dimensional Gaussian as spatial model) and its spectrum is fitted with a log-parabola model,

$$ \frac{dN}{dE} = N_0 \cdot \left(\frac{E}{E_0}\right)^{-\alpha-\beta \log(E/E_0)} \quad (3) $$

(with log the natural logarithm); (ii) 4FGL J1811.5–1925 is modelled as a point-like source and its spectrum is fitted with

<footnote>$^8$ We have used the 4FGL-DR2 catalogue as a basis for our analysis, but have checked that the region is modelled in the same way in the more recent 4FGL-DR3 catalogue (H.E.S.S. Collaboration 2022).</footnote>
<table><caption>**Table 2.** Best-fit parameter values for the H.E.S.S. 2-component model.</caption><thead><tr><th>Par. [unit]</th><th>Value</th></tr></thead><tbody><tr><td colspan="2">Component A (PL spectral model)</td></tr><tr><td>RA [deg]</td><td>272.551 ± 0.025<sub>stat</sub> ± 0.018<sub>sys</sub><br>(18<sup>h</sup>10<sup>m</sup>12<sup>s</sup> ± 6<sup>s</sup><sub>stat</sub> ± 4<sup>s</sup><sub>sys</sub>)</td></tr><tr><td>Dec. [deg]</td><td>-19.344 ± 0.023<sub>stat</sub> ± 0.013<sub>sys</sub><br>(-19°20.6' ± 1.4'<sub>stat</sub> ± 0.8'<sub>sys</sub>)</td></tr><tr><td>σ [deg]</td><td>0.622 ± 0.032<sub>stat</sub> ± 0.020<sub>sys</sub><sup>(a)</sup></td></tr><tr><td>e</td><td>0.824 ± 0.025<sub>stat</sub></td></tr><tr><td>φ [deg]</td><td>50.0 ± 3.1<sub>stat</sub></td></tr><tr><td>N<sub>0</sub> [10<sup>-12</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>]</td><td>8.42 ± 0.40<sub>stat</sub> ± 1.14<sub>sys</sub></td></tr><tr><td>Γ</td><td>2.239 ± 0.027<sub>stat</sub> ± 0.020<sub>sys</sub></td></tr><tr><td>E<sub>0</sub> [TeV]</td><td>1 (fixed)</td></tr><tr><td colspan="2">Component A (ECPL spectral model)</td></tr><tr><td>RA [deg]</td><td>272.554 ± 0.025<sub>stat</sub> ± 0.019<sub>sys</sub><br>(18<sup>h</sup>10<sup>m</sup>13<sup>s</sup> ± 6<sup>s</sup><sub>stat</sub> ± 5<sup>s</sup><sub>sys</sub>)</td></tr><tr><td>Dec. [deg]</td><td>-19.344 ± 0.021<sub>stat</sub> ± 0.012<sub>sys</sub><br>(-19°20.6' ± 1.3'<sub>stat</sub> ± 0.7'<sub>sys</sub>)</td></tr><tr><td>σ [deg]</td><td>0.613 ± 0.031<sub>stat</sub> ± 0.015<sub>sys</sub><sup>(a)</sup></td></tr><tr><td>e</td><td>0.820 ± 0.025<sub>stat</sub></td></tr><tr><td>φ [deg]</td><td>51.3 ± 3.1<sub>stat</sub></td></tr><tr><td>N<sub>0</sub> [10<sup>-12</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>]</td><td>9.05 ± 0.47<sub>stat</sub> ± 0.91<sub>sys</sub></td></tr><tr><td>Γ</td><td>1.90 ± 0.05<sub>stat</sub> ± 0.05<sub>sys</sub></td></tr><tr><td>E<sub>c</sub> [TeV]</td><td>12.7<sup>+2.7</sup><sub>-2.1</sub><sup>stat</sup> ± 1.2<sup>+2.6</sup><sub>-1.9</sub><sup>stat</sup></td></tr><tr><td>E<sub>0</sub> [TeV]</td><td>1 (fixed)</td></tr><tr><td colspan="2">Component B</td></tr><tr><td>RA [deg]</td><td>272.400 ± 0.010<sub>stat</sub><br>(18<sup>h</sup>09<sup>m</sup>36<sup>s</sup> ± 2.4<sup>s</sup><sub>stat</sub>)</td></tr><tr><td>Dec. [deg]</td><td>-19.406 ± 0.009<sub>stat</sub><br>(-19°24.4' ± 0.5'<sub>stat</sub>)</td></tr><tr><td>σ [deg]</td><td>0.0953 ± 0.0072<sub>stat</sub> ± 0.0034<sub>sys</sub></td></tr><tr><td>N<sub>0</sub> [10<sup>-12</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>]</td><td>0.95 ± 0.11<sub>stat</sub> ± 0.011<sub>sys</sub></td></tr><tr><td>Γ</td><td>1.98 ± 0.05<sub>stat</sub> ± 0.03<sub>sys</sub></td></tr><tr><td>E<sub>0</sub> [TeV]</td><td>1 (fixed)</td></tr></tbody></table>

**Notes.** For component A, we provide the best-fit values for the assumption of a PL spectral model and an ECPL spectral model. *σ*, *e*, and *φ* denote the 1-*σ* radius, eccentricity, and position angle of the Gaussian spatial model, respectively. *N*<sub>0</sub>, *Γ*, *E*<sub>0</sub>, and *E*<sub>c</sub> are parameters of the spectral models, as defined in Eqs. (1) and (2). Systematic uncertainties have been derived as described in Appendix B. <sup>(a)</sup>For the asymmetric component A, *σ* refers to the semi-major axis. The extents of the semi-minor axis compute to 0.353 ± 0.029<sub>stat</sub> and 0.351 ± 0.028<sub>stat</sub> for the PL and ECPL spectral model, respectively.

a power-law model (cf. Eq. (1)). We confirm with our analysis the presence of both sources, that is, we were not able to obtain a satisfactory fit with only one source, or a different choice of spatial models. In the following, we refer to the source models we obtained as J1810.3–1925e and J1811.5–1925, respectively, in distinction to the models provided in the 4FGL-DR2 catalogue. The best-fit parameter values of the models are summarised in Table 3. The systematic uncertainties on all model parameters except the flux normalisation *N*<sub>0</sub> are negligible compared to the statistical ones, and thus not quoted in the table. We acknowledge that our best-fit spectral model for J1810.3–1925e shows no significant curvature, but have decided to maintain the log-parabola model for consistency with the 4FGL-DR2 catalogue (which extended to lower energies than our analysis here).

<table><caption>**Table 3.** Best-fit parameter values for the Fermi-LAT data analysis.</caption><thead><tr><th>Par. [unit]</th><th>Value</th></tr></thead><tbody><tr><td colspan="2">J1810.3–1925e</td></tr><tr><td>RA [deg]</td><td>272.547 ± 0.033<sub>stat</sub><br>(18<sup>h</sup>10<sup>m</sup>11<sup>s</sup> ± 8<sup>s</sup><sub>stat</sub>)</td></tr><tr><td>Dec. [deg]</td><td>-19.397 ± 0.038<sub>stat</sub><br>(-19°23.8' ± 2.3'<sub>stat</sub>)</td></tr><tr><td>σ [deg]</td><td>0.317 ± 0.024<sub>stat</sub></td></tr><tr><td>N<sub>0</sub> [10<sup>-9</sup> GeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>]</td><td>2.86 ± 0.27<sub>stat</sub> ± 0.09<sub>sys</sub></td></tr><tr><td>α</td><td>2.53 ± 0.11<sub>stat</sub></td></tr><tr><td>β</td><td>-0.015 ± 0.043<sub>stat</sub></td></tr><tr><td>E<sub>0</sub> [GeV]</td><td>1.747 (fixed)</td></tr><tr><td colspan="2">J1811.5–1925</td></tr><tr><td>RA [deg]</td><td>272.874 ± 0.025<sub>stat</sub><br>(18<sup>h</sup>11<sup>m</sup>30<sup>s</sup> ± 6<sup>s</sup><sub>stat</sub>)</td></tr><tr><td>Dec. [deg]</td><td>-19.410 ± 0.027<sub>stat</sub><br>(-19°24.6' ± 1.6'<sub>stat</sub>)</td></tr><tr><td>N<sub>0</sub> [10<sup>-9</sup> GeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>]</td><td>0.016 ± 0.003<sub>stat</sub> ± 0.001<sub>sys</sub></td></tr><tr><td>Γ</td><td>2.42 ± 0.16<sub>stat</sub></td></tr><tr><td>E<sub>0</sub> [GeV]</td><td>7.48 (fixed)</td></tr></tbody></table>

**Notes.** *σ* is the 1-*σ* radius of the Gaussian spatial model. *N*<sub>0</sub>, *α*, *β*, *Γ*, and *E*<sub>0</sub> are parameters of the spectral models defined in Eqs. (3) and (1).

Removing J1810.3–1925e and J1811.5–1925 from the best-fit ROI model, we obtained the significance map shown in Fig. 5, panel (a). Panel (b) displays the significance map after adding J1811.5–1925 to the ROI model, while panel (c) shows the map with both source models restored. A comparison of the panels shows that the two fitted sources account for the majority of the emission around HESS J1809–193. The energy spectra extracted for J1810.3–1925e and J1811.5–1925 are displayed in Fig. 6. The obtained spectrum for J1810.3–1925e is in good agreement with that obtained by Araya (2018), considering that the region has been modelled differently (a single source with a disk spatial model) there.

The point-like source J1811.5–1925 is positionally coincident with PSR J1811–1925, which strongly suggests an association with this pulsar, as also listed in the 4FGL-DR2 catalogue. We therefore regard its emission as unrelated to HESS J1809–193. On the other hand, the best-fit position of J1810.3–1925e is close to PSR J1809–1917 and the two H.E.S.S. source components, suggesting a connection to HESS J1809–193. In particular, the fitted position and extent are very similar to those of the extended H.E.S.S. component (component A), as is evident from Fig. 5. In order to further explore the connection between the Fermi-LAT and H.E.S.S. data, we have also extracted energy spectra of the emission observed with Fermi-LAT using the best-fit spatial models of the two H.E.S.S. components as spatial templates (removing J1810.3–1925e from the model but retaining J1811.5–1925). The result is shown in Fig. 7. As expected due to its slightly larger spatial extent, the spectrum obtained for the template of component A is slightly above that of J1810.3–1925e. With the template of component B we obtained only flux upper limits, however, this is not a surprise given Fermi-LAT’s broadband sensitivity<fnref n="9" /> (dashed line in Fig. 7).

<footnote n="9">⁹ Broadband sensitivity curves for Fermi-LAT are available at https://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm. We have used the curve for Galactic coordinates *l* = 0° and *b* = 0°.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Significance maps above 1 GeV from the Fermi-LAT analysis. (a) With J1811.5–1925 and J1810.3–1925e removed from the best-fit model. (b) With J1810.3–1925e removed from the best-fit model. (c) With all sources included in the model. The light blue cross denotes the fitted position of J1811.5–1925, whereas the dark blue cross and dashed circle display the fitted position and $1-\sigma$ extent of J1810.3–1925e. For comparison, the components of the 2-component model fitted to the H.E.S.S. data are shown in green and purple (same as in Fig. 1). The grey dashed line marks the Galactic plane, while the coloured triangle markers denote the positions of PSR J1809–1917, PSR J1811–1925, and XTE J1810–197.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Fermi-LAT energy spectrum results for J1810.3–1925e and J1811.5–1925. We compare our spectra to that obtained by Araya (2018) for the entire region.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Fermi-LAT energy spectra obtained with H.E.S.S. model templates. The dashed grey line shows the 10-yr Fermi-LAT broadband sensitivity.</figcaption></figure>

# 4. Modelling

In this section, we present an interpretation of the observational results by means of modelling the primary cosmic-ray particle populations responsible for the observed γ-ray emission. We investigate two scenarios: (i) that the emission detected with H.E.S.S. is entirely attributed to the PWN of PSR J1809–1917, that is, of purely leptonic origin (Sect. 4.1); (ii) that there is an additional contribution to the emission from hadronic cosmic rays accelerated in one of the SNRs and interacting in the molecular clouds (Sect. 4.2).

## 4.1. Pulsar wind nebula scenario

We employed a one-zone PWN model, in which we performed a time-dependent modelling of the pulsar energy output, the

ambient magnetic field, and the injected electrons, following the approach outlined in Albert et al. (2021b). The parameters of the model are summarised in Table 4. The input parameters consist mostly of measured properties of PSR J1809–1917. For the pulsar braking index, however, we have adopted the canonical value of $n = 3$, assuming that this is more representative for the full history of the pulsar than the recent measurement of $n = 23.5$, which may be affected by undetected glitches of the pulsar (Parthasarathy et al. 2019, 2020). The electron injection spectrum follows a power law with index $-\alpha$ and an exponential cut-off at energy $E_c$. Its normalisation is proportional to $\theta \times \dot{E}$, that is, coupled to the time-dependent spin-down power. As specified in Gaensler & Slane (2006) and Venter & de Jager (2007), we took the time evolution of the pulsar period as $P(t) = P_0(1 + t/\tau_0)^{0.5}$, of the pulsar spin-down power as $\dot{E}(t) = \dot{E}_0(1 + t/\tau_0)^{-2}$, and of the magnetic field as $B(t) = B_0[1 + (t/\tau_0)^{0.5}]^{-1}$, where
<table><caption>**Table 4.** Parameters of the PWN model.</caption><thead><tr><th>Par.</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td colspan="3">Input parameters</td></tr><tr><td><i>d</i></td><td>Pulsar distance <sup>(a)</sup></td><td>3.3 kpc</td></tr><tr><td><i>E&#x0302;</i></td><td>Pulsar spin-down power <sup>(a)</sup></td><td>1.8 × 10<sup>36</sup> erg s<sup>-1</sup></td></tr><tr><td><i>τ</i><sub>c</sub></td><td>Pulsar characteristic age <sup>(a)</sup></td><td>51.4 kyr</td></tr><tr><td><i>P</i></td><td>Pulsar period <sup>(a)</sup></td><td>82.76 ms</td></tr><tr><td><i>P&#x0302;</i></td><td>Pulsar period derivative <sup>(a)</sup></td><td>2.55 × 10<sup>-14</sup> s s<sup>-1</sup></td></tr><tr><td><i>n</i></td><td>Pulsar braking index <sup>(b)</sup></td><td>3</td></tr><tr><td colspan="3">Adjusted parameters</td></tr><tr><td><i>θ</i></td><td>Electron power fraction</td><td>0.6</td></tr><tr><td><i>B</i></td><td>Magnetic field</td><td>4 µG</td></tr><tr><td><i>P</i><sub>0</sub></td><td>Pulsar birth period</td><td>50 ms</td></tr><tr><td><i>E</i><sub>c</sub></td><td>Cutoff energy</td><td>420 TeV</td></tr><tr><td><i>α</i></td><td>Injection spectrum index</td><td>2.0</td></tr><tr><td><i>τ</i><sub>young</sub></td><td>Age of young e<sup>-</sup></td><td>1.2 kyr</td></tr><tr><td><i>τ</i><sub>med</sub></td><td>Age of medium-age e<sup>-</sup></td><td>4.7 kyr</td></tr></tbody></table>

**Notes.** Pulsar parameters denote present-day values unless otherwise specified. The parameter values for the 'adjusted' parameters were obtained using an MCMC method, but should be regarded as indicative values rather than precise fit results (see main text). <sup>(a)</sup>Taken from Manchester et al. (2005). <sup>(b)</sup>Assumed value.

$\tau_0 = P_0/(2\dot{P}_0)$ is the initial spin-down time scale. We then computed the non-thermal emission from the injected electrons (i.e. synchrotron radiation and IC emission) employing the GAMERA library (Hahn 2015), which takes into account cooling losses of the electrons. For the IC target photon fields, we have used the model by Popescu et al. (2017)<fnref n="10" />. Finally, we have fitted the adjustable parameters of the model to the observed spectral energy distribution (SED) of HESS J1809–193, where the optimisation has been carried out using a Markov chain Monte Carlo (MCMC) method implemented in the EMCEE package (Foreman-Mackey et al. 2013). We note that some of the model parameters are correlated with each other or not well constrained by the available data. Therefore, we stress that while we have carried out an optimisation of the model parameters, the obtained values should not be regarded as measurements of the corresponding quantities, but rather as one possible combination of parameter values that yield a reasonable description of the observational data. The parameter values given in Table 4 are those that yielded the highest numerical probability, that is, the best fit to the data.

In the model, we invoked three 'generations' of electrons: (i) 'relic' electrons, which have been injected over the lifetime of the system ($\tau \approx 33$ kyr<fnref n="11" />) and are associated with the extended H.E.S.S. component (A); (ii) 'medium-age' electrons, which have been injected within the last $\tau_{med} \approx 4.7$ kyr and are associated with the compact H.E.S.S. component (B); (iii) 'young' electrons, which have been injected within the last $\tau_{young} \approx 1.2$ kyr and are associated with the X-ray nebula. In this picture, the relic electrons

are assumed to have escaped from the central region (which contains the X-ray PWN and the compact component B) at some instant in the past. For lack of evidence when this escape has occurred, the relic electrons are injected from the birth of the system until the medium-age electrons start to be injected. We note that, despite associating the different generations with different spatial regions, we have not performed a spatial modelling – the association is made in terms of the SED only.

In addition to the already presented H.E.S.S. spectra, we used in the fit the spectrum of the X-ray nebula as measured by Anada et al. (2010) with the *Suzaku* satellite between 2 and 10 keV. Since we associate only the most recently injected young electrons with the X-ray nebula, we integrated the measured flux in the immediate vicinity of the nebula only (regions 2, 3, 6, and 7 in Fig. 4 and Table 4 of Anada et al. 2010). Additionally, we derived an upper limit (at 95% confidence level) for the X-ray flux emitted by the medium-age electrons using the measured flux in regions 9–16 and applying a scaling factor that takes into account the difference in solid angle between these regions and the compact H.E.S.S. component associated with the medium-age electrons. The upper limit is not used in the fit and only serves as a sanity check for the model.

We show the obtained SEDs for the three generations of electrons in Fig. 8, together with the observed data. The model describes the spectra measured with H.E.S.S. and *Suzaku* well, and the predicted X-ray flux of the medium-age electrons does not exceed the *Suzaku* upper limit. The fit yields, for example, a moderate required present-day magnetic field of ~4 µG and a reasonable spectral index for the injection spectrum of ~2. Furthermore, a maximum electron energy of several hundred TeV is implied by the data. The total predicted $\gamma$-ray spectrum is also well compatible with the total flux from HESS J1809–193 as measured by HAWC (Goodman 2022).

The model fails, however, to explain the spectrum of the Fermi-LAT source J1810.3–1925e below ~10 GeV. This would require an additional IC component, emitted by electrons even older than the relic electrons. In this case, however, it would be expected that the emission of J1810.3–1925e exhibits a larger spatial extent than that of component A of HESS J1809–193, which is not the case. Alternatively, a hadronic component related to the SNR G011.0–00.0 could be invoked – this scenario will be discussed in more detail in Sect. 4.2.

The offset between component B and PSR J1809–1917 may be explained, for example, by the proper motion of the pulsar. Indeed, Klingler et al. (2018, 2020) have detected a northward proper motion of ~20–40 mas yr⁻¹, albeit not with high significance. This would imply a travel time between the best-fit position of component B and the current pulsar position of ~10–20 kyr. This is somewhat larger than our estimate of the age of the medium-age electrons associated with component B. However, considering that also an asymmetric crushing of the PWN by the SNR reverse shock can lead to a displacement between the PWN and the pulsar (Blondin et al. 2001), the scenario still appears feasible.

Having derived the expected age of the PWN system, we used our measurement of the size of the extended H.E.S.S. component to infer how fast the relic electrons associated to this component have diffused since their injection (see Fig. 9). We have assumed an energy-dependent diffusion coefficient

$$D = D_0 \left( \frac{E_e}{40 \text{ TeV}} \right)^\delta, \quad (4)$$

where $E_e$ is the electron energy, $D_0$ denotes the diffusion coefficient at a reference energy of 40 TeV, and $\delta$ specifies the energy

<footnote n="10">¹⁰ We note that the prediction of this large-scale model may not be very accurate in the specific region studied here. As our conclusions are based on order-of-magnitude estimates, however, they are unaltered even if the predicted radiation field densities are wrong by a factor of a few.</footnote>

<footnote n="11">¹¹ The 'true' age of the pulsar can be computed as $\tau = (P/((n-1)\dot{P})) \cdot (1 - (P_0/P)^{n-1})$ (Gaensler & Slane 2006). We note that this formula depends on the unknown pulsar birth period $P_0$, for which we have used the value suggested by our model optimisation (cf. Table 4). Other values of $P_0$ will lead to different estimates of the pulsar age.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** SED of HESS J1809−193 for the leptonic model. (a) Full energy range. (b) Zoom into high-energy regime. SED curves for the three assumed electron generations, obtained with GAMERA, are shown with dark grey lines. In panel (b), the light grey lines show individual solutions from the MCMC sampling, and thus give an indication of the statistical spread. The H.E.S.S. and Fermi-LAT data points have been derived in this work. The Suzaku X-ray data are from Anada et al. (2010), where the butterfly corresponds to the young electrons (dotted line) and the upper limit refers to the medium-age electrons (dashed line; see main text for details). Shown for comparison but not used in the fit are radio data for G011.0−00.0 (Brogan et al. 2006) and data points from HAWC (Goodman 2022).</figcaption></figure>

dependence of the diffusion. Using again the GAMERA library to derive the expected size of the relic electron component as a function of γ-ray energy, we determined the two parameters $D_0$ and $\delta$ by fitting the expected size to the observed size of component A of HESS J1809−193 (cf. Table 1) – noting again that the results of the fit are strongly model-dependent and should not be taken as a measurement. The best-fit diffusion coefficient of $D_0 \sim 1.1 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}$ appears reasonable and is of the same order of magnitude as the coefficient obtained for the Geminga halo by Abeysekara et al. (2017a). On the other hand, the observed data do not provide very strong constraints for $\delta$, with both Kolgoromov scaling ($\delta = 1/3$) and Bohm scaling ($\delta = 1$) consistent with the observations. While our simple estimate assumes a radially symmetric diffusion of the electrons, we note that the elongation of component A aligns with the asymmetric extension of the X-ray PWN, possibly hinting at a particular arrangement of the magnetic field in the region. Lastly, we point out that because the highest-energy relic electrons have cooled since they were injected, a cutoff to the corresponding γ-ray spectrum is expected to occur. The measured cutoff energy of ~13 TeV for component A of HESS J1809−193 is well in line with this prediction, as can be seen in Fig. 8b.

We conclude that the appearance of HESS J1809−193 is compatible with that of a halo of old electrons (component A) around the PWN (component B & X-ray emission). We also note that in terms of its X-ray-to-TeV luminosity ratio, PSR J1809−1917 fits well into the population of PWN (Kargaltsev et al. 2013).

## 4.2. Possible hadronic contributions

Given the presence of SNRs and molecular clouds in the vicinity of HESS J1809−193, we also need to consider the possibility that cosmic-ray nuclei accelerated at the SNR shock fronts and interacting hadronically in the molecular clouds are responsible for at least part of the observed γ-ray emission. Indeed, a mixed (leptonic+hadronic) scenario seems possible in principle: while we are not aware of distance estimates for

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Measured and predicted radius of the extended component of HESS J1809−193 (component A). As measured radius, we used the $1-\sigma$ extent of the semi-major axis of the elongated Gaussian spatial model for component A (cf. Table 1). The solid blue curve has been obtained by fitting the electron diffusion parameters ($D_0, \delta$; cf. Eq. (4)) to the measured data points; the dashed orange and dotted green curves show results for fixed values of $\delta$ as indicated in the legend.</figcaption></figure>

$G011.1+00.1$, existing distance estimates for $G011.0-00.0$ of $2.6 \text{ kpc}$ (Bamba et al. 2003), ~$3 \text{ kpc}$ (Castelletti et al. 2016), and $2.4 \pm 0.7 \text{ kpc}$ (Shan et al. 2018) seem broadly consistent with those for PSR J1809−1917 of $3.7 \text{ kpc}$ (Morris et al. 2002) and $3.27 \text{ kpc}$ (Parthasarathy et al. 2019). Furthermore, molecular gas is present throughout the region (cf. Appendix D), and in particular the dense molecular clouds found by Castelletti et al. (2016) and Voisin et al. (2019) seem to lie at distances compatible with that of $G011.0-00.0$, thus providing the required target material for cosmic-ray interactions. This has led Voisin et al. (2019) to propose that $G011.0-00.0$ is the host SNR of PSR J1809−1917. However, while the pulsar proper motion could be compatible with this scenario, the association is not firm (Klingler et al. 2018, 2020).
<table><caption>**Table 5.** Best-fit parameter values for the hadronic *pp* models.</caption><thead><tr><th>Par. [unit]</th><th>Value</th></tr></thead><tbody><tr><td colspan="2">Component A</td></tr><tr><td>N<sub>0</sub><sup>p,A</sup> [10<sup>34</sup> eV<sup>-1</sup>]</td><td>10.0 ± 1.8</td></tr><tr><td>Γ<sup>p,A</sup></td><td>1.76 ± 0.17</td></tr><tr><td>E<sub>c</sub><sup>p,A</sup> [TeV]</td><td>90<sup>+35</sup><sub>-30</sub></td></tr><tr><td>E<sub>0</sub><sup>p,A</sup> [TeV]</td><td>20 (fixed)</td></tr><tr><td colspan="2">Component B</td></tr><tr><td>N<sub>0</sub><sup>p,B</sup> [10<sup>34</sup> eV<sup>-1</sup>]</td><td>1.0 ± 0.4</td></tr><tr><td>Γ<sup>p,B</sup></td><td>1.34 ± 0.45</td></tr><tr><td>E<sub>c</sub><sup>p,B</sup> [TeV]</td><td>110<sup>+135</sup><sub>-50</sub></td></tr><tr><td>E<sub>0</sub><sup>p,B</sup> [TeV]</td><td>20 (fixed)</td></tr></tbody></table>

**Notes.** $N_0$, $\Gamma$, $E_0$, and $E_c$ are parameters of the ECPL spectral model (Eq. (2)), where the superscript *p*, *X* denotes that these are the parameters of the primary proton spectrum of component X = {A, B}, respectively. The quoted errors represent statistical uncertainties only.

Although the measured spectrum of the *Fermi*-LAT source J1810.3–1925e is comparatively soft, it could in principle be described (below ~10 GeV) using a hadronic model. However, it remains unclear why the spatial model of J1810.3–1925e coincides with that of the extended H.E.S.S. component in this case, as we would rather expect the emission to be more compact and centred on the positions of the molecular clouds. We also note that simultaneously modelling the emission of J1810.3–1925e and either of the two H.E.S.S. components in a purely hadronic scenario would require, in order to explain the transition from the steep spectrum of the former to the harder spectra of the latter, a spectral hardening in the primary cosmic-ray spectrum, for which there is no obvious explanation.

As presented in Sect. 4.1, both components of HESS J1809–193 can be modelled well within a leptonic scenario. Nevertheless, we have also explored the implications of either of the components being of hadronic origin. To this end, we have fitted a proton-proton (*pp*) model to both components, employing the NAIMA package (Zabalza 2015). The primary proton spectrum is described using an ECPL model (see Eq. (2)) and we have assumed a distance to the source of 3 kpc. We used the wrapper class for NAIMA models implemented in GAMMAPY, so that they could be fitted directly to the H.E.S.S. data sets (as opposed to fitting them to the extracted flux points only), using the same likelihood framework as before (cf. Sect. 2.1.2). The fit results are presented in Table 5, and the resulting spectra displayed in Fig. 10. The same spatial models as in the previous analysis (cf. Sect. 3.1) were assumed and compatible best-fit parameters were obtained for them.

The *pp* model for component A prefers a relatively hard spectral index of $\Gamma^{p,A} = 1.76 \pm 0.17$. Integrating the primary spectrum above 1 GeV yields a total required energy of $W^{p,A} \sim 3.2 \times 10^{50} (n/1 \text{ cm}^{-3})^{-1}$ erg, which represents – unless high ISM densities are invoked – a significant fraction of the canonically assumed kinetic energy released in a supernova explosion of ~$10^{51}$ erg (e.g. Ginzburg 1975). In this context, we note that while the FUGIN CO data indeed show the presence of molecular gas beyond the dense clouds discovered by Castelletti et al. (2016) and Voisin et al. (2019), they also indicate a gradient in the gas density across the extent of HESS J1809–193 (see the map in Appendix D). This gradient is not reflected in the

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** SED of HESS J1809–193 with hadronic (*pp*) models. The H.E.S.S. and *Fermi*-LAT data points have been derived in this work, the HAWC data points (shown for comparison only) are taken from Goodman (2022). The lines show the predicted $\gamma$-ray spectra for hadronic models fitted to each of the two H.E.S.S. components, respectively.</figcaption></figure>

observed $\gamma$-ray emission, which makes the interpretation of component A in a hadronic scenario challenging.

For component B, our fit yields an even harder proton spectral index of $\Gamma^{p,B} = 1.35 \pm 0.45$. The required energy in protons above 1 GeV is $W^{p,B} \sim 2.7 \times 10^{49} (n/1 \text{ cm}^{-3})^{-1}$ erg. Even considering that only part of the cosmic rays potentially accelerated by G011.0–00.0 will reach the dense molecular clouds, the high density of gas in the vicinity of component B in general (~1000 cm$^{-3}$, cf. Appendix D) makes this energy seem well affordable. An explanation of component B of HESS J1809–193 in a hadronic scenario therefore appears entirely reasonable, and would furthermore explain its offset from the position of PSR J1809–1917 and the X-ray PWN without the need of requiring, for example, a large proper motion velocity of the pulsar. It would appear natural in this case to associate component B also with the highest-energy emission up to 100 TeV measured with HAWC. However, as is evident from Fig. 10, the fitted cutoff energy $E_c^{p,B} \sim 110$ TeV of the primary proton spectrum leads to a too strong cutoff in the $\gamma$-ray spectrum, leaving the highest HAWC flux points unexplained. We have therefore repeated the fit of the hadronic models, adding a $\chi^2$ term that denotes the deviation between the sum of the predicted $\gamma$-ray fluxes of both H.E.S.S. components and the HAWC flux points to the total TS. In this case, we obtain slightly softer spectra ($\Gamma^{p,A} = 1.95 \pm 0.10$; $\Gamma^{p,B} = 1.56 \pm 0.22$) and slightly higher cutoff energies ($E_c^{p,A} = 140^{+80}_{-50}$ TeV; $E_c^{p,B} = 200^{+420}_{-130}$ TeV). These values are, however, consistent within uncertainties with those obtained in the previous fit, which demonstrates that it is possible to also explain the measured HAWC flux points within a hadronic scenario (or a mixed one, in which component A is leptonic and component B is hadronic).

Finally, we note that while the relatively hard primary proton spectra obtained for both components of HESS J1809–193 are not consistent with generic predictions of diffuse shock acceleration (Bell 2013), they are compatible with a scenario in which cosmic rays accelerated in a supernova remnant illuminate a nearby gas cloud (e.g. Gabici et al. 2009). On the other hand, there is also the possibility of a continuous wind of hadronic
cosmic rays powered by the pulsar (e.g. Gallant & Arons 1994; Amato et al. 2003), which may be an interesting scenario to explore for the case of PSR J1809−1917, as already noted by Voisin et al. (2019).

## 5. Conclusion

We have presented a new analysis of the γ-ray emission from HESS J1809−193, employing improved analysis techniques. For the first time, we were able to resolve the emission into two distinct components, which we have modelled with Gaussian spatial models. Component A appears extended and elongated, with a $1-\sigma$ semi-major and semi-minor axis of ~0.62° and ~0.35°, respectively, and exhibits a spectrum with a cutoff at ~13 TeV. Superimposed, component B appears symmetric and more compact with a $1-\sigma$ radius of ~0.1°, and shows a harder spectrum with no clear cutoff.

We have interpreted the results in a leptonic scenario, in which the γ-ray emission is due to high-energy electrons provided by the energetic pulsar PSR J1809−1917, which is known to power an X-ray PWN. The model is based on three generations of electrons, associated with component A, component B, and the X-ray PWN, respectively (going from old to recently injected electrons). The measured extent and spectrum of component A are compatible with a halo of old electrons that have escaped the PWN.

The presence of SNRs and molecular clouds within the region suggests that (part of) the γ-ray emission could also be of hadronic origin. Indeed, we found that both of the components of HESS J1809−193 can in principle be modelled within a hadronic scenario. However, a lack of correlation between the γ-ray emission of component A and the gas present in the region disfavours a hadronic interpretation for this component. Conversely, for component B, which is spatially coincident with the shell of the SNR G011.0−00.0 and several molecular clouds, this is a viable alternative explanation. The measurement of γ-ray emission up to 100 TeV with HAWC could be viewed as additional support for this interpretation. It would, however, leave the X-ray PWN without a counterpart at TeV energies (component A being associated with electrons injected long ago only), which would be unexpected when comparing with other PWN systems (Kargaltsev et al. 2013).

Our analysis of Fermi-LAT data has confirmed the presence of an extended source, J1810.3−1925e, that based on its location and morphology appears to be associated to component A of HESS J1809−193. However, the spectrum of J1810.3−1925e does not connect smoothly to that of component A of HESS J1809−193, implying the need for a spectral hardening around 100 GeV. While our presented model is not able to describe this feature, we note that the overall shape of the SED is reminiscent of that of another well-known PWN system, Vela X, which also exhibits a break at around 100 GeV (Tibaldo et al. 2018). However, multiple distinct components have not been resolved at TeV energies for this system yet. Furthermore, with its characteristic age of only ~10 kyr (Manchester et al. 2005) and a very low braking index of $n = 1.4$ (Lyne et al. 1996), the Vela pulsar has an evolution history quite different from PSR J1809−1917.

Another interesting PWN to compare to is HESS J1825−137, which is the prototype of an extended (~100 pc diameter) PWN that shrinks in size at high γ-ray energies (H.E.S.S. Collaboration 2019). The pulsar powering HESS J1825−137, PSR B1823−13 (PSR J1826−1334), is quite similar to PSR J1809−1917 in terms of spin-down power

($\dot{E} = 2.8 \times 10^{36}$ erg s$^{-1}$), period ($P = 101$ ms), and distance ($d = 3.6$ kpc), but may be slightly younger (characteristic age $\tau_c = 21.4$ kyr; Manchester et al. 2005). Comparing their γ-ray PWN, HESS J1809−193 is somewhat less extended than HESS J1825−137 and does not exhibit an energy-dependent morphology. On the other hand, HESS J1809−193 seems to be composed of two distinct components, whereas HESS J1825−137 can be modelled with a single component that decreases in extent with increasing energy. This may suggest that the PWN systems have evolved differently, for example due to differences in the density of the surrounding ISM, or due to a different evolution of the corresponding pulsar (e.g. Khangulyan et al. 2018 have proposed an unusually short birth period of $P_0 \sim 1$ ms for PSR B1823−13).

Finally, it is interesting that HESS J1809−193 shows characteristics very similar to those of HESS J1702−420 (H.E.S.S. Collaboration 2021): both have been resolved into a compact, hard-spectrum component surrounded by an extended, softer-spectrum component. This may in principle suggest a similar origin of the γ-ray emission; however, HESS J1702−420 is a ‘dark’ source that lacks an obvious counterpart at other wavelengths (see also Giunti et al. 2022), hampering a further comparison with HESS J1809−193.

While we are not able to draw definitive conclusions about the origin of the γ-ray emission of HESS J1809−193, our detailed and simultaneous characterisation of its morphology and spectrum is a big step towards understanding this source. Further observations, in particular with HAWC (Abeysekara et al. 2017b) as well as with the upcoming Cherenkov Telescope Array (CTA; Acharya et al. 2018) and Southern Wide-Field Gamma-Ray Observatory (SWGO; Huentemeyer et al. 2019), will be crucial in broadening our knowledge about HESS J1809−193 further.

**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, 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. This research made use of the ASTROPY (https://www.astropy.org; Astropy Collaboration 2013, 2018), MATPLOTLIB (https:// matplotlib.org; Hunter 2007), and CORNER (https:// corner.readthedocs.io; Foreman-Mackey 2016) software packages.

## References

* Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017a, Science, 358, 911

* Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017b, ApJ, 843, 39

* Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2020, Phys. Rev. Lett., 124, 021102

* Acharya, B. S., Agudo, I., Al Samarai, I., et al. 2018, Science with the Cherenkov Telescope Array (World Scientific Publishing)

* Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, A&A, 457, 899
Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007, A&A, **472**, 489

Aharonian, F., Bergström, L., & Dermer, C. 2013, Astrophysics at Very High Energies, Swiss Society for Astrophysics and Astronomy, Saas-Fee Advanced Course (Berlin Heidelberg: Springer-Verlag), 30

Albert, A., Alfaro, R., Alvarez, C., et al. 2021a, ApJ, **911**, L27

Albert, A., Alfaro, R., Alvarez, C., et al. 2021b, ApJ, **911**, 143

Alford, J. A. J., & Halpern, J. P. 2016, ApJ, **818**, 122

Amato, E., Guetta, D., & Blasi, P. 2003, A&A, **402**, 827

Anada, T., Bamba, A., Ebisawa, K., & Dotani, T. 2010, PASJ, **62**, 179

Araya, M. 2018, ApJ, **859**, 69

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, **558**, A33

Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, **156**, 123

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

Atwood, W., Albert, A., Baldini, L., et al. 2013, 2012 Fermi Symp. proc. - eConf C121028

Ballet, J., Burnett, T. H., Digel, S. W., & Lott, B. 2020, Fermi Large Area Telescope Fourth Source Catalog Data Release 2

Bamba, A., Ueno, M., Koyama, K., & Yamauchi, S. 2003, ApJ, **589**, 253

Becherini, Y., Djannati-Ataï, A., Marandon, V., Punch, M., & Pita, S. 2011, Astropart. Phys., **34**, 858

Bell, A. R. 2013, Astropart. Phys., **43**, 56

Berezinskii, V. S., Bulanov, S. V., Dogiel, V. A., et al. 1990, in Astrophysics of Cosmic Rays, ed. V. L. Ginzburg (North-Holland)

Berge, D., Funk, S., & Hinton, J. 2007, A&A, **466**, 1219

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

Blondin, J. M., Chevalier, R. A., & Frierson, D. M. 2001, ApJ, **563**, 806

Brogan, C. L., Gelfand, J. D., Gaensler, B. M., Kassim, N. E., & Lazio, T. J. W. 2006, ApJ, **639**, L25

Cash, W. 1979, ApJ, **228**, 939

Castelletti, G., Giacani, E., & Petriella, A. 2016, A&A, **587**, A71

Cristofari, P. 2021, Universe, **7**, 324

Deil, C., Zanin, R., Lefaucheur, J., et al. 2017, in Proc. 35th Int. Cosmic Ray Conf. (ICRC2017), 766

Deil, C., Wood, M., Hassan, T., et al. 2018, https://zenodo.org/record/4701492

Deil, C., Donath, A., Terrier, R., et al. 2020, https://zenodo.org/record/1409831

de Oña Wilhelmi, E., López-Coto, R., Amato, E., & Aharonian, F. 2022, ApJ, **930**, L2

Foreman-Mackey, D. 2016, JOSS, **1**, 24

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

Gabici, S., Aharonian, F. A., & Casanova, S. 2009, MNRAS, **396**, 1629

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

Gallant, Y. A., & Arons, J. 1994, ApJ, **435**, 230

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

Gillessen, S. 2004, Ph.D. Thesis, Ruprecht-Karls-Universität Heidelberg, Germany

Ginzburg, V. L. 1975, Phil. Trans. R. Soc. Lond. A, **277**, 463

Giunti, L., Acero, F., Khelifi, B., et al. 2022, A&A, **667**, A130

Goodman, J. 2022, Presentation at Gamma 2022 Conference

Green, D. A. 2019, J. Astrophys. Astron., **40**, 36

Hahn, J. 2015, in Proc. 34th Int. Cosmic Ray Conf. (ICRC2015), 917

Hahn, J., de los Reyes, R., Bernlöhr, K., et al. 2014, Astropart. Phys., **54**, 25

H.E.S.S. Collaboration (Abdalla, H., et al.) 2018a, A&A, **612**, A1

H.E.S.S. Collaboration (Abdalla, H., et al.) 2018b, A&A, **612**, A11

H.E.S.S. Collaboration (Abdalla, H., et al.) 2019, A&A, **621**, A116

H.E.S.S. Collaboration (Abdalla, H., et al.) 2020, ApJS, **247**, 33

H.E.S.S. Collaboration (Abdalla, H., et al.) 2021, A&A, **653**, A152

H.E.S.S. Collaboration (Abdalla, H., et al.) 2022, ApJS, **260**, 53

Hinton, J. A., Funk, S., Parsons, R. D., & Ohm, S. 2011, ApJ, **743**, L7

Holler, M., Berge, D., van Eldik, C., et al. 2015, in 34th International Cosmic Ray Conference (ICRC2015), 847

Huentemeyer, P., BenZvi, S., Dingus, B., et al. 2019, Bull. Am. Astron. Soc., **51**, 109

Hunter, J. D. 2007, Comput. Sci. Eng., **9**, 90

Kargaltsev, O., & Pavlov, G. G. 2007, ApJ, **670**, 655

Kargaltsev, O., Rangelov, B., & Pavlov, G. G. 2013, The Universe Evolution: Astrophysical and Nuclear Aspects, eds. I., Strakovsky, & L., Blokhintsev (Nova Science Publishers), 359

Khangulyan, D., Koldoba, A. V., Ustyugova, G. V., Bogovalov, S. V., & Aharonian, F. 2018, ApJ, **860**, 59

Klingler, N., Kargaltsev, O., Pavlov, G. G., & Posselt, B. 2018, ApJ, **868**, 119

Klingler, N., Yang, H., Hare, J., et al. 2020, ApJ, **901**, 157

Komin, N., Carrigan, S., Djannati-Ataï, A., et al. 2007, in Proc. 30th Int. Cosmic Ray Conf. (ICRC2007)

Li, T., & Ma, Y. 1983, ApJ, **272**, 317

Linden, T., Auchettl, K., Bramante, J., et al. 2017, Phys. Rev. D, **96**, 103016

Lyne, A. G., Pritchard, R. S., Graham-Smith, F., & Camilo, F. 1996, Nature, **381**, 497

Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, **129**, 1993

Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, **461**, 396

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

Morris, D. J., Hobbs, G., Lyne, A. G., et al. 2002, MNRAS, **335**, 275

Ohm, S., van Eldik, C., & Egberts, K. 2009, Astropart. Phys., **31**, 383

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

Parthasarathy, A., Shannon, R. M., Johnston, S., et al. 2019, MNRAS, **489**, 3810

Parthasarathy, A., Johnston, S., Shannon, R. M., et al. 2020, MNRAS, **494**, 2012

Popescu, C. C., Yang, R., Tuffs, R. J., et al. 2017, MNRAS, **470**, 2539

Renaud, M., Hoppe, S., Komin, N., et al. 2008, AIP Conf. Proc., **1085**, 285

Shan, S. S., Zhu, H., Tian, W. W., et al. 2018, ApJS, **238**, 35

Sudoh, T., Linden, T., & Hooper, D. 2021, J. Cosmol. Astropart. Phys., **08**, 010

Tibaldo, L., Zanin, R., Faggioli, G., et al. 2018, A&A, **617**, A78

Umemoto, T., Minamidani, T., Kuno, N., et al. 2017, PASJ, **69**, 78

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

Voisin, F., Rowell, G. P., Burton, M. G., et al. 2019, PASA, **36**, e014

Wilks, S. S. 1938, Ann. Math. Stat., **9**, 60

Wood, M., Caputo, R., Charles, E., et al. 2017, in Proc. 35th Int. Cosmic Ray Conf. (ICRC2017), 824

Zabalza, V. 2015, in Proc. 34th Int. Cosmic Ray Conf. (ICRC2015), 922

Zarić, D., Green, D., Strzys, M., & Vovk, I. 2021, in 37th Int. Cosmic Ray Conf. (ICRC2021), 818, https://icrc2021-venue.desy.de/media/downloadAttachment/key/5db238ec30557e415d0854ce31f9abd6/maid/78

¹ Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2, Ireland

² Max-Planck-Institut für Kernphysik, PO Box 103980, 69029 Heidelberg, Germany

³ Landessternwarte, Universität Heidelberg, Königstuhl, 69117 Heidelberg, Germany

⁴ Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands

⁵ Laboratoire Leprince-Ringuet, École Polytechnique, CNRS, Institut Polytechnique de Paris, 91128 Palaiseau, France

⁶ University of Namibia, Department of Physics, Private Bag 13301, Windhoek 10005, Namibia

⁷ Centre for Space Research, North-West University, Potchefstroom 2520, South Africa

⁸ DESY, 15738 Zeuthen, Germany

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

¹⁰ Université de Paris, CNRS, Astroparticule et Cosmologie, 75013 Paris, France

¹¹ Department of Physics and Electrical Engineering, Linnaeus University, 351 95 Växjö, Sweden

¹² Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany

¹³ Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université de Paris, 92190 Meudon, France

¹⁴ Sorbonne Université, Université Paris-Diderot, Sorbonne Paris-Cité, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Énergies, LPNHE, 4 place Jussieu, 75252 Paris, France

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

¹⁶ Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Erwin-Rommel-Str. 1, 91058 Erlangen, Germany

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

¹⁸ Université Savoie Mont Blanc, CNRS, Laboratoire d'Annecy de Physique des Particules – IN2P3, 74000 Annecy, France

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

²⁰ Université Bordeaux, CNRS, LP2I Bordeaux, UMR 5797, 33170 Gradignan, France

²¹ School of Physical Sciences, University of Adelaide, Adelaide 5005, Australia

²² Aix-Marseille Université, CNRS/IN2P3, CPPM, Marseille, France
23 School of Science, Western Sydney University, Locked Bag 1797, Penrith South DC, NSW 2751, Australia

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

25 Institut für Astro- und Teilchenphysik, Leopold-Franzens- Universität Innsbruck, 6020 Innsbruck, Austria

26 Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149, 22761 Hamburg, Germany

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

28 Institute of Astronomy, Faculty of Physics, Astronomy and Infor- matics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland

29 Nicolaus Copernicus Astronomical Center, Polish Academy of Sci- ences, ul. Bartycka 18, 00-716 Warsaw, Poland

30 Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, 72076 Tübingen, Germany

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

32 Yerevan Physics Institute, 2 Alikhanian Brothers St., 375036 Yere- van, Armenia

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

34 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

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

36 Department of Physics, Rikkyo University, 3-34-1 Nishi-Ikebukuro, Toshima-ku, Tokyo 171-8501, Japan
# Appendix A: Fit of hadronic background model

We describe in this appendix the fit of the hadronic background model to the analysed observation runs. This procedure is necessary to ensure a valid description of the background for all runs.

We adjusted the background model by running the three-dimensional likelihood fit without any source components included, whereby we excluded regions around known γ-ray sources from the fit (as indicated in Fig. A.2). For each observation, we fitted the overall background normalisation ($\phi_{BG}$) and a spectral tilt parameter ($\delta_{BG}$), which modifies the predicted background rate $R_{BG}$ at energy $E$ as

$$R'_{BG} = R_{BG} \cdot (E/E_0)^{-\delta_{BG}}, \qquad (A.1)$$

with reference energy $E_0 = 1$ TeV, to correct for small inaccuracies of the spectral shape of the background model. Figure A.1 displays distributions of the fitted parameter values. In Fig. A.2, we show the resulting significance map obtained with all observations, which we have smoothed using a top hat kernel with 0.07° radius, which corresponds approximately to the size of the point spread function of H.E.S.S. for this analysis. We computed the significance for each spatial pixel following Li & Ma (1983), assuming a negligible uncertainty of the predicted number of background events, which is justified when considering the entire energy range of the analysis, as we did here. This is also known as the ‘Cash’ statistic (see Cash 1979).

While significant γ-ray emission is clearly visible in each of the exclusion regions, no significant deviation from the predicted background is present outside these regions. This picture is confirmed by the distribution of significance values in all map pixels, which we show in Fig. A.3. That the distribution for all pixels outside the exclusion regions closely follows a Gaussian distribution with unity width – the expectation for purely statistical fluctuations – indicates that we have achieved a very good description of the hadronic background. Assuming the excess in width above unity is due to a systematic effect that scales the background rate by a constant factor, the observed width of $\sigma = 1.041$ implies a level of background systematics of ~2% for the studied data set.

<figure><img src="image_12.png" /><figcaption>**Fig. A.1.** Distributions of fitted background model parameters $\phi_{BG}$ and $\delta_{BG}$. The mean ($\mu$) and standard deviation ($\sigma$) of the distributions are specified in the panels.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. A.2.** Significance map after fit of hadronic background model. We used an oversampling radius of 0.07° for smoothing. Besides HESS J1809–193, the 6° × 6° RoI contains the H.E.S.S. sources HESS J1804–216, HESS J1808–204 and HESS J1813–178. The grey dashed line marks the Galactic plane, and the black dashed circles show exclusion regions as used in the fit of the hadronic background model.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. A.3.** Significance distributions for hadronic background model fit. The black histogram shows the distribution of significance values of all spatial pixels of the map in Fig. A.2. The dark grey histogram is for all pixels outside the exclusion regions around HESS J1804–216, HESS J1808–204, and HESS J1813–178 (i.e. contains only the excess from HESS J1809–193), and the light grey histogram is for all pixels outside any of the exclusion regions. The blue line shows the result of fitting a Gaussian distribution to the light grey histogram, the fitted width of the Gaussian is given in the upper right corner of the plot.</figcaption></figure>
## Appendix B: Estimation of systematic uncertainties

We list in Table B.1 the parameters that we have used to vary the IRFs in order to estimate systematic uncertainties for our fit results. The variation of the global energy scale is implemented such that the energy axes of the effective area, energy dispersion, and point spread function IRFs are scaled by a factor $\phi_E$, simulating the effect of a wrongly calibrated energy scale. We have sampled the applied factors $\phi_E$ from a Gaussian distribution with 10% width, which approximately corresponds to a shift of the energy scale as may be expected, for example, from atmospheric variations (Hahn 2014).

For the background models, we have considered three different systematic variations: besides the overall normalisation ($\phi_{BG}$) and spectral tilt parameter ($\delta_{BG}$) that were already introduced in Appendix A, we have applied in addition a linear gradient across the field of view, with amplitude $A_{BG}^{grad}$ and direction angle $\alpha_{BG}^{grad}$. To choose the magnitude of the variation, we have compared the (already adjusted) background model to the observed data. The adopted values reflect the variation of the parameters that may still reasonably be expected for the entire data set after the adjustment described in Appendix A. A gradient of the background in the field of view may, for example, arise due to unmodelled diffuse $\gamma$-ray emission. All systematic parameters are assumed to affect the entire data set (i.e. all observation runs) in the same way.

The procedure to derive the systematic errors is then the following:

1. Sample random values for the variation parameters from the distributions indicated in Table B.1.

2. Apply systematic variations to the IRFs of all data sets.

3. Generate pseudo data sets with randomised observed counts by sampling for each bin from a Poisson distribution with mean equal to the number of events predicted by the background model and best-fit source models, using the varied IRFs to compute the latter.

4. Re-perform likelihood fit, using the original, unmodified IRFs of the data sets.

5. Repeat steps (a)–(d) 2500 times.

6. Determine the spread of the distributions of fitted source parameters, and compute systematic errors by subtracting quadratically their statistical uncertainties.

Step 6 of the above procedure is illustrated in Fig. B.1, where we show two example parameter distributions, corresponding to the amplitude $N_0$ and spectral index $\Gamma$ of source model component A (cf. Sect. 3.1). For both parameters, the spread of the distribution (indicated by the red dashed lines) exceeds the expectation from purely statistical fluctuations (grey band), which demonstrates that the considered systematic effects affect the uncertainties of these parameters.

Finally, we illustrate the correlation between the source model parameters (again exemplary those of source component A, cf. Sect. 3.1) and the systematic variation parameters in Fig. B.2.<fnref n="12" /> As intuitively makes sense, the strongest correlation is observed between the energy scale shift parameter $\phi_E$ and the source model normalisation $N_0$. The parameters $\phi_{BG}$ and $\delta_{BG}$ are – as in our standard fit procedure – re-adjusted in the fit of the pseudo data sets, hence no strong correlation with the source model parameters is expected. The parameters of the spatial part of the source model (e.g. its R.A. and Dec. coordinate) are affected mostly by the assumed background model gradient, and a modulation of these parameters with the angle $\alpha_{BG}^{grad}$ is observed.

While the correlations between the systematic variation parameters and the source model parameters are sometimes weak, they do lead to a slight broadening of the source model parameter distributions in most cases (see e.g. the distribution of the spectral index $\Gamma$ shown in Fig. 1b), and thus enable the estimation of a systematic uncertainty. We specify the resulting systematic uncertainties for all parameters along with the general fit results in Sect. 3.1, Table 2. For parameters for which the distribution was not broadened by the systematic effects considered here, we specify no systematic error – this is the case for the eccentricity ($e$) and position angle ($\phi$) parameters of source component A, as well as for the fitted R.A. and Dec. coordinate of source component B. However, we note that these parameters may be affected by other systematic effects neglected here. In particular, the fitted source positions are subject to a systematic uncertainty of the pointing position of the H.E.S.S. telescopes, which is of the order of 10'' – 20'' (Gillessen 2004).

<table><caption>**Table B.1.** Parameter variations for systematic uncertainty estimation.</caption><thead><tr><th>Par.</th><th>Variation</th><th>Description</th></tr></thead><tbody><tr><td colspan="3" style="text-align: center;">Global energy scale</td></tr><tr><td>φ<sub>E</sub></td><td>Gaussian (μ = 1, σ = 0.1)</td><td>Shift of energy scale</td></tr><tr><td colspan="3" style="text-align: center;">Background model variations</td></tr><tr><td>φ<sub>BG</sub></td><td>Gaussian (μ = 1, σ = 0.01)</td><td>Background model normalisation</td></tr><tr><td>δ<sub>BG</sub></td><td>Gaussian (μ = 0, σ = 0.02)</td><td>Background model spectral tilt</td></tr><tr><td>A<sup>grad</sup><sub>BG</sub></td><td>Gaussian (μ = 1, σ = 0.01)</td><td>Amplitude of background gradient (in deg<sup>-1</sup>)</td></tr><tr><td>α<sup>grad</sup><sub>BG</sub></td><td>Uniform (0° – 360°)</td><td>Direction of background gradient</td></tr></tbody></table>

**Notes.** For a detailed explanation of the parameters, see text.

<footnote n="12">¹² In the figure, we show the absolute value of $A_{BG}^{grad}$, and have flipped $\alpha_{BG}^{grad}$ by 180° if $A_{BG}^{grad}$ was negative.</footnote>
<figure><img src="image_15.png" /><figcaption>**Fig. B.1.** Source parameter distributions for two example parameters. (a) Flux normalisation $N_0$ of source component A. (b) Power-law index $\Gamma$ of source component A. The black dashed line and grey band denote the best-fit value and statistical uncertainty of this parameter, also specified in the top right corner. The red line shows the fit of a Gaussian to the histogram, where the fitted mean and width are indicated in the top left corner. The red dashed lines indicate the $1-\sigma$ width of the fitted Gaussian, which represents the total (i.e. statistical and systematic) uncertainty on this parameter. The derived systematic error, obtained by subtracting quadratically the statistical error from the fitted total width, is also stated in the figure.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. B.2.** Correlation plots between the source model parameters of source component A (rows) and the systematic variation parameters (columns). The value specified in the top right corner of each panel is Pearson's correlation coefficient.</figcaption></figure>
# Appendix C: Significance maps in energy bands

To check the agreement between a fitted model and the observed data, it can be illustrative to investigate residual significance maps in separate energy ranges. Furthermore, it is possible to search for indications of energy-dependent source morphology by fitting a source model within restricted energy ranges, and observing how the fitted source parameters vary between the different ranges. For these purposes, we have defined four mutually exclusive energy bands, with lower boundaries at 0.27 TeV, 0.75 TeV, 2.1 TeV, and 5.6 TeV.

## C.1. Fit of 1-component model in energy bands

As detailed in the main part of the paper, the 1-component model does not yield a satisfactory description of the observed H.E.S.S. data. This has led us to adopt the 2-component model, which provides a much better fit. We investigate in this section the possibility that a similarly good fit can be obtained by allowing the parameters of the Gaussian spatial model of the 1-component model to vary with energy. To this end, we have repeated the likelihood fit in the four different energy bands defined above, whereby we kept the spectral index Γ fixed to its best-fit value from the fit across all energies (Γ = 2.184). In Fig. C.1, we show the resulting residual significance maps for all four energy bands, with the best-fit position and 1-σ radius of the fitted component overlaid. The maps demonstrate that even when the parameters of the spatial model are allowed to vary with energy, the 1-component model cannot describe the data well – this is most

evident from the energy band between 0.75 TeV and 2.1 TeV, where strong residuals are still visible.

## C.2. Residual significance maps for 2-component model in energy bands

We provide in Fig. C.2 residual significance maps for the 2-component model, computed for the four different energy bands defined above. We note that these maps have been computed based on the 2-component model fitted across the full energy range (i.e. between 0.27 TeV and 100 TeV). The absence of strong residuals in any of the bands indicates a good agreement between the fitted model and the observed data, at all energies.

## C.3. Fit of 2-component model in energy bands

While the 2-component model yields a satisfactory description of the data across the entire energy range (cf. previous section), we also checked for a possible energy-dependent morphology of component A, by repeating the fit in the same energy bands as above. In these fits, we fixed the parameters of component B, as well as the spectral index of component A. The resulting residual significance maps are displayed in Fig. C.3. While the centre position and shape of the elongated Gaussian model vary slightly between the different energy ranges, those changes are well within the statistical uncertainties. In particular, as summarised in Table 1 in the main text, the fitted extent of component A shows no significant variation with energy.
<figure><img src="image_17.png" /></figure>

**Fig. C.1.** Residual significance maps for the 1-component model, fitted in energy bands. We used an oversampling radius of 0.07° for smoothing. Energy bands are specified in the upper right corner of each panel. The blue circle markers and dashed lines display the best-fit position and 1-sigma extent of the Gaussian models. The grey dashed line marks the Galactic plane, while the white dashed circles show regions excluded from the analysis, and the black triangle marker denotes the position of PSR J1809-1917.
<figure><img src="image_18.png" /></figure>

**Fig. C.2.** Residual significance maps for the 2-component model. We used an oversampling radius of 0.07° for smoothing. Energy bands are specified in the upper right corner of each panel. The green and purple circle marker and dashed line display the best-fit position and 1-sigma extent of the Gaussian models of component A and B, respectively. The grey dashed line marks the Galactic plane, white dashed circles shows regions excluded from the analysis, and the black triangle marker denotes the position of PSR J1809-1917.
<figure><img src="image_19.png" /></figure>

**Fig. C.3.** Residual significance maps for the 2-component model, fitted in energy bands. We used an oversampling radius of 0.07° for smoothing. Energy bands are specified in the upper right corner of each panel. The green and purple circle marker and dashed line display the best-fit position and 1-sigma extent of the Gaussian models of component A and B, respectively. The grey dashed line marks the Galactic plane, white dashed circles show regions excluded from the analysis, and the black triangle marker denotes the position of PSR J1809-1917.
**Appendix D: Fugin CO map**

We show in Fig. D.1 a map of the ¹²CO (J=1–0) emission, which traces molecular hydrogen gas, taken from the Fugin survey (Umemoto et al. 2017). The emission is shown for an interval in velocity with respect to the local standard of rest of v=16–27km s⁻¹ (following Castelletti et al. 2016), which corresponds to a distance of ~3 kpc. The map has been smoothed with a Gaussian kernel of 0.5' width. Using an CO-to-H₂ conversion factor of 1.5 × 10²⁰ cm⁻² / (K km s⁻¹), and accounting for an additional 20% of He, we obtained a gas density of ~1000 cm⁻³ in the vicinity of component B of HESS J1809–193.

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

**Fig. D.1.** Map of Fugin ¹²CO (J=1–0) line emission in the region surrounding PSR J1809–1917 (Umemoto et al. 2017). Contour lines at 40 K km s⁻¹ and 80 K km s⁻¹ are shown in black and dark grey, respectively. The lines at 40 K km s⁻¹ are also displayed in Fig. 2(b). The black triangle marker denotes the position of PSR J1809–1917, light blue circles show the positions of SNRs, and the grey dashed line marks the Galactic plane. The position and 1-σ extent of component A and component B of HESS J1809–193 are displayed in green and purple, respectively.