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

TL;DR

This study resolves gamma-ray emission around pulsar PSR J1809-1917 into extended and compact components, suggesting a halo of escaped electrons and complex origins involving the pulsar wind nebula, supernova remnant, and molecular clouds.

Contribution

First detailed spatial and spectral analysis of HESS J1809-193 revealing a halo of escaped electrons and multiple emission components.

Findings

Extended component exhibits a spectral cutoff at ~13 TeV.

Compact component is located near the pulsar with no clear spectral cutoff.

Extended gamma-ray emission is confirmed by Fermi-LAT analysis.

Abstract

Context. HESS J1809$-$193 is an unassociated very-high-energy $\gamma$-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 $\gamma$-ray emission up to energies of $\sim$100 TeV with the HAWC observatory has led to renewed interest in HESS J1809$-$193. Aims. We aim to understand the origin of the $\gamma$-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 that exhibits a spectral cut-off at $\sim$13 TeV, and a compact component that is located close to PSR J1809$-$1917 and shows no clear spectral cut-off. The Fermi-LAT analysis also revealed extended $\gamma$-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. (abridged)

Tables

Table 1. Table 1: Extent of H.E.S.S. component A in energy bands.

$E_{\min }-E_{\max }$	$E_{\mathrm{mean}}$	${\sigma }_{\mathrm{major}}$	${\sigma }_{\mathrm{minor}}$
(TeV)	(TeV)	(deg)	(deg)
$0.27-0.75$	0.43	$0.69\pm 0.10$	$0.41\pm 0.10$
$0.75-2.1$	1.2	$0.62\pm 0.05$	$0.31\pm 0.04$
$2.1-5.6$	3.2	$0.62\pm 0.06$	$0.35\pm 0.05$
$>5.6$	9.6	$0.61\pm 0.11$	$0.46\pm 0.13$

Table 2. Table 2: Best-fit parameter values for the H.E.S.S. 2-component model. For component A, we provide the best-fit values for the assumption of a PL spectral model and an ECPL spectral model.

Par. [unit]	Value
Component A (PL spectral model)
R.A. [deg]	$272.551\pm {0.025}_{\mathrm{stat}}\pm {0.018}_{\mathrm{sys}}$
	(${18}^h{10}^m{12}^s\pm 6_{\mathrm{stat}}^s\pm 4_{\mathrm{sys}}^s$)
Dec. [deg]	$-19.344\pm {0.023}_{\mathrm{stat}}\pm {0.013}_{\mathrm{sys}}$
	($-{19}^{\circ }{20.6}^{\prime }\pm {1.4}_{\mathrm{stat}}^{\prime }\pm {0.8}_{\mathrm{sys}}^{\prime })$
$\sigma$ [deg]	$0.622\pm {0.032}_{\mathrm{stat}}\pm {0.020}_{\mathrm{sys}}$a𝑎aa𝑎a$a$For the asymmetric component A, $\sigma$ refers to the semi-major axis. The extents of the semi-minor axis compute to $0.353\pm {0.029}_{\mathrm{stat}}$ and $0.351\pm {0.028}_{\mathrm{stat}}$ for the PL and ECPL spectral model, respectively.
$e$	$0.824\pm {0.025}_{\mathrm{stat}}$
$\varphi$ [deg]	$50.0\pm {3.1}_{\mathrm{stat}}$
$N_0\left[{10}^{-12}{\mathrm{TeV}}^{-1}{\mathrm{cm}}^{-2}{s}^{-1}\right]$	$8.42\pm {0.40}_{\mathrm{stat}}\pm {1.14}_{\mathrm{sys}}$
$\Gamma$	$2.239\pm {0.027}_{\mathrm{stat}}\pm {0.020}_{\mathrm{sys}}$
$E_0$ [TeV]	1 (fixed)
Component A (ECPL spectral model)
R.A. [deg]	$272.554\pm {0.025}_{\mathrm{stat}}\pm {0.019}_{\mathrm{sys}}$
	(${18}^h{10}^m{13}^s\pm 6_{\mathrm{stat}}^s\pm 5_{\mathrm{sys}}^s$)
