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

,


Introduction
The H.E.S.S. Galactic Plane Survey (HGPS; Abdalla et al. 2018b) 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 multiwavelength counterparts, a large fraction of the sources remain without firm association.Among the firmly identified sources, 1.8 × 10 36 erg s −1 , d ∼ 3.3 kpc) (Manchester et al. 2005), which powers an X-ray PWN with an extension of ∼3 (Kargaltsev & Pavlov 2007;Anada et al. 2010;Klingler et al. 2018Klingler et al. , 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 1 -most likely provided by one of the pulsars -that up-scatter photons from ambient radiation fields to γ-ray energies via the Inverse Compton (IC) process ('leptonic scenario'), or as being due to interactions of high-energy cosmic-ray nuclei -accelerated, for example, at SNR shock fronts -within nearby molecular clouds ('hadronic scenario').In the HESS J1809−193 discovery paper (Aharonian et al. 2007), as well as in two follow-up studies presented shortly afterwards (Komin et al. 2007;Renaud et al. 2008), the authors have shown that the PWN surrounding PSR J1809−1917 can naturally explain the observed γ-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 γ-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.

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

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(Deil et al. , 2020) (v0.17)) (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 sourcefree 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.

Likelihood analysis
In the likelihood analysis, the best-fit models have been obtained by minimising the quantity −2 log(L), where L = i P(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 simulations2 (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 with normalisation N 0 , spectral index Γ, and reference energy E 0 , as well as a power law with exponential cut-off (ECPL), where E c is the cut-off 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 in Appendix A).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 (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. 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 nested3 models, the preference of one model over another one can be computed from the 'test statistic', TS = −2 log(L 0 /L 1 ), which -in the limit of sufficient statistics, and far enough from parameter boundaries -follows a χ 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 φ 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).

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 mismodelling 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 de- 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 tail 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.

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 years and 5 months of data, taken between August 4, 2008 andDecember 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 Fermitools4 version 2.2.0 and Fermipy5 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−1925eprovided 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 • × 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 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 γray emission. 6n 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−1925immediately 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 analysis7 .

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 (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.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σ), 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  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 ∼7 from the position of PSR J1809−1917 and its surrounding X-ray nebula, indicated by the brown contour lines (Anada et al. 2010).On the other hand, component B lies with its centre point directly on the western edge of the SNR G011.0+00.0, and is furthermore spatially coincident with dense molecular clouds as indicated by 12 CO (J=3-2) observations by the James Clerk Maxwell Telescope (JCMT; Castelletti et al. 2016) and CS observations by the Mopra telescope (Voisin et al. 2019).Moreover, the contour lines from the FUGIN 12 CO (J=1-0) survey (Umemoto et al. 2017) illustrate that molecular gas is present throughout the region.In Appendix D, we provide a map of the FUGIN 12 CO data with the two components of HESS J1809−193 overlaid.
Finally, we show in Fig. 4 the energy spectra and flux points obtained for the two components of HESS J1809−193, and compare these to previously obtained results from the literature.We note that, when comparing the fitted PL models for component A and B (solid lines), the spectrum of component B appears somewhat harder than that of component 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 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 (Γ = 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 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 • radius, and thus part of the larger-scale emission has been missed.
We summarise the best-fit parameter values found for the 2component 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.

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:8 (i) 4FGL J1810.3−1925e is modelled as an extended source (using a two-dimensional Gaussian as spatial model) and its spectrum is fitted with a logparabola model, (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−1925eshows no significant curvature, but have decided to maintain the logparabola model for consistency with the 4FGL-DR2 catalogue (which extended to lower energies than our analysis here).Removing J1810.3−1925e and J1811.5−1925from the bestfit ROI model, we obtained the significance map shown in Fig. 5, Notes.σ, e, and φ denote the 1-σ radius, eccentricity, and position angle of the Gaussian spatial model, respectively.N 0 , Γ, E 0 , and E c are parameters of the spectral models, as defined in Eqs. ( 1) and ( 2).Systematic uncertainties have been derived as described in Appendix B. (a) For the asymmetric component A, σ refers to the semi-major axis.
The extents of the semi-minor axis compute to 0.353 ± 0.029 stat and 0.351 ± 0.028 stat for the PL and ECPL spectral model, respectively.
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 posi- Notes.σ is the 1-σ radius of the Gaussian spatial model.N 0 , α, β, Γ, and E 0 are parameters of the spectral models defined in Eqs. ( 3) and ( 1).
tion 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−1925efrom 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 sensitivity9 (dashed line in Fig. 7).

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

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 am-   5−1925 (this work) disk model (Araya 2018) 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.bient 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(Parthasarathy et al. , 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 0 (1+t/τ 0 ) 0.5 , of the pulsar spin-down power as Ė(t) = Ė0 (1 + t/τ 0 ) −2 , and of the  magnetic field as B(t) = B 0 [1+(t/τ 0 ) 0.5 ] −1 , where τ 0 = P 0 /(2 Ṗ0 ) is the initial spin-down time scale.We then computed the nonthermal 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)  10 .Finally, we have fitted the adjustable parameters of the model to the observed spectral energy distribution (SED) of HESS J1809−193, where the optimisation has been carried out using a Markov chain Monte Carlo (MCMC) 10 We note that the prediction of this large-scale model may not be very accurate in the specific region studied here.As our conclusions are based on order-of-magnitude estimates, however, they are unaltered even if the predicted radiation field densities are wrong by a factor of a few.
Article number, page 8 of 23 F. Aharonian et al.: HESS J1809−193: a halo of escaped electrons around a pulsar wind nebula?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 (τ ≈ 33 kyr11 ) 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 mea-sured 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 ∼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).
The model fails, however, to explain the spectrum of the Fermi-LAT source J1810.3−1925ebelow ∼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−1925eexhibits 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. (2018Klingler et al. ( , 2020) ) have detected a northward proper motion of ∼20-40 mas yr −1 , albeit not with high significance.This would imply a travel time between the bestfit 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 where E e is the electron energy, D 0 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 0 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 0 ∼ 1.1 × 10 28 cm 2 s −1 appears reasonable and is of the same order of magnitude as the coefficient obtained for the Geminga halo by Abeysekara et al. (2017a).On the other hand, the observed data do not provide very strong constraints for δ, with both Kolgoromov scaling (δ = 1/3) and Bohm scaling (δ = 1) Article number, page 9 of 23 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).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 (D 0 , δ; 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.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 γ-ray spectrum is expected to occur.The measured cut-off 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).

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) and3.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(Klingler et al. , 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−1925ecoincides 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 Article number, page 10 of 23 F. Aharonian et al.: HESS J1809−193: a halo of escaped electrons around a pulsar wind nebula?Table 5. Best-fit parameter values for the hadronic pp models.

Par. [unit] Value
Component A N p,A 0 [10 34 eV −1 ] 10.0 ± 1.8 Γ p,A 1.76 ± 0.17 Notes.N 0 , Γ, E 0 , and E c are parameters of the ECPL spectral model (Eq.2), where the superscript p, X denotes that these are the parameters of the primary proton spectrum of component X = {A, B}, respectively.The quoted errors represent statistical uncertainties only.
the steep spectrum of the former to the harder spectra of the latter, a spectral hardening in the primary cosmic-ray spectrum, for which there is no obvious explanation.As presented in Sect.4.1, both components of HESS J1809−193 can be modelled well within a leptonic scenario.Nevertheless, we have also explored the implications of either of the components being of hadronic origin.To this end, we have fitted a proton-proton (pp) model to both components, employing the Naima package (Zabalza 2015).The primary proton spectrum is described using an ECPL model (see Eq. 2) and we have assumed a distance to the source of 3 kpc.We used the wrapper class for Naima models implemented in Gammapy, so that they could be fitted directly to the H.E.S.S. data sets (as opposed to fitting them to the extracted flux points only), using the same likelihood framework as before (cf.Sect.2.1.2).The fit results are presented in Table 5, and the resulting spectra displayed in Fig. 10.The same spatial models as in the previous analysis (cf.Sect.3.1) were assumed and compatible best-fit parameters were obtained for them.
The pp model for component A prefers a relatively hard spectral index of Γ 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 50 (n/1 cm −3 ) −1 erg, which represents -unless high ISM densities are invoked -a significant fraction of the canonically assumed kinetic energy released in a supernova explosion of ∼10 51 erg (e.g.Ginzburg 1975).In this context, we note that while the FUGIN CO data indeed show the presence of molecular gas beyond the dense clouds discovered by Castelletti et al. (2016) and Voisin et al. (2019), they also indicate a gradient in the gas density across the extent of HESS J1809−193 (see the map in Appendix D).This gradient is not reflected in the 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 49 (n/1 cm −3 ) −1 erg.Even considering that only part of the cosmic rays potentially accelerated by G011.0−00.0 will reach the dense molecular clouds, the high density of gas in the vicinity of component B in general (∼1000 cm −3 , cf.Appendix D) makes this energy seem well affordable.An explanation of component B of HESS J1809−193 in a hadronic scenario therefore appears entirely reasonable,   2022).The lines show the predicted γ-ray spectra for hadronic models fitted to each of the two H.E.S.S. components, respectively.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 cut-off energy E p,B c ∼ 110 TeV of the primary proton spectrum leads to a too strong cut-off in the γ-ray spectrum, leaving the highest HAWC flux points unexplained.We have therefore repeated the fit of the hadronic models, adding a χ 2 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 (E p,A c = 140 +80 −50 TeV; E p,B c = 200 +420 −130 TeV).These values are, however, consistent within uncertainties with those obtained in the previous fit, which demonstrates that it is possible to also explain the measured HAWC flux points within a hadronic scenario (or a mixed one, in which component A is leptonic and component B is hadronic).
Finally, we note that while the relatively hard primary proton spectra obtained for both components of HESS J1809−193 are not consistent with generic predictions of diffuse shock acceleration (Bell 2013), they are compatible with a scenario in which cosmic rays accelerated in a supernova remnant illuminate a nearby gas cloud (e.g.Gabici et al. 2009).On the other hand, there is also the possibility of a continuous wind of hadronic cosmic rays powered by the pulsar (e.g.Gallant & Arons 1994;Amato et al. 2003), which may be an interesting scenario to explore for the case of PSR J1809−1917, as already noted by Voisin et al. (2019).

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 Article number, page 11 of 23 A&A proofs: manuscript no.paper_hess_j1809_193 a 1-σ semi-major and semi-minor axis of ∼ 0.62 • and ∼ 0.35 • , respectively, and exhibits a spectrum with a cut-off 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 cut-off.
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−1925edoes 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 (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 ( Ė = 2.8 × 10 36 erg s −1 ), 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 energydependent 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 ∼ 1 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 γ-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; Abreu et al. 2019), will be crucial in further broadening our knowledge about HESS J1809−193.
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.

Fig. 2 .
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 fromCastelletti 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 .
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.1(b), while the blue, filled histogram corresponds to the map in Fig.1(c).The black, dashed line displays a Gaussian distribution with a mean µ = 0 and width σ = 1.

Fig. 4 .
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), Abdalla et al. (2018b), Zarić et al. (2021), Abeysekara et al. (2020), and Goodman (2022).

