© The Authors 2023

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 (Ė = 6.4 × 10³⁶ erg s⁻¹, d  ∼  5 kpc) and PSR J1809−1917 (Ė = 1.8 × 10³⁶ erg s⁻¹, d  ∼  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 (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¹ – 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. 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 −2log(ℒ), where ℒ = ∏_iP(n_i | μ_i), and P(n_i | μ_i) is the Poisson probability of observing n_i events in bin i, given a predicted number of events μ_i from the background model, and γ-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² (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

(1)

with normalisation N₀, spectral index Γ, and reference energy E₀, as well as a power law with exponential cutoff (ECPL),

(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 γ-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 γ-ray emission of HESS J1809−193 by adding source components to the model prediction. For nested³ models, the preference of one model over another one can be computed from the ‘test statistic’, TS = −2log(ℒ₀/ℒ₁), which – in the limit of sufficient statistics, and far enough from parameter boundaries – follows a χ² 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 ϕ₀ 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 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 analysis. For a more general estimate of the systematic uncertainties of H.E.S.S., we refer to Aharonian et al. (2006), where systematic 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⁴ version 2.2.0 and FERMIPY⁵ 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⁶.

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

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

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.

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-σ 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σ), and in fact no residual deviations except for those expected from statistical fluctuations can be made out. In the following, we 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 ¹²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 ¹²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 ¹²CO data with the two components of HESS J1809−193 overlaid.

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.

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 μ = 0 and width σ = 1.

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 (Γ = 1.98 ± 0.05 vs. Γ = 2.24 ± 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 (Γ = 1.90 ± 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σ. For component B, on the other hand, we found no significant preference for a cutoff to the spectrum.

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

The spectrum of component A, which dominates the total emission except at the highest energies, is well compatible with 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⁸: (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,

(3)

(with log the natural logarithm); (ii) 4FGL J1811.5−1925 is modelled as a point-like source and its spectrum is fitted with 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₀ 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).

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.

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

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.

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⁹ (dashed line in Fig. 7).

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.

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 −α and an exponential cutoff at energy E_c. Its normalisation is proportional to θ × Ė, 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₀(1 + t/τ₀)^0.5, of the pulsar spin-down power as Ė(t) = Ė₀(1+t/τ₀)⁻², and of the magnetic field as B(t) = B₀[1 + (t/τ₀)^0.5]⁻¹, where τ₀ = P₀/(2Ṗ₀) 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)¹⁰. 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 (τ ≈ 33 kyr¹¹) and are associated with the extended H.E.S.S. component (A); (ii) ‘medium-age’ electrons, which have been injected within the last τ_med ≈ 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 τ_young ≈ 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 γ-ray spectrum is also well compatible with the total flux from HESS J1809−193 as measured by HAWC (Goodman 2022).

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

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

(4)

Fig. 9. Measured and predicted radius of the extended component of HESS J1809−193 (component A). As measured radius, we used the 1-σ 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 (D0, δ; cf. Eq. (4)) to the measured data points; the dashed orange and dotted green curves show results for fixed values of δ as indicated in the legend.

where E_e is the electron energy, D₀ denotes the diffusion coefficient at a reference energy of 40 TeV, and δ specifies the energy 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₀ and δ 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₀ ∼ 1.1 × 10²⁸ cm² s⁻¹ 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 δ, with both Kolgoromov scaling (δ = 1/3) and Bohm scaling (δ = 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 G011.1+00.1, existing distance estimates for G011.0−00.0 of 2.6 kpc (Bamba et al. 2003), ∼3 kpc (Castelletti et al. 2016), and 2.4 ± 0.7 kpc (Shan et al. 2018) seem broadly consistent with those for PSR J1809−1917 of 3.7 kpc (Morris et al. 2002) and 3.27 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 lead 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).

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.

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 γ-ray spectra for hadronic models fitted to each of the two H.E.S.S. components, respectively.

The pp model for component A prefers a relatively hard spectral index of Γ^p, A = 1.76 ± 0.17. Integrating the primary spectrum above 1 GeV yields a total required energy of W^p, A ∼ 3.2 × 10⁵⁰ (n/1 cm⁻³)⁻¹ 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⁵¹ 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 observed γ-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 Γ^p, B = 1.35 ± 0.45. The required energy in protons above 1 GeV is W^p, B ∼ 2.7 × 10⁴⁹ (n/1 cm⁻³)⁻¹ 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⁻³, 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 of the primary proton spectrum leads to a too strong cutoff in the γ-ray spectrum, leaving the highest HAWC flux points unexplained. We have therefore repeated the fit of the hadronic models, adding a χ² term that denotes the deviation between the sum of the predicted γ-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 (Γ^p, A = 1.95 ± 0.10; Γ^p, B = 1.56 ± 0.22) and slightly higher cutoff energies (; ). 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-σ 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-σ 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 (Ė = 2.8 × 10³⁶ erg s⁻¹), period (P = 101 ms), and distance (d = 3.6 kpc), but may be slightly younger (characteristic age τ_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₀ ∼ 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.

Acknowledgments

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

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 (ϕ_BG) and a spectral tilt parameter (δ_BG), which modifies the predicted background rate R_BG at energy E as

(A.1)

with reference energy E₀ = 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).

Fig. A.1. Distributions of fitted background model parameters ϕBG and δBG. The mean (μ) and standard deviation (σ) of the distributions are specified in the panels.

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.

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 σ = 1.041 implies a level of background systematics of ∼2% for the studied data set.

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.

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 ϕ_E, simulating the effect of a wrongly calibrated energy scale. We have sampled the applied factors ϕ_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 (ϕ_BG) and spectral tilt parameter (δ_BG) that were already introduced in Appendix A, we have applied in addition a linear gradient across the field of view, with amplitude and direction angle . 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 γ-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₀ and spectral index Γ 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.

Fig. B.1. Source parameter distributions for two example parameters. (a) Flux normalisation N0 of source component A. (b) Power-law index Γ 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-σ 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.

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.¹² As intuitively makes sense, the strongest correlation is observed between the energy scale shift parameter ϕ_E and the source model normalisation N₀. The parameters ϕ_BG and δ_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 is observed.

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.

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 Γ 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 (ϕ) 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).

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.

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, white dashed circles shows regions excluded from the analysis, and the black triangle marker denotes the position of PSR J1809−1917.

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.

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.

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.

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

Fig. D.1. Map of FUGIN 12CO (J = 1–0) line emission in the region surrounding PSR J1809−1917 (Umemoto et al. 2017). Contour lines at 40 K km s−1 and 80 K km s−1 are shown in black and dark grey, respectively. The lines at 40 K km s−1 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.

All Tables

All Figures

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.

In the text

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.

In the text

	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 μ = 0 and width σ = 1.
In the text

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

In the text

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

In the text

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

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

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

In the text

Fig. 9. Measured and predicted radius of the extended component of HESS J1809−193 (component A). As measured radius, we used the 1-σ 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 (D0, δ; cf. Eq. (4)) to the measured data points; the dashed orange and dotted green curves show results for fixed values of δ as indicated in the legend.

In the text

	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 γ-ray spectra for hadronic models fitted to each of the two H.E.S.S. components, respectively.
In the text

	Fig. A.1. Distributions of fitted background model parameters ϕBG and δBG. The mean (μ) and standard deviation (σ) of the distributions are specified in the panels.
In the text

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

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.

In the text

Fig. B.1. Source parameter distributions for two example parameters. (a) Flux normalisation N0 of source component A. (b) Power-law index Γ 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-σ 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.

In the text

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

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, white dashed circles shows regions excluded from the analysis, and the black triangle marker denotes the position of PSR J1809−1917.

In the text

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.

In the text

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 shows regions excluded from the analysis, and the black triangle marker denotes the position of PSR J1809−1917.

In the text

Fig. D.1. Map of FUGIN 12CO (J = 1–0) line emission in the region surrounding PSR J1809−1917 (Umemoto et al. 2017). Contour lines at 40 K km s−1 and 80 K km s−1 are shown in black and dark grey, respectively. The lines at 40 K km s−1 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.

In the text