Dec. [deg]	$-19.344\pm {0.021}_{\mathrm{stat}}\pm {0.012}_{\mathrm{sys}}$
	($-{19}^{\circ }{20.6}^{\prime }\pm {1.3}_{\mathrm{stat}}^{\prime }\pm {0.7}_{\mathrm{sys}}^{\prime })$
$\sigma$ [deg]	$0.613\pm {0.031}_{\mathrm{stat}}\pm {0.015}_{\mathrm{sys}}$a𝑎aa𝑎a$a$For the asymmetric component A, $\sigma$ refers to the semi-major axis. The extents of the semi-minor axis compute to $0.353\pm {0.029}_{\mathrm{stat}}$ and $0.351\pm {0.028}_{\mathrm{stat}}$ for the PL and ECPL spectral model, respectively.
$e$	$0.820\pm {0.025}_{\mathrm{stat}}$
$\varphi$ [deg]	$51.3\pm {3.1}_{\mathrm{stat}}$
$N_0\left[{10}^{-12}{\mathrm{TeV}}^{-1}{\mathrm{cm}}^{-2}{s}^{-1}\right]$	$9.05\pm {0.47}_{\mathrm{stat}}\pm {0.91}_{\mathrm{sys}}$
$\Gamma$	$1.90\pm {0.05}_{\mathrm{stat}}\pm {0.05}_{\mathrm{sys}}$
$E_c$ [TeV]	${12.7}_{-2.1}^{+2.7}{|}_{\mathrm{stat}}{}_{-1.9}{}^{+2.6}|_{\mathrm{sys}}^{}$
$E_0$ [TeV]	1 (fixed)
Component B
R.A. [deg]	$272.400\pm {0.010}_{\mathrm{stat}}$
	(${18}^h{09}^m{36}^s\pm {2.4}_{\mathrm{stat}}^s$)
Dec. [deg]	$-19.406\pm {0.009}_{\mathrm{stat}}$
	($-{19}^{\circ }{24.4}^{\prime }\pm {0.5}_{\mathrm{stat}}^{\prime })$
$\sigma$ [deg]	$0.0953\pm {0.0072}_{\mathrm{stat}}\pm {0.0034}_{\mathrm{sys}}$
$N_0\left[{10}^{-12}{\mathrm{TeV}}^{-1}{\mathrm{cm}}^{-2}{s}^{-1}\right]$	$0.95\pm {0.11}_{\mathrm{stat}}\pm {0.011}_{\mathrm{sys}}$
$\Gamma$	$1.98\pm {0.05}_{\mathrm{stat}}\pm {0.03}_{\mathrm{sys}}$
$E_0$ [TeV]	1 (fixed)

Table 3. Table 3: Best-fit parameter values for the Fermi -LAT data analysis.

Par. [unit]	Value
J1810.3$-$1925e
R.A. [deg]	$272.547\pm {0.033}_{\mathrm{stat}}$
	(${18}^h{10}^m{11}^s\pm 8_{\mathrm{stat}}^s$)
Dec. [deg]	$-19.397\pm {0.038}_{\mathrm{stat}}$
	($-{19}^{\circ }{23.8}^{\prime }\pm {2.3}_{\mathrm{stat}}^{\prime })$
$\sigma$ [deg]	$0.317\pm {0.024}_{\mathrm{stat}}$
$N_0\left[{10}^{-9}{\mathrm{GeV}}^{-1}{\mathrm{cm}}^{-2}{s}^{-1}\right]$	$2.86\pm {0.27}_{\mathrm{stat}}\pm {0.09}_{\mathrm{sys}}$
$\alpha$	$2.53\pm {0.11}_{\mathrm{stat}}$
$\beta$	$-0.015\pm {0.043}_{\mathrm{stat}}$
$E_0$ [GeV]	1.747 (fixed)
J1811.5$-$1925
R.A. [deg]	$272.874\pm {0.025}_{\mathrm{stat}}$
	(${18}^h{11}^m{30}^s\pm 6_{\mathrm{stat}}^s$)