Fig. 5 .
Fig. 5. Significance maps above 1 GeV from the Fermi-LAT analysis.(a) With J1811.5−1925 and J1810.3−1925eremoved from the best-fit model.(b) With J1810.3−1925eremoved 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. 7 .
Fig. 7. Fermi-LAT energy spectra obtained with H.E.S.S. model templates.The dashed grey line shows the 10-year Fermi-LAT broadband sensitivity.

Fig. 8 .
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 fromAnada 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).
Fig.9.Measured and predicted radius of the extended component of HESSJ1809−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. Table1).The solid blue curve has been obtained by fitting the electron diffusion parameters (D 0 , δ; 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.

Fig. 10 .
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 fromGoodman (2022).The lines show the predicted γ-ray spectra for hadronic models fitted to each of the two H.E.S.S. components, respectively.
Article number, page 2 of 23 F. Aharonian et al.: HESS J1809−193: a halo of escaped electrons around a pulsar wind nebula?

Table 1
the fitted extent of component A in the four employed energy bands.

Table 1 .
Extent of H.E.S.S. component A in energy bands.Notes.E min and E max denote the lower and upper boundary of the energy band, respectively, E mean the weighted mean energy.σ major and σ minor denote the 1-sigma extent of the semi-major and semi-minor axis of the elongated Gaussian spatial model for H.E.S.S. component A, respectively.See Appendix C.3 for further details.

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.

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

Table 4 .
Parameters of the PWN modelNotes.Pulsar parameters denote present-day values unless otherwise specified.The parameter values for the 'adjusted' parameters were obtained using an MCMC method, but should be regarded as indicative values rather than precise fit results (see main text).