Dec. [deg]	$-19.410\pm {0.027}_{\mathrm{stat}}$
	($-{19}^{\circ }{24.6}^{\prime }\pm {1.6}_{\mathrm{stat}}^{\prime })$
$N_0\left[{10}^{-9}{\mathrm{GeV}}^{-1}{\mathrm{cm}}^{-2}{s}^{-1}\right]$	$0.016\pm {0.003}_{\mathrm{stat}}\pm {0.001}_{\mathrm{sys}}$
$\Gamma$	$2.42\pm {0.16}_{\mathrm{stat}}$
$E_0$ [GeV]	7.48 (fixed)

Table 4. Table 4: Parameters of the PWN model

Par.	Description	Value
Input parameters
$d$	pulsar distancea𝑎aa𝑎a$a$Taken from Manchester et al. (2005).	$3.3\mathrm{kpc}$
$\dot{E}$	pulsar spin-down powera𝑎aa𝑎a$a$Taken from Manchester et al. (2005).	$1.8\text{×}{10}^{36}\mathrm{erg}{s}^{-1}$
${\tau }_c$	pulsar characteristic agea𝑎aa𝑎a$a$Taken from Manchester et al. (2005).	$51.4\mathrm{kyr}$
$P$	pulsar period a𝑎aa𝑎a$a$Taken from Manchester et al. (2005).	$82.76\mathrm{ms}$
$\dot{P}$	pulsar period derivativea𝑎aa𝑎a$a$Taken from Manchester et al. (2005).	$2.55\text{×}{10}^{-14}s{s}^{-1}$
$n$	pulsar braking indexb𝑏bb𝑏b$b$Assumed value.	3
Adjusted parameters
$\theta$	electron power fraction	$0.6$
$B$	magnetic field	$4\mu G$
$P_0$	pulsar birth period	$50\mathrm{ms}$
$E_c$	cut-off energy	$420\mathrm{TeV}$
$\alpha$	injection spectrum index	$2.0$
${\tau }_{\mathrm{young}}$	age of young e-	$1.2\mathrm{kyr}$
${\tau }_{\mathrm{med}}$	age of medium-age e-	$4.7\mathrm{kyr}$

Table 5. Table 5: Best-fit parameter values for the hadronic p ​ p 𝑝 𝑝 pp models.

Par. [unit]	Value
Component A
$N_0^{p,A}\left[{10}^{34}{\mathrm{eV}}^{-1}\right]$	$10.0\pm 1.8$
${\Gamma }^{p,A}$	$1.76\pm 0.17$
$E_c^{p,A}$ [TeV]	${90}_{-30}^{+35}$
$E_0^{p,A}$ [TeV]	20 (fixed)
Component B
$N_0^{p,B}\left[{10}^{34}{\mathrm{eV}}^{-1}\right]$	$1.0\pm 0.4$
${\Gamma }^{p,B}$	$1.34\pm 0.45$
$E_c^{p,B}$ [TeV]	${110}_{-50}^{+135}$
$E_0^{p,B}$ [TeV]	20 (fixed)

Table 6. Table 6: Parameter variations for systematic uncertainty estimation.

Par.	Variation	Description
Global energy scale
${\varphi }_E$	Gaussian ($\mu =1,\sigma =0.1$)	Shift of energy scale
Background model variations
${\varphi }_{\mathrm{BG}}$	Gaussian ($\mu =1,\sigma =0.01$)	Background model normalisation
${\delta }_{\mathrm{BG}}$	Gaussian ($\mu =0,\sigma =0.02$)	Background model spectral tilt
$A_{\mathrm{BG}}^{\mathrm{grad}}$	Gaussian ($\mu =1,\sigma =0.01$)	Amplitude of background gradient (in deg-1)
${\alpha }_{\mathrm{BG}}^{\mathrm{grad}}$	Uniform ($0^{\circ }-{360}^{\circ }$)	Direction of background gradient

Equations

$$\frac{\mathrm{d}N}{\mathrm{d}E}=N_0\cdot {\left(\frac{E}{E_0}\right)}^{-\Gamma },$$

$$\frac{\mathrm{d}N}{\mathrm{d}E}=N_0\cdot {\left(\frac{E}{E_0}\right)}^{-\Gamma }\cdot \exp \left(-\frac{E}{E_c}\right),$$

$$\frac{\mathrm{d}N}{\mathrm{d}E}=N_0\cdot {\left(\frac{E}{E_0}\right)}^{-\alpha -\beta \log (E/E_0)}$$

$$D=D_{0}\left(\frac{E_{e}}{$40\text{\,}\mathrm{TeV}$}\right)^{\delta}\,,$$

$$R_{\mathrm{BG}}^{\prime }=R_{\mathrm{BG}}\cdot (E/E_0{)}^{-{\delta }_{\mathrm{BG}}},$$

Full text

Aims. We aim to understand the origin of the $\gamma$-ray emission of HESS J1809$-$193.

Methods. We analysed $93.2\text{\,}\mathrm{h}$ of data taken on HESS J1809$-$193 above $0.27\text{\,}\mathrm{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\text{\,}\mathrm{yr}$ of Fermi-LAT data above $1\text{\,}\mathrm{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-$\sigma$ semi-major / semi-minor axis of $\sim$0.62∘ / $\sim$0.35∘) that exhibits a spectral cut-off at \sim$$13\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}, and a compact component (modelled as a symmetric Gaussian with a 1-$\sigma$ radius of $\sim$0.1∘) that is located close to PSR J1809$-$1917 and shows no clear spectral cut-off. The Fermi-LAT analysis also revealed extended $\gamma$-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.

1 Introduction

The H.E.S.S. Galactic Plane Survey (HGPS; Abdalla et al., 2018b) has revealed a large number of Galactic $\gamma$-ray sources in the very-high-energy (VHE; $E>$100\text{,}\mathrm{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 ($\dot{E}=$6.4\text{\times}{10}^{36}\text{,}\mathrm{erg}\text{,}{\mathrm{s}}^{-1}, $d\sim$5\text{\,}\mathrm{kpc}) and, noteworthy in particular, PSR J1809$-$1917 ($\dot{E}=$1.8\text{\times}{10}^{36}\text{,}\mathrm{erg}\text{,}{\mathrm{s}}^{-1}, $d\sim$3.3\text{\,}\mathrm{kpc}) (Manchester et al., 2005), which powers an X-ray PWN with an extension of $\sim$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 $\gamma$-ray emission as originating from high-energy electrons111In this paper, we use the term ‘electrons’ to refer to both electrons and positrons. – most likely provided by one of the pulsars – that up-scatter photons from ambient radiation fields to $\gamma$-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 have shown that the PWN surrounding PSR J1809$-$1917 can naturally explain the observed $\gamma$-ray emission in a leptonic scenario, and thus 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 $\gamma$-ray emission from HESS J1809$-$193 with the High Altitude Water Cherenkov Observatory (HAWC) above $56\text{\,}\mathrm{TeV}$ (and quite possibly above $100\text{\,}\mathrm{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 $\gamma$-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 $\gamma$-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\text{\,}\mathrm{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 $\gamma$ 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\text{\,}\mathrm{m}$ above sea level. The original array, referred to as ‘HESS-I’, was installed in 2000-2003 and comprised four telescopes with $12\text{\,}\mathrm{m}$-diameter mirrors (CT1-4), arranged in a square layout with $120\text{\,}\mathrm{m}$ side length (Aharonian et al., 2006). The array was completed in 2012 by a fifth telescope, CT5, featuring a $28\text{\,}\mathrm{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 $\gamma$ rays in the energy range between \sim$$0.1\text{\,}\mathrm{TeV} and \sim$$100\text{\,}\mathrm{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\text{\,}\mathrm{min}$ duration. Selecting runs that encompass HESS J1809$-$193 within $\sim$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\text{\,}\mathrm{h}$. This represents a significant increase with respect to the previous dedicated publications on HESS J1809$-$193, which used $25\text{\,}\mathrm{h}$ (Aharonian et al., 2007), $32\text{\,}\mathrm{h}$ (Komin et al., 2007), and $41\text{\,}\mathrm{h}$ (Renaud et al., 2008) of data.

In the data analysis, we have selected $\gamma$ 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 $\gamma$ rays by several orders of magnitude and they cannot be rejected completely in the event selection without a too severe loss of $\gamma$-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\text{\,}\mathrm{TeV}$ for all observations, while the lowest threshold, obtained for $\sim$10% of the observations, is $0.27\text{\,}\mathrm{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^{\circ}$) 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 simulations222Technically, 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. (Bernlöhr, 2008). As spatial source models we have used 2-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

[TABLE]

with normalisation $N_{0}$, spectral index $\Gamma$, and reference energy $E_{0}$, as well as a power law with exponential cut-off (ECPL),

[TABLE]

where $E_{c}$ is the cut-off energy.

We carried out the likelihood fit in a region of interest (RoI) of 6° $\times$ 6°, centred on HESS J1809$-$193 (see Fig. 12 in Appendix A). For the binning of our data, we used spatial pixels of 0.02° $\times$ 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 (Abdalla et al., 2018b), and HESS J1808$-$204 (Abdalla et al., 2018a), which we have masked in the fit using circular exclusion regions (cf. Fig. 12).

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

models, the preference of one model over another one can be computed from the ‘test statistic’, $\mathrm{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 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 $\gamma$ rays in the energy range from \sim$$20\text{\,}\mathrm{MeV} to several hundred GeV (Atwood et al., 2009). Here, we analysed 12 years 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\text{\,}\mathrm{GeV}$, we restricted the analysis to events above this energy. To suppress $\gamma$ rays originating from the Earth’s limb, we furthermore excluded events with zenith angles above $90^{\circ}$. The data were analysed using Fermitools444https://fermi.gsfc.nasa.gov/ssc/data/analysis/software version 2.2.0 and Fermipy 555https://fermipy.readthedocs.io version 1.1.5 (Wood et al., 2017).

The analysis was carried out in a region of interest (ROI) of $10^{\circ}\times 10^{\circ}$, centred on the nominal position of 4FGL J1810.3$-$1925e provided in the 4FGL-DR2 catalogue (Abdollahi et al., 2020; Ballet et al., 2020). The events and exposure maps were binned using spatial bins of $0.1^{\circ}\times 0.1^{\circ}$ size and five bins per decade in energy. $\gamma$ ray source models were then fitted to the data with a likelihood fit as described in Section 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 $\gamma$-ray emission.666We 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.

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^{\circ}$ 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 analysis777We have followed the procedure outlined at https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/Aeff_Systematics.html.

3 Results

3.1 H.E.S.S. results

We show in Fig. 1(a) the residual significance map after the fit of the hadronic background model, while the map in Fig. 2(a) displays the deduced flux of $\gamma$ rays above the threshold energy of $0.27\text{\,}\mathrm{TeV}$. Extended $\gamma$-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^{\circ}$, 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. 1(a), and in Fig. 1(b) 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. Fig. 1(a) 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. 1(c), 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 will 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. 2(b), we provide a detail 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 $\sim$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 12CO ($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 12CO ($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 12CO 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 component A and B (solid lines), the spectrum of component B appears somewhat harder than that of component A ($\Gamma=1.98\pm 0.05$ vs. $\Gamma=2.24\pm 0.03$). However, the flux upper limits obtained above energies of \sim$$20\text{\,}\mathrm{TeV} for component A seem to indicate the presence of a cut-off 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 \sim$$5\text{\,}\mathrm{TeV} (dashed line in Fig. 4). With respect to the PL spectral model, the ECPL model is preferred with a statistical significance of $\sim$8$\sigma$. For component B, on the other hand, we found no significant preference for a cut-off to the spectrum.

The spectrum of component A, which dominates the total emission except at the highest energies, is well compatible with that published in Abdalla et al. (2018b), 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^{\circ}$ 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 (Abdollahi et al., 2020; Ballet et al., 2020) lists two sources that are located in the immediate vicinity of HESS J1809$-$193:101010We 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 (Abdollahi et al., 2022). (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,

[TABLE]

(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 will 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_{0}$ 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.

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 sensitivity121212Broadband sensitivity curves for Fermi-LAT are available at https://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance. We have used the curve for Galactic coordinates $l=0^{\circ}$ and $b=0^{\circ}$. (dashed line in Fig. 7).

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 $\gamma$-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 (Section 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 (Section 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 $\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)131313We 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.. 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 life time of the system ($\tau\approx$33\text{,}\mathrm{kyr}151515The ‘true’ age of the pulsar can be computed as $\tau=(P/((n-1)\cdot\dot{P}))\cdot(1-(P_{0}/P)^{n-1})$ (Gaensler & Slane, [2006](#bib.bib42)). 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](#S4.T4)). Other values of $P_{0}$ will lead to different estimates of the pulsar age. ) and are associated with the extended H.E.S.S. component (A); (ii) ‘medium-age’ electrons, which have been injected within the last $\tau_{\mathrm{med}}\approx$4.7\text{\,}\mathrm{kyr} and are associated with the compact H.E.S.S. component (B); (iii) ‘young’ electrons, which have been injected within the last $\tau_{\mathrm{young}}\approx$1.2\text{,}\mathrm{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\text{\,}\mathrm{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 / 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 \sim$$4\text{\,}\mathrm{\SIUnitSymbolMicro G} and a reasonable spectral index for the injection spectrum of $\sim$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 \sim$$10\text{\,}\mathrm{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 Section 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 \sim$$2040\text{\,}\mathrm{mas}\text{\,}{\mathrm{yr}}^{-1}, 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 \sim$$1020\text{\,}\mathrm{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

[TABLE]

where $E_{e}$ is the electron energy, $D_{0}$ denotes the diffusion coefficient at a reference energy of $40\text{\,}\mathrm{TeV}$, and $\delta$ 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 $\gamma$-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\text{\times}{10}^{28}\text{,}{\mathrm{cm}}^{2}\text{,}{\mathrm{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](#bib.bib7)). 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 cut-off to the corresponding $\gamma$-ray spectrum is expected to occur. The measured cut-off energy of $\sim13\text{,}\mathrm{TeV}$for component A of HESS J1809$-$193 is well in line with this prediction, as can be seen in Fig. 8.

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 $\gamma$-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\text{\,}\mathrm{kpc}$ (Bamba et al., 2003), \sim$$3\text{\,}\mathrm{kpc} (Castelletti et al., 2016), and $2.4\pm 0.7\,$\mathrm{kpc}$$ (Shan et al., 2018) seem broadly consistent with those for PSR J1809$-$1917 of $3.7\text{\,}\mathrm{kpc}$ (Morris et al., 2002) and $3.27\text{\,}\mathrm{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 \sim$$10\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}) 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\text{\,}\mathrm{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,\mathrm{A}}=1.76\pm 0.17$. Integrating the primary spectrum above $1\text{\,}\mathrm{GeV}$ yields a total required energy of $W^{p,\mathrm{A}}\sim$3.2\text{\times}{10}^{50}$\,(n/$1\text{,}{\mathrm{cm}}^{-3}$)^{-1}\,$\mathrm{erg}, which represents – unless high ISM densities are invoked – a significant fraction of the canonically assumed kinetic energy released in a supernova explosion of $\sim{10}^{51}\text{,}\mathrm{erg} (e.g. Ginzburg, [1975](#bib.bib46)). 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](#bib.bib33)) and Voisin et al. ([2019](#bib.bib81)), they also indicate a gradient in the gas density across the extent of HESS J1809-193 (see the map in Appendix [D](#A4)). This gradient is not reflected in the 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,\mathrm{B}}=1.35\pm 0.45$. The required energy in protons above $1\text{\,}\mathrm{GeV}$ is $W^{p,\mathrm{B}}\sim$2.7\text{\times}{10}^{49}$\,(n/$1\text{,}{\mathrm{cm}}^{-3}$)^{-1}\,$\mathrm{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 ($\sim1000\text{,}{\mathrm{cm}}^{-3}, cf. Appendix [D](#A4)) 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\text{,}\mathrm{TeV}measured with HAWC. However, as is evident from Fig. [10](#S4.F10), the fitted cut-off energyE_{c}^{p,\mathrm{B}}\sim110\text{\,}\mathrm{TeV}$$ of the primary proton spectrum leads to a too strong cut-off 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,\mathrm{A}}=1.95\pm 0.10$;$\Gamma^{p,\mathrm{B}}=1.56\pm 0.22$) and slightly higher cut-off energies ($E_{c}^{p,\mathrm{A}}=140_{-50}^{+80},\mathrm{TeV}$$; E_{c}^{p,\mathrm{B}}=200_{-130}^{+420},$\mathrm{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 $\gamma$-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 $\sim 0.62^{\circ}$ and $\sim 0.35^{\circ}$, respectively, and exhibits a spectrum with a cut-off at \sim$$13\text{\,}\mathrm{TeV}. Superimposed, component B appears symmetric and more compact with a 1-$\sigma$ radius of $\sim 0.1^{\circ}$, and shows a harder spectrum with no clear cut-off.

We have interpreted the results in a leptonic scenario, in which the $\gamma$-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 $\gamma$-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 $\gamma$-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 $\gamma$-ray emission up to $100\text{\,}\mathrm{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\text{\,}\mathrm{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\text{\,}\mathrm{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 \sim$$10\text{\,}\mathrm{k}\mathrm{y}\mathrm{r} (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 (\sim$$100\text{\,}\mathrm{pc} diameter) PWN that shrinks in size at high $\gamma$-ray energies (Abdalla et al., 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\text{\times}{10}^{36}\text{,}\mathrm{erg}\text{,}{\mathrm{s}}^{-1}), period ($P=$101\text{\,}\mathrm{ms}), and distance ($d=$3.6\text{,}\mathrm{kpc}), but may be slightly younger (characteristic age $\tau_{c}=$21.4\text{\,}\mathrm{kyr}) (Manchester et al., 2005). Comparing their $\gamma$-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\text{,}\mathrm{ms}$$ for PSR B1823$-$13).

Finally, it is interesting that HESS J1809$-$193 shows characteristics very similar to those of HESS J1702$-$420 (Abdalla et al., 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 $\gamma$-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 $\gamma$-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; Abreu et al., 2019), will be crucial in further broadening our knowledge about HESS J1809$-$193.

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 Astropy171717https://www.astropy.org (Robitaille et al., 2013; Price-Whelan et al., 2018), Matplotlib181818https://matplotlib.org (Hunter, 2007), and Corner191919https://corner.readthedocs.io (Foreman-Mackey, 2016) software packages.

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 3D likelihood fit without any source components included, whereby we excluded regions around known $\gamma$-ray sources from the fit (as indicated in Fig. 12). For each observation, we fitted the overall background normalisation ($\phi_{\mathrm{BG}}$) and a spectral tilt parameter ($\delta_{\mathrm{BG}}$), which modifies the predicted background rate $R_{\mathrm{BG}}$ at energy $E$ as

[TABLE]

with reference energy $E_{0}=1\,\mathrm{TeV}$, to correct for small inaccuracies of the spectral shape of the background model. Figure 11 displays distributions of the fitted parameter values. In Fig. 12, we show the resulting significance map obtained with all observations, which we have smoothed using a top hat kernel with $0.07^{\circ}$ 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 $\gamma$-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. 13. 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 $\sim$2% for the studied data set.

Appendix B Estimation of systematic uncertainties

We list in Table 6 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 et al. 2014).

For the background models, we have considered three different systematic variations: besides the overall normalisation ($\phi_{\mathrm{BG}}$) and spectral tilt parameter ($\delta_{\mathrm{BG}}$) that were already introduced in Appendix A, we have applied in addition a linear gradient across the field of view, with amplitude $A^{\mathrm{grad}}_{\mathrm{BG}}$ and direction angle $\alpha^{\mathrm{grad}}_{\mathrm{BG}}$. 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. (a)

Sample random values for the variation parameters from the distributions indicated in Table 6. 2. (b)

Apply systematic variations to the IRFs of all data sets. 3. (c)

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

Re-perform likelihood fit, using the original, unmodified IRFs of the data sets. 5. (e)

Repeat steps (a)–(d) $2500$ times. 6. (f)

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

Step (f) of the above procedure is illustrated in Fig. 14, 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. 15.202020 In the figure, we show the absolute value of $A^{\mathrm{grad}}_{\mathrm{BG}}$, and have flipped $\alpha^{\mathrm{grad}}_{\mathrm{BG}}$ by $180^{\circ}$ if $A^{\mathrm{grad}}_{\mathrm{BG}}$ was negative.

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_{\mathrm{BG}}$ and $\delta_{\mathrm{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^{\mathrm{grad}}_{\mathrm{BG}}$ 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. 14), 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^{\prime\prime}-20^{\prime\prime}$ (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\text{\,}\mathrm{TeV}$, $0.75\text{\,}\mathrm{TeV}$, $2.1\text{\,}\mathrm{TeV}$, and $5.6\text{\,}\mathrm{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 $\Gamma$ fixed to its best-fit value from the fit across all energies ($\Gamma=2.184$). In Fig. 16, we show the resulting residual significance maps for all four energy bands, with the best-fit position and 1-$\sigma$ 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\text{\,}\mathrm{TeV}$ and $2.1\text{\,}\mathrm{TeV}$, where strong residuals are still visible.

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

We provide in Fig. 17 residual significance maps for the 2-component model, computed for the four different energy bands defined above. Note that these maps have been computed based on the 2-component model fitted across the full energy range (i.e. between $0.27\text{\,}\mathrm{TeV}$ and $100\text{\,}\mathrm{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. 18. 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.

Appendix D FUGIN CO map

We show in Fig. 19 a map of the 12CO ($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=$1627\text{,}\mathrm{km}\text{,}{\mathrm{s}}^{-1} (following Castelletti et al. [2016](#bib.bib33)), which corresponds to a distance of $\sim3\text{,}\mathrm{kpc}$. The map has been smoothed with a Gaussian kernel of$0.5\arcmin{}$width. Using an CO-to-H2 conversion factor of$1.5\times 10^{20},\mathrm{cm}^{-2},/,(\mathrm{K},\mathrm{km},\mathrm{s}^{-1})$, and accounting for an additional 20$\sim$$1000\text{,}{\mathrm{cm}}^{-3}$in the vicinity of component B of HESS J1809$-$193.