# A deep spectromorphological study of the γ-ray emission surrounding the young massive stellar cluster Westerlund 1

F. Aharonian, H. Ashkar, M. Backes, V. Barbosa Martins, Y. Becherini, D.

Berge, B. Bi, M. Böttcher, M. de Bony de Lavergne, F. Bradascio, et al.

► To cite this version:

F. Aharonian, H. Ashkar, M. Backes, V. Barbosa Martins, Y. Becherini, et al.. A deep spectromorphological study of the γ-ray emission surrounding the young massive stellar cluster Westerlund 1. Astron.Astrophys., 2022, 666, pp.A124. 10.1051/0004-6361/202244323. hal-03771735

HAL Id: hal-03771735

https://hal.science/hal-03771735v1

Submitted on 12 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
A deep spectromorphological study of the γ-ray emission
surrounding the young massive stellar cluster Westerlund 1

F. Aharonian<sup>1,2</sup>, H. Ashkar<sup>3</sup>, M. Backes<sup>4,5</sup>, V. Barbosa Martins<sup>6</sup>, Y. Becherini<sup>7,8</sup>, D. Berge<sup>6,9</sup>, B. Bi<sup>10</sup>, M. Böttcher<sup>5</sup>, M. de Bony de Lavergne<sup>11</sup>, F. Bradascio<sup>12</sup>, R. Brose<sup>1</sup>, F. Brun<sup>12</sup>, T. Bulik<sup>13</sup>, C. Burger-Scheidlin<sup>1</sup>, F. Cangemi<sup>14</sup>, S. Caroff<sup>14</sup>, S. Casanova<sup>15</sup>, M. Cerruti<sup>7</sup>, T. Chand<sup>5</sup>, S. Chandra<sup>5</sup>, A. Chen<sup>16</sup>, O. Chibueze<sup>5</sup>, P. Cristofari<sup>17</sup>, J. Damascene Mbarubucyeye<sup>6</sup>, A. Djannati-Atai<sup>7</sup>, J.-P. Ernenwein<sup>18</sup>, K. Feijen<sup>19</sup>, G. Fichet de Clairfontaine<sup>17</sup>, G. Fontaine<sup>3</sup>, S. Funk<sup>20</sup>, S. Gabici<sup>7</sup>, Y. A. Gallant<sup>21</sup>, S. Ghafourizadeh<sup>22</sup>, G. Giavitto<sup>6</sup>, L. Giunti<sup>7,12</sup>, D. Glawion<sup>20</sup>, J. F. Glicenstein<sup>12</sup>, P. Goswami<sup>5</sup>, M.-H. Grondin<sup>23</sup>, L. K. Härer<sup>2</sup>, M. Haupt<sup>6</sup>, J. A. Hinton<sup>2</sup>, M. Hörbe<sup>24</sup>, W. Hofmann<sup>2</sup>, T. L. Holch<sup>6</sup>, M. Holler<sup>25</sup>, D. Horns<sup>26</sup>, M. Jamrozy<sup>27</sup>, V. Joshi<sup>20</sup>, I. Jung-Richardt<sup>20</sup>, E. Kasai<sup>4</sup>, K. Katarzyński<sup>28</sup>, U. Katz<sup>20</sup>, B. Khélifi<sup>7</sup>, W. Kluźniak<sup>29</sup>, Nu. Komin<sup>16</sup>, K. Kosack<sup>12</sup>, D. Kostunin<sup>6</sup>, G. Kukec Mezek<sup>8</sup>, R. G. Lang<sup>20</sup>, S. Le Stum<sup>18</sup>, A. Lemière<sup>7</sup>, M. Lemoine-Goumard<sup>23</sup>, J.-P. Lenain<sup>14</sup>, F. Leuschner<sup>10</sup>, T. Lohse<sup>9</sup>, A. Luashvili<sup>17</sup>, I. Lypova<sup>22</sup>, J. Mackey<sup>1</sup>, J. Majumdar<sup>6</sup>, D. Malyshev<sup>20</sup>, V. Marandon<sup>2</sup>, P. Marchegiani<sup>16</sup>, A. Marcowith<sup>21</sup>, G. Martí-Devesa<sup>25</sup>, R. Marx<sup>22</sup>, G. Maurin<sup>11</sup>, M. Meyer<sup>26</sup>, A. Mitchell<sup>20,2</sup>, R. Moderski<sup>29</sup>, L. Mohrmann<sup>2,*</sup>, A. Montanari<sup>12</sup>, E. Moulin<sup>12</sup>, J. Muller<sup>3</sup>, T. Murach<sup>6</sup>, K. Nakashima<sup>20</sup>, M. de Naurois<sup>3</sup>, A. Nayerhoda<sup>15</sup>, J. Niemiec<sup>15</sup>, S. Ohm<sup>6,*</sup>, L. Olivera-Nieto<sup>2</sup>, E. de Ona Wilhelmi<sup>6</sup>, M. Ostrowski<sup>27</sup>, S. Panny<sup>25</sup>, M. Panter<sup>2</sup>, R. D. Parsons<sup>9</sup>, G. Peron<sup>2</sup>, D. A. Prokhorov<sup>30</sup>, G. Pühlhofer<sup>10</sup>, M. Punch<sup>7,8</sup>, A. Quirrenbach<sup>22</sup>, R. Rauth<sup>25,*</sup>, P. Reichherzer<sup>12</sup>, A. Reimer<sup>25</sup>, O. Reimer<sup>25</sup>, M. Renaud<sup>21</sup>, B. Reville<sup>2</sup>, F. Rieger<sup>2</sup>, G. Rowell<sup>19</sup>, B. Rudak<sup>29</sup>, E. Ruiz-Velasco<sup>2</sup>, V. Sahakian<sup>31</sup>, H. Salzmann<sup>10</sup>, D. A. Sanchez<sup>11</sup>, A. Santangelo<sup>10</sup>, M. Sasaki<sup>20</sup>, F. Schüssler<sup>12</sup>, H. M. Schutte<sup>5</sup>, U. Schwanke<sup>9</sup>, J. N. S. Shapopi<sup>4</sup>, A. Specovius<sup>20,*</sup>, S. Spencer<sup>24</sup>, Ł. Stawarz<sup>27</sup>, R. Steenkamp<sup>4</sup>, S. Steinmassl<sup>2</sup>, C. Steppa<sup>32</sup>, I. Sushch<sup>5</sup>, H. Suzuki<sup>33</sup>, T. Takahashi<sup>34</sup>, T. Tanaka<sup>33</sup>, R. Terrier<sup>7</sup>, C. Thorpe-Morgan<sup>10</sup>, M. Tsirou<sup>2</sup>, N. Tsiuji<sup>35</sup>, R. Tuffs<sup>2</sup>, T. Unbehaun<sup>20</sup>, C. van Eldik<sup>20</sup>, B. van Soelen<sup>36</sup>, M. Vecchi<sup>30</sup>, J. Veh<sup>20</sup>, C. Venter<sup>5</sup>, J. Vink<sup>30</sup>, S. J. Wagner<sup>22</sup>, R. White<sup>2</sup>, A. Wierzcholska<sup>15</sup>, Y. Wun Wong<sup>20</sup>, M. Zacharias<sup>17,5</sup>, D. Zargaryan<sup>1</sup>, A. A. Zdziarski<sup>29</sup>, S. J. Zhu<sup>6</sup>, S. Zouari<sup>7</sup>, N. Żywucka<sup>5</sup> (H.E.S.S. Collaboration), R. Blackwell<sup>19</sup>, C. Braiding<sup>37</sup>, M. Burton<sup>38</sup>, K. Cubuk<sup>38</sup>, M. Filipović<sup>39</sup>, N. Tothill<sup>39</sup>, and G. Wong<sup>39</sup>

(Affiliations can be found after the references)

Received 22 June 2022 / Accepted 21 July 2022

ABSTRACT

**Context.** Young massive stellar clusters are extreme environments and potentially provide the means for efficient particle acceleration. Indeed, they are increasingly considered as being responsible for a significant fraction of cosmic rays (CRs) that are accelerated within the Milky Way. Westerlund 1, the most massive known young stellar cluster in our Galaxy, is a prime candidate for studying this hypothesis. While the very-high-energy γ-ray source HESS J1646–458 has been detected in the vicinity of Westerlund 1 in the past, its association could not be firmly identified.

**Aims.** We aim to identify the physical processes responsible for the γ-ray emission around Westerlund 1 and thus to understand the role of massive stellar clusters in the acceleration of Galactic CRs better.

**Methods.** Using 164 h of data recorded with the High Energy Stereoscopic System (H.E.S.S.), we carried out a deep spectromorphological study of the γ-ray emission of HESS J1646–458. We furthermore employed H I and CO observations of the region to infer the presence of gas that could serve as target material for interactions of accelerated CRs.

**Results.** We detected large-scale (~2° diameter) γ-ray emission with a complex morphology, exhibiting a shell-like structure and showing no significant variation with γ-ray energy. The combined energy spectrum of the emission extends to several tens of TeV, and it is uniform across the entire source region. We did not find a clear correlation of the γ-ray emission with gas clouds as identified through H I and CO observations.

**Conclusions.** We conclude that, of the known objects within the region, only Westerlund 1 can explain the majority of the γ-ray emission. Several CR acceleration sites and mechanisms are conceivable and discussed in detail. While it seems clear that Westerlund 1 acts as a powerful particle accelerator, no firm conclusions on the contribution of massive stellar clusters to the flux of Galactic CRs in general can be drawn at this point.

**Key words.** acceleration of particles – radiation mechanisms: non-thermal – shock waves – stars: massive – gamma rays: general – galaxies: star clusters: individual: Westerlund 1

<footnote>* Corresponding authors; e-mail: contact.hess@hess-experiment.eu</footnote>
# 1. Introduction

Young massive stellar clusters are environments of copious star formation, and they typically host a large number of very massive stars (Portegies Zwart et al. 2010). For this reason, they have long been considered as potential sites of cosmic-ray (CR) acceleration (Parizot et al. 2004). The acceleration may take place at shock fronts of supernova remnants (SNRs), which may collide with the strong winds of massive stars in the cluster (see e.g. Bykov et al. 2020, and references therein), or at the termination shock of the superbubble that is excelled by the combined stellar winds of the cluster (e.g. Bykov 2014; Gupta et al. 2018; Morlino et al. 2021). Massive clusters form from correspondingly massive molecular clouds, which are not very common in the Milky Way, but often found in starburst galaxies (Fujii & Portegies Zwart 2016). Nevertheless, the notion that massive star clusters are responsible for the majority of hadronic CRs<fnref n="1" /> that are accelerated within our Galaxy represents a viable alternative to the long-standing ‘SNR paradigm’, in which (isolated) SNRs are the primary acceleration sites (Portegies Zwart et al. 2010; Aharonian et al. 2019; Morlino et al. 2021).

Through interaction with ambient gas and radiation fields, high-energy hadronic CRs produce high-energy $\gamma$ rays, which provides strong motivation for the search for $\gamma$-ray emission from massive stellar clusters (this has already been realised long ago, see e.g. Cesarsky & Montmerle 1983). Indeed, a bubble, or ‘cocoon’ in the Cygnus X star-forming region has been detected in Fermi-LAT data in the ~1 GeV–few 100 GeV energy range (Ackermann et al. 2011). Subsequently, Fermi-LAT has detected diffuse $\gamma$-ray emission in the same energy range around a number of other massive stellar clusters in the Milky Way (Yang & Aharonian 2017; Yang et al. 2018; Yang & Wang 2020; Sun et al. 2020a,b). Searches at higher energies (i.e. at 1 TeV and above) with ground-based instruments have also led to several detections: (i) the Cygnus region harbours several sources of TeV-energy $\gamma$ rays (Abdo et al. 2007; Abeysekara et al. 2018), and it has recently been detected up to energies of hundreds of TeV (Abeysekara et al. 2021; Cao et al. 2021); (ii) the young stellar cluster Westerlund 2 within the star formation region RCW 49, which hosts with WR 20a an extraordinarily massive binary star system (Aharonian et al. 2007; Abramowski et al. 2011); (iii) the young stellar cluster Westerlund 1 (Abramowski et al. 2012), which is introduced in more detail below; (iv) the super bubble 30 Dor C, whose detection is particularly noteworthy due to its distant location in the Large Magellanic Cloud (Abramowski et al. 2015); (v) the stellar cluster Cl* 1806–20, which contains both a luminous blue variable candidate, LBV 1806–20, and a magnetar, SGR 1806–20 (Abdalla et al. 2018a). However, the $\gamma$-ray emission can be directly linked to the stellar clusters in only some of the above cases, the precise CR acceleration sites are unidentified as of yet, and none of the detections constitutes unequivocal evidence for the acceleration of hadronic CRs due to the respective clusters. The assertion of the latter point is complicated by the fact that high-energy $\gamma$ rays can be produced by CRs via two competing processes. Besides their production in the decay of neutral pions (and other short-lived particles), produced in turn in interactions of hadronic CRs with ambient matter – the ‘hadronic scenario’ – they may also be created through the inverse Compton (IC) process, in which high-energy electrons

<table><caption>**Table 1.** Parameter values for Westerlund 1 assumed in this work.</caption><thead><tr><th>Par.</th><th>Description</th><th>Value</th><th>Ref.</th></tr></thead><tbody><tr><td>d</td><td>Distance from Earth</td><td>3.9 kpc</td><td>(1), (2)</td></tr><tr><td>τ</td><td>Cluster age</td><td>4 Myr</td><td>(3), (4)</td></tr><tr><td>L<sub>w</sub></td><td>Kinetic luminosity of cluster wind</td><td>10<sup>39</sup> erg s<sup>-1</sup></td><td>(5)</td></tr><tr><td>v<sub>w</sub></td><td>Velocity of cluster wind</td><td>3000 km s<sup>-1</sup></td><td>(6)</td></tr></tbody></table>

**References.** (1) Kothes & Dougherty (2007); (2) Davies & Beasor (2019); (3) Clark et al. (2005); (4) Brandner et al. (2008); (5) Muno et al. (2006a); (6) Morlino et al. (2021).

and positrons<fnref n="2" /> can up-scatter low-energy photons from ambient radiation fields to TeV energies – the ‘leptonic scenario’. These two scenarios can only be distinguished by carrying out detailed spectromorphological studies of the $\gamma$-ray emission, and combining the results with those obtained at other wavelengths. In this article, we present such a study for the young massive stellar cluster Westerlund 1.

Westerlund 1, named after its discoverer Bengt Westerlund (Westerlund 1961) and located at RA (2000) = 16h47m04.0s, Dec (2000) = -45°51'04.9" (Brandner et al. 2008), is the most massive known young stellar cluster in the Milky Way, with an estimated mass of around $10^5 M_{\odot}$ (Clark et al. 2005; Brandner et al. 2008; Portegies Zwart et al. 2010). It hosts a rich population of evolved massive stars, including significant fractions of all known Galactic Yellow Hypergiants (Clark et al. 2005) and Wolf-Rayet stars (Crowther et al. 2006). The half-mass radius of the cluster is approximately 1 pc (Brandner et al. 2008). Many estimates for the age of the cluster and its distance from Earth have been put forward in the past. Most age estimates agree with an age of 3–5 Myr (Clark et al. 2005; Crowther et al. 2006; Brandner et al. 2008), although the single-age paradigm has been questioned recently after finding that the observed luminosities of cool supergiants are more consistent with an age of ~10 Myr (Beasor et al. 2021). Early distance estimates, using various techniques, find distances of around 4 kpc (Clark et al. 2005; Crowther et al. 2006; Kothes & Dougherty 2007; Brandner et al. 2008). Recently, data from the Gaia spacecraft were used to obtain new distance estimates. While most of them are compatible with the old estimates (Davies & Beasor 2019; Rate et al. 2020; Beasor et al. 2021; Negueruela et al. 2022), closer distances of ~2.7 kpc have also been obtained (Aghakhanloo et al. 2020, 2021). Clark et al. (2019) have questioned the reliability of Gaia (DR2) data in the Westerlund 1 field altogether, rendering the new estimates somewhat uncertain. For this article, we adopted an age of 4 Myr and a distance of 3.9 kpc, as these values are compatible with the majority of published results. Additionally, we will need in the course of this paper estimates for the properties of the collective cluster wind of Westerlund 1, which is formed as a superposition of the strong winds of the massive stars in the cluster. Only few estimates can be found in the literature; we assumed here for the kinetic luminosity of the wind $L_w \sim 10^{39}$ erg s$^{-1}$ (Muno et al. 2006a) and for the wind velocity $v_w \sim 3000$ km s$^{-1}$ (Morlino et al. 2021) as typical values. All parameter values for Westerlund 1 assumed in this work are summarised in Table 1.

Westerlund 1 has been studied extensively in the X-ray domain. Observations with the Chandra telescope have revealed diffuse hard X-ray emission from the core of the cluster

<footnote n="1">¹ Here and in the following, the term ‘hadronic cosmic ray’ refers to cosmic-ray nuclei, as opposed to cosmic-ray electrons and positrons, for example.</footnote>

<footnote n="2">² Hereafter, we use the term ‘electrons’ to refer to both electrons and positrons.</footnote>
(Muno et al. 2006a), which was later identified as likely being of thermal origin with XMM-Newton observations (Kavanagh et al. 2011). Additionally, Muno et al. (2006b) have identified an X-ray magnetar, CXOU J164710.2–455216, which supposedly was created in the explosion of a very massive (>40 $M_\odot$) progenitor star (Clark et al. 2008; Belczynski & Taam 2008). Interestingly, CXOU J164710.2–455216 is the only known remnant of a stellar explosion within Westerlund 1.

Moving to larger spatial scales (i.e. beyond the bounds of the cluster itself), an analysis of Fermi-LAT data between 3 and 300 GeV by Ohm et al. (2013) revealed extended γ-ray emission in the vicinity of Westerlund 1. With more data accumulated since then, the latest Fermi-LAT source catalogue (4FGL-DR2, Abdollahi et al. 2020; Ballet et al. 2020) lists six sources within 1.1° from the cluster centre. Besides the stellar cluster and its members, several objects that are located at relatively small angular separations from Westerlund 1 could potentially be contributing to the observed γ-ray emission. This includes two energetic ($E > 2 \times 10^{35}$ erg s⁻¹) pulsars, PSR J1648–4611 and PSR J1650–4601 (Manchester et al. 2005), as well as the low-mass X-ray binary (LMXB) 4U 1642–45 (Forman et al. 1978). On the other hand, it is quite possible that some of the six sources listed in the 4FGL catalogue share a common physical origin, and were separated into distinct components only because the true source morphology is very complex.

Finally, at even higher energies, Abramowski et al. (2012) detected a large, extended (~2° diameter) emission region centred on Westerlund 1, named HESS J1646–458, between 0.45 and ~20 TeV with the H.E.S.S. experiment. Based on the properties of the emission and taking into account multi-wavelength data, the authors found Westerlund 1 to be the most likely explanation of the γ-ray emission in a single-source scenario, but were unable to draw definitive conclusions based on the data set available at the time.

Since then, the exposure collected with H.E.S.S. on HESS J1646–458 has almost quintupled, in large part thanks to a dedicated observation campaign in 2017. This, together with recent advances in analysis techniques, enabled a new, detailed study of the γ-ray emission surrounding Westerlund 1, which we present here. In Sect. 2.1, we introduce the H.E.S.S. data set and provide a description of the data analysis. The results of the H.E.S.S. data analysis are given in Sect. 2.2. Besides H.E.S.S. data, we have also analysed data from HI and CO observations in the vicinity of Westerlund 1 and present the results in Sect. 2.3. A detailed discussion of the results, considering multiple explanations for the observed γ-ray emission, is presented in Sect. 3, before we summarise our findings and provide an outlook in Sect. 4.

## 2. Observations and data analysis

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

H.E.S.S. is an array of five imaging atmospheric Cherenkov telescopes (IACTs), located in the Southern Hemisphere in Namibia (23°16′18″ S, 16°30′00″ E) at an altitude of 1800 m above sea level (Aharonian et al. 2006; Holler et al. 2015). The array comprises four 12 m-diameter telescopes (CT1-4) that are arranged on a 120 m-side square and began operation in late 2003. A fifth telescope (CT5), with 28 m diameter, was added in the centre of the array in 2012. The telescopes detect the Cherenkov light produced in atmospheric air showers initiated by primary γ rays, where the main background consists of showers caused by hadronic CR. With the central telescope included, the array

is sensitive to γ rays in the energy range between ~0.1 TeV and ~100 TeV; with CT5 alone thresholds as low as 20 GeV have been achieved in studies of pulsed emission (Abdalla et al. 2018b).

The H.E.S.S. data set for HESS J1646–458 comprises 362 observation runs after quality selection, taken over the course of more than 13 years between June 18, 2004 and October 11, 2017. The runs amount to a total observation time of 164.2 h, which represents a significant increase with respect to the previous publication (Abramowski et al. 2012) (33.8 h). We note that not all of the observations have targeted HESS J1646–458 directly; some have been taken as part of surveys, and some were primarily targeted at the nearby sources HESS J1640–465 (Abramowski et al. 2014a) and HESS J1641–463 (Abramowski et al. 2014b), leading to a gradient in exposure across the HESS J1646–458 region.

Only data from the four smaller telescopes (CT1-4) were considered in the analysis presented here. For best performance, we restricted the maximum zenith angle of the analysed observation runs to <60° and the maximum offset angle between the reconstructed direction of events and the telescope pointing direction to <2°. With this selection, an energy threshold of 0.37 TeV was achieved in the final analysis. γ-like events were selected using the method described in Ohm et al. (2009) and their energy and arrival direction were reconstructed with the algorithm presented in Parsons & Hinton (2014). Subsequently, we converted our data to the FITS-based data format described in Deil et al. (2018) and performed the high-level data analysis using the GAMMAPY package (Deil et al. 2017, 2020; v0.17). All findings were confirmed with two cross-check analyses: one based on a completely independent calibration and data analysis chain (de Naurois & Rolland 2009), and one based on the same calibration and event reconstruction algorithms as those used in the main analysis, but carried out with the CTOOLS package (Knödlseder et al. 2016; v.1.6.3); the latter analysis is documented in Specovius (2021). Furthermore, results of another, intermediate analysis of the data set, which has inspired parts of the analysis presented here, can be found in Zorn (2019).

In the high-level analysis, we employed a concept that has only recently been established for the analysis of IACT data: a 3-dimensional likelihood analysis, in which the data can be modelled simultaneously in two spatial dimensions and as a function of energy (Mohrmann et al. 2019). In this method, contrary to more established ones, the residual background of CR-induced air shower events ('hadronic' background) for a given observation run is not directly estimated from source-free regions within the observed field itself, but rather provided by a background model. This model was constructed from archival observations and subsequently adjusted to the analysed observations following the procedure outlined in Mohrmann et al. (2019). The 3-dimensional likelihood analysis method is especially suited for the analysis of complex source regions and largely extended sources, and is thus a suitable choice for the analysis of HESS J1646–458.

As a first step in the analysis, separate energy thresholds were determined for each observation run, requiring that the energy reconstruction bias is below 10% and that the background model is used only above its validity threshold (Mohrmann et al. 2019). Aiming for sufficient exposure across the entire region down to the lowest energies, a minimal energy threshold of 0.37 TeV was enforced. We subsequently computed 3-dimensional maps of the observed number of events, predicted background, and exposure, comprised of spatial pixels of 0.02° × 0.02° and an energy axis of 16 bins per decade of
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Significance maps after background subtraction. The position of Westerlund 1 is marked by the black star symbol; the grey, dashed line shows the Galactic plane. **a:** map for the entire 6° × 6° region of interest, smoothed with a 0.07° top-hat kernel. The final exclusion map is shown in black, earlier iterations in blue, green, purple, and orange. Locations of previously detected sources that are not connected to HESS J1646-458 are indicated by black, open symbols. **b:** map with detail view of the emission surrounding Westerlund 1, smoothed with a 0.22° top-hat kernel. The colour scale is saturated at the maximum observed significance value associated with the HESS J1646-458 region. Contour lines corresponding to a significance of 4, 8, and 12σ are shown in blue. Signal regions a–p used for spectrum extraction are overlaid (black), as are segments 1–5 for the computation of radial profiles (white dashed).</figcaption></figure>

energy, with the 6° × 6° region of interest centred on Westerlund 1. For each observation, we adjusted the background model to the observed data by fitting a global normalisation and a spectral tilt parameter, taking into account only regions in the field of view that are free of γ-ray emission. In order to safely exclude regions with γ-ray emission from the background fit, an iterative procedure as described in Abdalla et al. (2018c) was employed to generate an exclusion map.

## 2.2. H.E.S.S. analysis results

### 2.2.1. Maps and radial profiles

We show in Fig. 1 the resulting residual significance maps after background subtraction, where the significance was computed following Li & Ma (1983). For the map in Fig. 1a, a top-hat smoothing with a kernel of radius 0.07° has been employed. The corresponding distribution of significance values – for all pixels and those outside the exclusion map, displayed in black – is presented in Fig. 2. As expected for the case of purely statistical fluctuations, the distribution for pixels outside the exclusion map follows very closely that of a Gaussian distribution with unit width, indicating a good description of the hadronic background. Figure 1b shows a significance map obtained with a correlation radius of 0.22°. Overlaid are 16 square ‘signal regions’ (of size 0.45° × 0.45° each), labelled a–p, that cover the γ-ray emission of HESS J1646-458, as well as 5 circular segments. These signal regions and segments are used in the further characterisation of the γ-ray emission (see below).

Having obtained a satisfying description of the residual hadronic background, we computed flux maps displaying the excess γ-ray emission (see Fig. 3). Focusing on the larger-scale structure of the emission of HESS J1646-458, a top-hat smoothing with a kernel of radius 0.22° – the same value as already used in Abramowski et al. (2012) – has been applied for the maps in

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Significance entry distribution. The black histogram corresponds to all pixels of the map shown in Fig. 1a with non-zero entries, the grey histogram to all pixels outside of the final exclusion map. The green line represents the fit of a Gaussian distribution to the grey histogram, the best-fit mean and width are $μ = -0.043 ± 0.005$ and $σ = 1.008 ± 0.005$, respectively. A Gaussian distribution with mean $μ = 0$ and width $σ = 1$ is shown by the orange, dashed line for comparison.</figcaption></figure>

panels a, c, and d. Additionally, we show in panel b a flux map smoothed with a Gaussian kernel of 0.07° width, which corresponds approximately to the size of the point-spread function for the analysis configuration employed here.

Very strong γ-ray emission is observed from the known, nearby sources HESS J1640-465 and HESS J1641-463. Turning to HESS J1646-458, we observe that its spatial morphology is very complex. Notably, the emission is not peaked at the position of Westerlund 1, but rather exhibits a structure resembling that of a shell, surrounding the stellar cluster. This global structure is present in all of the displayed maps, and thus does not seem to vary with γ-ray energy. On top of the large-scale
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Flux maps of the HESS J1646-458 region. The position of Westerlund 1 is marked by the black star symbol; the grey, dashed line shows the Galactic plane. Coloured symbols indicate objects listed in the legend in panel a. Dark grey square markers denote positions of sources from the 4FGL-DR2 catalogue (Abdollahi et al. 2020; Ballet et al. 2020), where those sources that are still significant ($\sqrt{\text{TS}} > 3$) above 30 GeV are shown with a diamond marker (◇). Grey circles labelled ‘A’ and ‘B’ mark regions defined in Abramowski et al. (2012); region ‘C’ (at RA $16^{h}49^{m}4.8^{s}$, Dec $-46^{o}06'00''$) is newly defined here. The white circle marker indicates the coordinate with respect to which the radial profiles in Figs. 4 and 9a have been computed. The scale bar denotes a projected distance of 40 pc, for the nominal distance to Westerlund 1 of 3.9 kpc. The maps are for different energy thresholds, indicated at the bottom of each panel. The maps in panels a, c, and d were computed using a 0.22° top-hat smoothing kernel, while the map in panel b was computed using a 0.07° Gaussian smoothing kernel. Colour scales are saturated at the maximum observed flux value associated with the HESS J1646-458 region. Contour lines shown in blue are at flux levels of $F = (12.5/20/27.5) \times 10^{-9} \text{ cm}^{-2} \text{ s}^{-1} \text{ sr}^{-1}$ for panels a and b, at $F = (3/5.5/8) \times 10^{-9} \text{ cm}^{-2} \text{ s}^{-1} \text{ sr}^{-1}$ for panel c, and at $F = (1/1.5) \times 10^{-9} \text{ cm}^{-2} \text{ s}^{-1} \text{ sr}^{-1}$ for panel d.</figcaption></figure>

structure, we identify peaks of the emission in the circular regions labelled ‘A’ and ‘B’, which correspond to regions with enhanced emission already found in Abramowski et al. (2012), confirming these findings. Additionally, we find a peak – visible in particular in the flux map above 4.9 TeV – in region ‘C’, which encompasses the two energetic pulsars PSR J1648-4611 and PSR J1650-4601, as well as an emission region extending beyond the shell-like structure, to the east of region C. For future reference, we provide the following source names for these regions: HESS J1645-455 (region A), HESS J1647-465 (region B), HESS J1649-460 (region C), and HESS J1652-462 (emission east of region C). We stress, however, that we have found – as will be detailed in the course of this paper – no indications that the γ-ray emission in these regions is of a different

origin than that of the rest of the emission, and that the regions should therefore not be interpreted as distinct sources. Rather, the regions have been labelled in order to ease the discussion of the source morphology.

Because HESS J1646-458 is located along the Galactic plane, and towards the inner Galaxy, it is safe to assume that diffuse γ-ray emission contributes to the observed signal to some degree. This diffuse emission is produced by CRs that propagate freely within the Galactic disc, and can be due to bremsstrahlung or IC emission of CR electrons, or interactions of hadronic CRs with gas. Due to its diffuse nature, the diffuse γ-ray emission from the Galaxy is challenging to measure directly, and while it has been detected over large scales in the TeV energy range (e.g. Abramowski et al. 2014c; Amenomori et al. 2021), these
<table><caption>**Table 3.** Spectral analysis results for signal regions.</caption><thead><tr><th>Signal region</th><th>Excess events</th><th>Significance</th><th>Significance (E > 4.9 TeV)</th><th>φ<sub>0</sub> (10<sup>-13</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>Γ</th><th>√ΔTS</th></tr></thead><tbody><tr><td>a</td><td>396.1</td><td>5.3σ</td><td>0.9σ</td><td>3.76 ± 0.66</td><td>2.71 ± 0.18</td><td>5.9</td></tr><tr><td>b</td><td>454.9</td><td>5.6σ</td><td>1.7σ</td><td>4.34 ± 0.64</td><td>2.53 ± 0.13</td><td>7.5</td></tr><tr><td>c</td><td>901.8</td><td>10.3σ</td><td>2.8σ</td><td>6.33 ± 0.58</td><td>2.49 ± 0.08</td><td>12.3</td></tr><tr><td>d</td><td>1014.0</td><td>10.8σ</td><td>7.7σ</td><td>6.66 ± 0.58</td><td>2.20 ± 0.06</td><td>16.1</td></tr><tr><td>e</td><td>430.7</td><td>4.7σ</td><td>2.9σ</td><td>2.84 ± 0.51</td><td>2.35 ± 0.12</td><td>6.7</td></tr><tr><td>f</td><td>648.9</td><td>7.7σ</td><td>4.0σ</td><td>4.60 ± 0.64</td><td>2.33 ± 0.11</td><td>10.0</td></tr><tr><td>g</td><td>1238.5</td><td>13.5σ</td><td>6.0σ</td><td>7.41 ± 0.54</td><td>2.45 ± 0.07</td><td>16.1</td></tr><tr><td>h</td><td>1409.2</td><td>14.5σ</td><td>4.6σ</td><td>8.14 ± 0.54</td><td>2.50 ± 0.06</td><td>17.3</td></tr><tr><td>i</td><td>653.4</td><td>9.0σ</td><td>4.0σ</td><td>6.65 ± 0.71</td><td>2.41 ± 0.09</td><td>11.4</td></tr><tr><td>j</td><td>1229.0</td><td>14.0σ</td><td>6.8σ</td><td>9.07 ± 0.63</td><td>2.39 ± 0.06</td><td>17.7</td></tr><tr><td>k</td><td>1246.4</td><td>13.2σ</td><td>3.6σ</td><td>7.73 ± 0.54</td><td>2.48 ± 0.06</td><td>16.5</td></tr><tr><td>l</td><td>1405.7</td><td>14.1σ</td><td>6.3σ</td><td>7.95 ± 0.54</td><td>2.51 ± 0.06</td><td>16.9</td></tr><tr><td>m</td><td>469.5</td><td>6.8σ</td><td>1.7σ</td><td>5.40 ± 0.73</td><td>2.56 ± 0.13</td><td>8.2</td></tr><tr><td>n</td><td>415.4</td><td>5.1σ</td><td>3.5σ</td><td>3.49 ± 0.62</td><td>2.33 ± 0.13</td><td>7.4</td></tr><tr><td>o</td><td>1259.2</td><td>14.1σ</td><td>5.9σ</td><td>8.23 ± 0.57</td><td>2.39 ± 0.06</td><td>17.7</td></tr><tr><td>p</td><td>996.7</td><td>10.5σ</td><td>4.0σ</td><td>6.29 ± 0.55</td><td>2.36 ± 0.07</td><td>14.7</td></tr></tbody></table>

**Notes.** See Fig. 1b for the definition of the signal regions. Significance values were computed following Li & Ma (1983), assuming a perfect knowledge of the background. $\phi_0$ and $\Gamma$ are the best-fit parameter values of the power-law fit for each region (cf. Eq. (1)). $\sqrt{\Delta TS}$ denotes the square root of the difference in test statistic ($TS = -2 \ln(\mathcal{L})$) between the best-fit power-law model and the null hypothesis (corresponding to $\phi_0 = 0$).

The spectra in the signal regions are remarkably similar to each other, both in terms of the fitted power-law indices (which show no dependence on the separation of the region from the centre, see Fig. 6) and the shape of the spectra as indicated by the extracted flux points (see Fig. 5). The only significant deviation is observed in region ‘d’, where the fitted power-law index deviates from the average of all other regions by ~4σ. We have not been able to identify an issue in the analysis procedure that could explain this deviation, and conclude that it either indicates that the spectrum of the emission in this region is indeed harder, or that it is an unexpectedly large statistical fluctuation. The similarity of the spectra supports our previous finding of a lack of energy-dependent morphology of HESS J1646–458, and motivates the extraction of a combined energy spectrum.

We computed combined flux points for HESS J1646–458 by adding up the flux points from all 16 square regions, where we have used the best-fit flux for each point and region also in cases where an upper limit was previously derived. The result is shown in Fig. 7. Besides statistical uncertainties, the displayed error bars contain a systematic uncertainty arising from the applied background model. The systematic uncertainty is of the same magnitude as the statistical one at the lowest energies considered here, and quickly becomes negligible at higher energies. The combined flux points clearly show that the γ-ray emission of HESS J1646–458 extends to at least several tens of TeV. A comparison with the spectra of the individual signal regions (blue lines in Fig. 7) shows that the total spectrum is not dominated by any single signal region across the entire energy range.

In order to characterise the combined spectrum, we fitted several spectral models to the derived flux points. A simple PL model (cf. Eq. (1)) does not lead to a satisfactory fit ($p = 0.06\%$). The solid orange line in Fig. 7 shows the result for a power law with exponential cut-off (ECPL),

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

for which we obtained (with $E_0 = 1$ TeV kept fixed in the fit) $\phi_0 = (1.00 \pm 0.03) \times 10^{11}$ TeV$^{-1}$ cm$^{-2}$ s$^{-1}$, $\Gamma = 2.30 \pm 0.04$, and $E_c = (44^{+17}_{-11})$ TeV. This corresponds to a total γ-ray luminosity of HESS J1646–458 between the threshold energy of 0.37 TeV and 100 TeV of $L_γ \sim 9 \times 10^{34}$ (d/3.9 kpc)$^2$ erg s$^{-1}$, where d is the assumed distance to the source. We note that, while the ECPL model yields an acceptable fit ($p = 6.3\%$), the high-energy flux points do not provide a clear indication of an exponential cut-off to the spectrum. Thus, while the energy spectrum of HESS J1646–458 clearly extends to several tens of TeV, its maximum energy cannot be determined reliably with the analysis presented here, and may conceivably lie beyond 100 TeV. However, the last flux point should not be regarded as a clear indication of an upward trend in the spectrum, as it may be afflicted by unknown systematic uncertainties in the high-energy response of the system, which is difficult to calibrate.

Assuming the γ-ray emission to be generated in collisions of CR protons with ambient matter, we also fitted a primary proton spectral model (of the same form as defined in Eq. (2)) to the γ-ray flux points, employing the NAIMA software package (Zabalza 2015) for this task. For the parameters of the primary proton spectrum we obtained a normalisation (at $E_0 = 1$ TeV) of $\phi_0^p = (1.28 \pm 0.17) \times 10^{38}$ (d/3.9 kpc)$^2$ (n/1 cm$^{-3}$)$^{-1}$ eV$^{-1}$ (with n the assumed density of the ambient matter), a spectral index of $\Gamma_p = 2.33 \pm 0.06$, and a cut-off energy of $E_c^p = (400^{+250}_{-130})$ TeV; the dashed green line in Fig. 7 displays the corresponding γ-ray spectrum. The 95% confidence level lower limit on the proton spectrum cut-off energy is $E_c^p > 214$ TeV. Extrapolating the proton spectrum down to an energy of 1 GeV, the required energy in primary protons is $W_p \sim 6 \times 10^{51}$ (d/3.9 kpc)$^2$ (n/1 cm$^{-3}$)$^{-1}$ erg.

In a similar manner, now adopting with NAIMA a leptonic framework in which the γ-ray emission is produced through IC scattering of CR electrons, we also fitted a primary electron spectrum to the HESS J1646–458 flux points (again assuming an ECPL model). Besides the cosmic microwave background, we used as target radiation fields an infrared field ($T_{IR} = 26$ K,
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Comparison of signal region spectra. All spectra were divided (at a reference energy of 1 TeV) by a reference power-law spectrum with spectral index $\Gamma = 2.41$, corresponding to the weighted average over all signal regions. Upper limits are at 95% confidence level, and only two upper limits after the last significant (i.e. $>2\sigma$) flux point are shown.</figcaption></figure>

$\rho_{\text{IR}} = 0.74 \, \text{eV cm}^{-3}$) and an optical field ($T_{\text{opt}} = 2400 \, \text{K}$, $\rho_{\text{opt}} = 1.4 \, \text{eV cm}^{-3}$) as predicted by the Popescu et al. (2017) model, as well as an additional field representing diffuse star light from the stellar cluster. For the latter, we assumed $T_{\text{SC}} = 40\,000 \, \text{K}$ (Crowther et al. 2006) and derived<fnref n="3" /> an energy density

<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Spectral index $\Gamma$ for the signal regions a–p, as a function of angular separation between the centre point of each region and the centre point of the total emission (white circle in Fig. 3). The red line and band denote the weighted mean and uncertainty across all regions, $\langle \Gamma \rangle = 2.41 \pm 0.02$, respectively.</figcaption></figure>

of $\rho_{\text{SC}} \sim 30 \, \text{eV cm}^3$ at a distance of 34 pc from the cluster, which approximately corresponds to the distance at which the radial $\gamma$-ray excess profile peaks (cf. Fig. 4). In this case, the best-fit parameters are $\phi_0^e = (4.7 \pm 0.5) \times 10^{35} (d/3.9 \, \text{kpc})^2 \, \text{eV}^{-1}$, $\Gamma_e = 2.97 \pm 0.07$, and $E_c^e = (180_{-70}^{+200}) \, \text{TeV}$, with a 95% confidence level lower limit on $E_c^e$ of 87 TeV – the resulting $\gamma$-ray spectrum is shown by the red, dashed-dotted line in Fig. 7. Electrons down to energies of $\sim 0.4 \, \text{TeV}$ contribute to the $\gamma$-ray emission detected with H.E.S.S. Assuming that the spectrum of primary electrons extends down to 0.1 TeV, we obtained a total required energy in primary electrons of $W_e \sim 7.2 \times 10^{48} (d/3.9 \, \text{kpc})^2 \, \text{erg}$. Dividing the required energy by the (energy-dependent) cooling time due to IC scattering (Hinton & Hofmann 2009) off the different target fields then yields an estimate for the minimum total required power of $L_e > 4.1 \times 10^{35} (d/3.9 \, \text{kpc})^2 \, \text{erg s}^{-1}$. The required power is larger if cooling due to synchrotron emission plays a sizeable role, for example, in a magnetic field with $B = 10 \, \mu\text{G}$, $L_e > 1.7 \times 10^{36} (d/3.9 \, \text{kpc})^2 \, \text{erg s}^{-1}$.

The chosen analysis method and tools in principle also enable a three-dimensional modelling of the $\gamma$-ray emission detected with H.E.S.S., that is, to decompose the emission into several components with separate energy spectra and morphological models. However, owing to the complex structure of the emission, a rather complicated model with multiple distinct components is required to obtain a satisfactory description, and its interpretation depends strongly on the chosen models for each of the components. Considering furthermore the similarity of energy spectra extracted in the signal regions, we do not present such a modelling here.

## 2.3. Analysis of radio line data

Attributing the $\gamma$-ray emission to interactions of accelerated cosmic-ray nuclei requires sufficiently dense target material. The presence of such target material can be estimated using radio observations, in particular of the 21 cm H I emission line – indicating neutral, atomic hydrogen – and of the CO ($J = 1-0$) transition, which is a tracer for dense clouds of molecular hydrogen, H$_2$ (Heyer & Dame 2015). We have therefore used the H I Southern Galactic Plane Survey (SGPS; McClure-Griffiths et al. 2005) and the Mopra Southern Galactic Plane CO Survey

<footnote n="3">³ To derive the energy density, we used $\rho_{\text{SC}} = L/(4\pi r^2 c)$, where L is the total cluster luminosity and r the distance from the cluster. For a wind efficiency $\eta = M v_w / (L/c) \approx 1$ (Vink & Gräfener 2012), $L = 2 (v_w/c)^{-1} L_w \sim 2 \times 10^{41} \, \text{erg s}^{-1}$ for the parameter values assumed in this work. Adding up the luminosities of OB supergiant and Yellow Hypergiant stars listed in Clark et al. (2005) yields $L \sim 6 \times 10^{40} \, \text{erg s}^{-1}$, which is consistent with our estimate. A slight reduction of the IC emission due to the cluster radiation field is expected because of its anisotropy, which we have not taken into account. However, since the contribution of this component is suppressed due to the Klein-Nishina effect, the modification of the total IC emission is negligible.</footnote>
<figure><img src="image_7.png" /></figure>

**Fig. 7.** Combined energy spectrum. The black data points correspond to the entire emission of HESS J1646−458; the solid orange, dashed green, and dashed-dotted red lines show the result of fitting a power law with exponential cut-off (ECPL), a hadronic (pp) model, and a leptonic (IC) model, respectively, to these points. Fitted power-law models for each region a–p are displayed by solid blue lines (with darker shades indicating closer proximity to Westerlund 1), while the dashed blue line denotes their sum. All power-law spectra are plotted up to 100 TeV for better visibility, however, the observed γ-ray excess is not significant up to this energy in any of the sub-regions. *Bottom panel:* ratio to the ECPL model; note that the last flux point (with a ratio to the ECPL model of ~3.7 ± 1.5) is beyond the vertical axis scale.

(Braiding et al. 2018) to investigate the amount of hydrogen gas in the vicinity of Westerlund 1. The analysis was repeated using the CO data from the (lower-resolution) survey by Dame et al. (2001) instead of the Mopra CO data, leading to consistent results.

The radio data analysis is hampered by the uncertainty on the distance to Westerlund 1. We show in Fig. 8 the H I and CO maps for an interval in velocity with respect to the local standard of rest of $v = [-60, -50]$ km s⁻¹, which corresponds to the distance of 3.9 kpc that we adopted for this paper (Kothes & Dougherty 2007). As some measurements indicate smaller distances, maps for two correspondingly chosen alternative velocity intervals, $v = [-48.5, -38.5]$ km s⁻¹ ($d \approx 3.3$ kpc) and $v = [-37, -27]$ km s⁻¹ ($d \approx 2.7$ kpc), are shown in Appendix B.

We find that the gas indicated by the radio observations at a distance of ~3.9 kpc shows no spatial correlation with the γ-ray emission that we observe with H.E.S.S. In fact, both the H I and CO maps indicate a particularly low atomic and molecular gas density in the circular regions B and C, which are bright in γ rays. Using an H I intensity-mass conversion factor of $X_{\text{H I}} = 1.823 \times 10^{18}$ cm⁻²/(K km s⁻¹; Rohlfss & Wilson 2004), we obtain for a circular region with radius 1.1°, centred on Westerlund 1, a total enclosed mass as indicated from H I of $M_{\text{H I, Wd1}} = 1.3 \times 10^5 M_{\odot}$. This translates into an average density

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

**Fig. 8.** Maps showing H I emission (McClure-Griffiths et al. 2005) (top panel) and $^{12}$CO emission (Braiding et al. 2018) (bottom panel) in the Westerlund 1 region. Both maps display the emission for an interval in velocity with respect to the local standard of rest of $v = [-60, -50]$ km s⁻¹, which approximately corresponds to a distance of 3.9 kpc. The position of Westerlund 1 is marked by the white star symbol and the grey, dashed line shows the Galactic plane. The transparent, white circle marker denotes the centre point with respect to which the radial CR density profiles in Fig. 9a have been computed; the dashed white line displays a circle with radius 1° – up to which the profiles have been computed – around this point. The red lines are contour lines of the flux map shown in Fig. 3a. Regions A, B, and C are the same as in Fig. 3.

of $n_{\text{H I, Wd1}} = 3.2 \text{ cm}^{-3}$<fnref n="4" />. Similarly, from the CO data we get<fnref n="5" /> $M_{\text{CO, Wd1}} = 4.3 \times 10^5 M_{\odot}$ and $n_{\text{CO, Wd1}} = 10.5 \text{ cm}^{-3}$, where $n_{\text{CO}}$ is the equivalent density for atomic hydrogen and can thus be directly compared to $n_{\text{H I}}$. We stress, however, that in partic-

<footnote n="4">⁴ Abramowski et al. (2012) derived, for a similar region, a much smaller value of $n_{\text{H I}} = 0.24 \text{ cm}^{-3}$. We attribute this to the usage of an erroneous formula in that paper.</footnote>

<footnote n="5">⁵ Due to the more indirect nature of the estimate, the CO-to-H₂ conversion factor, $X_{\text{CO}}$, is less well constrained than $X_{\text{H I}}$. Here we used $X_{\text{CO}} = 1.5 \times 10^{20}$ cm⁻²/(K km s⁻¹), which Ackermann et al. (2012) indicate as an appropriate value for the galactocentric radius of Westerlund 1, $R \approx 4.7$ kpc for a distance of 3.9 kpc (see their Fig. 25). This value is also within the range of (1.4–2.6) × $10^{20}$ cm⁻²/(K km s⁻¹) recommended by Bolatto et al. (2013). We have neglected the possible contribution of $^{4}$He, which could increase the mass estimate by ~25%.</footnote>
ular the molecular material indicated by the CO observations is not distributed homogeneously inside this region, but rather concentrated in smaller-scale clouds. For instance, in the CO cloud located in the northern part of the H.E.S.S. emission region (cf. Fig. 8), we find – assuming a spherical distribution of the gas – a density of $n_{CO,cloud} \sim 190 \text{ cm}^{-3}$.

Following Aharonian et al. (2019), we also derived a radial profile of the CR density, assuming that the γ-ray emission is produced in interactions of hadronic CRs with the gas indicated by the H I and CO data. Having used the measured γ-ray flux above our threshold energy of 0.37 TeV to compute the profiles, they indicate the density of CRs above an energy about ten times higher, that is, above ~4 TeV. Density profiles for all three considered velocity intervals are shown in Fig. 9, where we have used for the profiles in panel a the same shifted centre point as for the radial excess profiles (cf. Fig. 4), and for the profiles in panel b – for comparison with Aharonian et al. (2019) – the position of Westerlund 1 as centre point. Regarding first Fig. 9a, for the interval corresponding to the nominal distance of $d \approx 3.9 \text{ kpc}$ (blue curve with circle markers), we observe a distinct peak at a radial distance from the centre of ~0.4°, corresponding to the peak observed in the excess profiles (cf. Fig. 4) at about the same distance. When choosing the cluster position as centre point, the peak is smeared out, leading to a plateau at small radii. The profile is then compatible with that derived in Aharonian et al. (2019), and showing a gradual decline of the density moving away from the cluster. We stress, however, that the γ-ray emission is not radially symmetric around Westerlund 1, rendering the interpretation of radial profiles computed with respect to this position difficult. Our profiles for the alternative velocity intervals (i.e. distances; orange and green curves with triangle markers in Fig. 9) exhibit approximately the same shape as the profile for the nominal distance, but with a less distinct peak, and lower overall density. This corresponds to the higher density of hydrogen gas observed for these intervals (cf. Figs. B.1 and B.2). Finally, we note that the shapes of the obtained radial density profiles should be interpreted with care, owing to systematic uncertainties associated with the determination of the gas distributions. For instance, it is conceivable that a considerable amount of the CO molecules have been photodissociated by the intense ultraviolet radiation from Westerlund 1, implying that the CO line emission would no longer be an accurate tracer of molecular hydrogen gas (e.g. Wolfire et al. 2010).

## 3. Discussion

In this section, we will discuss in turn the various possible sites for the acceleration of CRs that could be responsible for the observed γ-ray emission. Evidently, we can include in this discussion only those objects that have been identified through observations at other wavelengths. While we consider the association of the γ-ray emission with one of these objects most likely, so far undiscovered objects (e.g. pulsars or SNRs) could also contribute to it.

### 3.1. Objects not connected to Westerlund 1

#### 3.1.1. 4U 1642–45

The LMXB 4U 1642–45 lies close to the peak in γ-ray emission observed in region A (cf. Fig. 3). However, despite this positional coincidence, we consider an association of 4U 1642–45 with HESS J1646–458 – or even the emission observed just in region A – highly unlikely: LMXBs are not known emit-

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

**Fig. 9.** Cosmic-ray density profiles above 4 TeV for different velocity intervals. The profiles in the upper plot have been computed with respect to the centre point shown by the circle marker in Fig. 8 (the same as for the excess profiles displayed in Fig. 4), whereas the profiles in the lower plot have been computed with respect to the position of Westerlund 1. The error bars denote the (statistical) uncertainty of the measured γ-ray flux only, and in particular do not reflect systematic uncertainties related to the gas distribution. The CR density at Earth, indicated for comparison, was computed from the all-particle CR spectrum as modelled in Breuhaus et al. (2022) (model ‘A’). The data points taken from Aharonian et al. (2019), shown in the lower plot, have been derived using a velocity interval of $v = [-60, -50] \text{ km s}^{-1}$, and should thus be compared to the blue curve with circle markers. The red, solid line, from the same paper, shows a profile corresponding to a $1/r$-distribution of the CRs, where $r$ is the distance to the cluster. Both show the density above a slightly higher threshold energy of 10 TeV.

ters of γ rays in the TeV energy range, in fact, Kantzas et al. (2022) have recently shown that even the upcoming Cherenkov Telescope Array (CTA) will be able to detect LMXB outbursts only under favourable circumstances. Moreover, the observed emission is spatially extended, which would be unexpected if it is produced inside a transient jet. Lastly, Abramowski et al. (2012) reported no temporal variations of the observed γ-ray flux, again disfavouring an association of HESS J1646–458 with 4U 1642–45. While the deviation in spectral index observed in signal region ‘d’ (cf. Sect. 2.2.2), which partly coincides with the circular region A, is intriguing, the lack of a plausible association leads us to conclude that it is likely a statistical fluctuation, or due to an unidentified, hard-spectrum source contributing to
the emission in region d (but not to the entire region A, as the signal regions c, g, and h, which also overlap with region A, do not show deviating spectral indices).

### 3.1.2. PSR J1648-4611 and PSR J1650-4601

The two energetic pulsars PSR J1648-4611 (with spin-down power $\dot{E} = 2.1 \times 10^{35}$ erg s$^{-1}$) and PSR J1650-4601 ($\dot{E} = 2.9 \times 10^{35}$ erg s$^{-1}$; Manchester et al. 2005) are located within region C (cf. Fig. 3), where we observe enhanced γ-ray emission. Because pulsar wind nebulae (PWNe) represent a large fraction of known Galactic γ-ray sources (e.g. Abdalla et al. 2018c), it seems likely that high-energy electrons provided by one (or both) of the two pulsars contribute to the γ-ray emission observed in region C. This remains true even though neither of the pulsar locations fully coincides with the peak of the emission, as it is not uncommon that γ-ray PWNe are observed offset from their respective pulsars. The detection of diffuse, hard X-ray emission from the vicinity of PSR J1648-4611 by Sakai et al. (2013) adds further support for this scenario. However, while the 4FGL-DR2 catalogue (Abdollahi et al. 2020; Ballet et al. 2020) contains sources associated with PSR J1648-4611 and PSR J1650-4601, their γ-ray emission is detected as pulsed and exhibits a very steep spectrum above 10 GeV, implying that these sources are directly connected to the pulsars rather than their putative PWNe.

On the other hand, viewing the entire emission of HESS J1646-458 as resulting from one of the two pulsars would imply an unusually large PWN. For instance, at the distance of PSR J1648-4611 of 5.7 kpc (Kramer et al. 2003), the emission region spans ~200 pc, which is twice the size of the largest known γ-ray PWN, HESS J1825-137 (Abdalla et al. 2019). For such an extended nebula, one would expect a considerable loss of energy due to synchrotron cooling of the electrons when they propagate towards the edges of the nebula, leading to softer γ-ray spectra in these regions (or, equivalently, energy-dependent morphology). Because we observe very similar energy spectra across the entire source, and no energy-dependent morphology, we conclude that a PWN powered by PSR J1648-4611 or PSR J1650-4601 cannot explain the entire γ-ray emission.

## 3.2. Distinct acceleration sites within Westerlund 1

Having established that other known objects in the region cannot explain the bulk of γ-ray emission from HESS J1646-458, we assume next that the CRs producing the emission are accelerated at one or multiple sites within the stellar cluster Westerlund 1. Various scenarios can be considered, and will be discussed in this section.

### 3.2.1. CXOU J164710.2-455216

At first sight, the magnetar CXOU J164710.2-455216 – the only known stellar remnant inside the cluster – may be suspected to power a γ-ray PWN that could potentially be associated to HESS J1646-458. However, its measured period and period derivative ($P = 10.6$ s, $\dot{P} = 9.2 \times 10^{-13}$ s s$^{-1}$; Israel et al. 2007) imply a rotational spin-down power of only $\dot{E} = 3 \times 10^{31}$ erg s$^{-1}$, which is orders of magnitude lower than for any of the pulsars associated with PWN detected at TeV energies with H.E.S.S. (Abdalla et al. 2018d). Even though the measured X-ray luminosity of $L_X \approx 3 \times 10^{33}$ erg s$^{-1}$ (Muno et al. 2006b) exceeds the rotational spin-down power, implying another source of energy (presumably connected to the high magnetic field of the magnetar), it still appears very unlikely that CXOU J164710.2-455216 would be able to sustain the observed γ-ray emission. Additionally, as is the case for PSR J1648-4611 and PSR J1650-4601, the spatial extent of HESS J1646-458 and the lack of energy-dependent morphology disfavour an association of the γ-ray emission with CXOU J164710.2-455216.

### 3.2.2. Supernova remnants

The existence of CXOU J164710.2-455216 implies that at least one supernova (SN) explosion took place within Westerlund 1. However, given the abundance of massive stars in Westerlund 1, and its age of several Myr, it seems certain that many more SNe have occurred already. Indeed, Muno et al. (2006a) have argued that the number of SNe during the last ~1 Myr could be as high as 80–150, attributing the lack of identified SNRs to a cavity in the interstellar medium (ISM), excavated by the stellar cluster. Without detailed knowledge about the progenitor mass of CXOU J164710.2-455216, or the SN rate in Westerlund 1, it is not straightforward to estimate the energy output that SNRs in the stellar cluster may provide. Nevertheless, assuming a canonical value for the kinetic energy released per SN of $10^{51}$ erg – possibly more in the case of CXOU J164710.2-455216, if its progenitor was indeed as heavy as $40 M_{\odot}$ (Muno et al. 2006b) – it seems plausible that one or several SNRs in Westerlund 1 could be responsible for the γ-ray emission in terms of the required energetics. For instance, in a hadronic scenario, assuming a density of the ambient matter of $5 m_{\text{H}} \text{ cm}^{-3}$ gives a required energy in protons of $1.2 \times 10^{51}$ erg (cf. Sect. 2.2.2). It would take about 10 SNRs to reach this energy if the conversion efficiency from kinetic energy into CRs is ~10%.

### 3.2.3. SN-wind and wind-wind interactions

Another, related, possibility for the acceleration of CRs inside Westerlund 1 are interactions of SN shocks with winds of massive stars in the cluster, or interactions between the winds of several stars. Indeed, for SN shocks interacting with fast stellar winds, the efficiency for CR acceleration can be as high as 30% (Bykov et al. 2020), and CR acceleration up to $\ge 40$ PeV has been conjectured (Bykov et al. 2015). Colliding winds of massive stars are also known to produce non-thermal emission (Eichler & Usov 1993; Reimer et al. 2006), and several searches for γ-ray emission from colliding wind binaries have been performed in the past (e.g. Werner et al. 2013; Pshirkov 2016; Martí-Devesa & Reimer 2021). A well-known example is provided by the colliding wind binary η Car, which has indeed been detected up to ~100 GeV with the Fermi-LAT (Tavani et al. 2009; Abdo et al. 2010), and even up to ~400 GeV with H.E.S.S. (Abdalla et al. 2020). These considerations strengthen the conclusion that γ-ray emission at the level observed with H.E.S.S. can in principle be produced by CRs accelerated at shock fronts within Westerlund 1.

However, for the scenario of a central CR source, it is important to also take into account the propagation of the CRs into the region where we observe γ-ray emission with H.E.S.S.. In this case, the non-observation of a peak in the γ-ray emission at the position of Westerlund 1 and the large extent of HESS J1646-458 essentially rule out the leptonic scenario. In a hadronic scenario, the complex morphology of HESS J1646-458 may in principle be attributable to the distribution of target material – although a clear correlation of the γ-ray emission with gas clouds as indicated by H I and CO
observations is lacking, and the inferred CR density does not peak at the centre for any of the considered distances to the source (cf. Fig. 9a), as would be expected for a steady CR injection there (see Aharonian et al. 2019, but also Bhadra et al. 2022). Possible options to alleviate this problem could be the presence of ‘dark’ gas that is not traced by HI or CO (e.g. Wolfire et al. 2010), or the assumption that the CRs were provided by a recent impulsive event (e.g. the SN explosion of the CXOU J164710.2–455216 progenitor star and/or other recent SNe), rather than being injected quasi-continuously over the lifetime of the cluster. Another relevant constraint comes from the maximum energy of the observed γ rays, which implies the presence of CRs with energies of several hundred TeV throughout the emission region. If the acceleration sites are located exclusively within the compact cluster, particles must pass through the wind zone. Taking reasonable limits on the diffusion properties, adiabatic losses in the radial wind are unavoidable. CRs that were injected at the cluster and have propagated to a distance R within the wind region would therefore need to be produced with maximum energy ($R/R_{Wd1}$)$^{2/3}$ times larger (Longair 1992), where $R_{Wd1} \sim 1$ pc is the radius of Westerlund 1. For the nominal distance of 3.9 kpc, we observe a peak in the γ-ray emission at R $\sim$ 34 pc, implying a need of multi-PeV CRs within Westerlund 1. The H.E.S.S. observations thus provide a valuable constraint for theoretical models of particle acceleration processes within a compact stellar cluster (Bykov et al. 2020).

### 3.3. Acceleration by Westerlund 1 as a whole

Finally, there is the possibility that effects due to the entire stellar cluster provide the means for efficient CR acceleration. In particular, we consider in the following two scenarios related to the collective cluster wind, in which the CR acceleration predominantly takes place outside of the actual stellar cluster itself.

#### 3.3.1. Turbulence in a superbubble

Due to the powerful collective cluster wind, as well as many SN explosions, massive young stellar clusters are thought to create ‘superbubbles’, that is, large cavities in the ISM, extending much beyond the boundaries of the cluster itself. In the shocked medium inside the bubble, strong magnetohydrodynamic (MHD) turbulences provide the conditions for particle acceleration via the second-order Fermi mechanism (e.g. Bykov 2014; Vieu et al. 2022). With a maximum proton energy of 200 TeV and a source extent of $\sim$100 pc, the Hillas criterion (Hillas 1984) implies for an average turbulent fluid velocity $u = 100 \text{ km s}^{-1}$ a minimum magnetic field strength of $\sim$13 µG in this scenario. While the acceleration time scales are much longer compared to the case of acceleration at shocks inside the cluster, the process could – under favourable circumstances – generate CRs with up to PeV energies during several Myr, the typical life time of young clusters (Bykov et al. 2020). For an adiabatically expanding wind, and a density of the ambient ISM $\rho_0$, the radius of the superbubble is given by $R_{SB} = 0.76 (L_w/\rho_0)^{1/5} \tau^{3/5}$ (Weaver et al. 1977; Koo & McKee 1992), or $R_{SB} \sim 256 (\rho_0/1 \text{ m}_H \text{ cm}^{-3})^{-1/5}$ pc for our assumptions (cf. Table 1). Adopting, for example, $\rho_0 = 5 \text{ m}_H \text{ cm}^{-3}$, we obtain $R_{SB} \sim 185$ pc – a value more than two orders of magnitude larger than the half-mass radius of the cluster. A superbubble with a size of this order around Westerlund 1 has not yet been revealed at other wavelengths. Because its dimensions also exceed the extent of the γ-ray emission detected with H.E.S.S., an asso-

ciation of this emission with the entire superbubble seems disfavoured. However, the assumption of a homogeneous medium is an oversimplification, and bubbles in a structured and possibly clumpy medium may have significantly different cooling rates and dynamics (e.g. Chu 2008). Moreover, more detailed models of superbubble evolution indicate that the simple estimate for their radius given above often over-predicts their true size (Yadav et al. 2017; Vieu et al. 2022). If this is also the case for Westerlund 1, smaller-scale structures as for example the bubble-like feature ‘B3’ in HI data reported by Kothes & Dougherty (2007) could be associated with the cluster. Nevertheless, even in this case a connection between the superbubble and the γ-ray emission is not obvious (but conceivable, considering also that the γ-ray emission need not be uniform from the superbubble volume).

#### 3.3.2. Termination shock of collective cluster wind

Another possible site for the acceleration of CRs due to Westerlund 1, but outside the stellar cluster itself, is the termination shock of the collective cluster wind. The termination shock forms where the pressure of the outgoing wind equals that of the ISM. Recently, Morlino et al. (2021) have proposed that termination shocks of collective stellar cluster winds may be efficient sites of particle acceleration, and demonstrated that CRs with PeV energies could be produced in powerful clusters like Westerlund 1. Considering again as above an adiabatic expansion, the radius of the termination shock is given by $R_{TS} = 0.92 (L_w/\rho_0)^{3/10} v_w^{-1/2} \tau^{2/5}$ (Koo & McKee 1992), or $R_{TS} \sim 51 (\rho_0/1 \text{ m}_H \text{ cm}^{-3})^{-3/10}$ pc with our adopted parameter values. Inserting $\rho_0 = 5 \text{ m}_H \text{ cm}^{-3}$ yields a radius of $R_{TS} \sim 32$ pc. Notably, this value coincides well with the radial distance of $\sim$34 pc at which we observe a peak in the γ-ray excess profiles (cf. Sect. 2.2.1 and Fig. 4). The scenario of particle acceleration at the cluster wind termination shock therefore provides a natural explanation for the shell-like structure exhibited by the γ-ray emission detected with H.E.S.S., and can furthermore reproduce its radial extent under reasonable assumptions. The apparent asymmetry of the shell with respect to the position of Westerlund 1 could be caused, for example, by a gradient in the density of the surrounding ISM, or by a SN that occurred within the cluster towards the direction of the asymmetry.

Adopting the hadronic scenario, we find that the termination shock model is also viable in terms of the required energetics: with $L_w \sim 10^{39}$ erg s$^{-1}$ and $\tau \sim 4$ Myr, the total available energy is $\sim 1.3 \times 10^{53}$ erg, which in principle suffices to explain the required energy in CR protons, $W_p \sim 6 \times 10^{51} (n/1 \text{ cm}^{-3})^{-1}$ erg – although, since the cooling time for protons exceeds the cluster lifetime, the acceleration process would need to be rather efficient, or the target density high. Furthermore, the energetics argument presupposes that the CRs can be confined within the γ-ray emission region over a significant fraction of the full cluster lifetime, which is not straightforward. For instance, adopting Bohm diffusion, the diffusion length for protons is $L \sim \sqrt{6Dt} \sim 81 (E/100 \text{ TeV})^{1/2} (B/10 \text{ µG})^{-1/2} (t/1 \text{ Myr})^{1/2}$ pc (Chandrasekhar 1943), where we have neglected projection effects. Hence, even for slow diffusion, to confine protons with energy 200 TeV – our lower limit for the cut-off energy of the primary proton spectrum in a hadronic scenario – within a region of radius $\sim$50 pc for only 1 Myr already requires a rather large magnetic field strength of $B \sim 50$ µG. Nevertheless, taking into account, for example, an additional smearing due to the transformation from protons to γ rays, the observed γ-ray morphology
can be reproduced with a not-too-extreme assumption on the magnetic field, provided that the protons do not diffuse too fast. However, as already mentioned in Sect. 3.2.3, the hadronic scenario is challenged further by the absence of a correlation of the γ-ray emission with the gas observed in the region.

Interestingly, the finding that the expected location of the termination shock coincides with the shell-like structure of the measured γ-ray emission also renders an interpretation within the leptonic scenario possible. This is because the geometry of the acceleration site naturally explains the relatively large extent of the emission region and its rather complex structure, which are otherwise not easy to accommodate in a leptonic scenario. The scenario is also feasible energetically; even in the presence of a 10 µG magnetic field, the required power of $L_e \sim 1.7 \times 10^{36}$ erg s⁻¹ (cf. Sect. 2.2.2) is easily provided by the cluster wind if the acceleration efficiency for electrons is of order 0.1%. However, as the accelerated electrons emit synchrotron radiation, the scenario is subject to constraints from observations at the corresponding wavelengths. For example, from measurements by the Planck satellite at a frequency of 30 GHz<fnref n="6" />, at which the radiation is dominated by synchrotron emission (Akrami et al. 2020), we infer an average intensity within 1° around Westerlund 1 of 0.55 MJy sr⁻¹. For a magnetic field of 10 µG, electrons with energies around 0.01 TeV emit synchrotron radiation at 30 GHz. Assuming that the electron spectrum extends down to these energies, the predicted intensity at 30 GHz is ~0.3 MJy sr⁻¹. Considering furthermore that only part of the emission detected with Planck originates from the vicinity of Westerlund 1, we conclude that the Planck measurements imply either a magnetic field strength smaller than 10 µG or a low-energy cut-off of the primary electron spectrum. Finally, it is worth noting in this regard that for another superbubble detected at TeV energies, 30 Dor C in the Large Magellanic Cloud, a synchrotron shell has been detected using X-ray measurements, and a leptonic scenario was found to be favoured to explain the TeV γ-ray emission (Kavanagh et al. 2019, and references therein).

### 3.4. Escape of particles from the emission region

So far we have assumed that particles accelerated in or around Westerlund 1 are confined over the lifetime of the cluster. If a significant fraction of accelerated particles can escape there are several important consequences: (i) γ-ray emission would be expected outside the system, in the case of hadronic CRs in particular in molecular clouds; (ii) in the (most likely) case of energy-dependent escape, the inferred spectrum of particles within the system would be softened with respect to the injection spectrum by the energy-dependent escape probability; (iii) the total energy requirements would be increased. There is little evidence for (i) in the maps shown in Fig. 3 – with the possible exception of the emission region east of region ‘C’, which, however, does not coincide with a molecular cloud at the nominal distance of Westerlund 1 (cf. Fig. 8). The inferred injection indices for electrons and protons of $\Gamma_e \sim 3.0^7$ and $\Gamma_p \sim 2.3$ seem broadly consistent with acceleration theory (e.g. Bell 2013), so there seems to be no indication for (ii), although a modest variation of the spectrum is tolerable within the precision of the

H.E.S.S. measurement. Finally, as the total energy requirement is already a challenge in most acceleration scenarios under the assumption of confinement, there also seems to be little room for (iii). Thus, while not entirely inconceivable, we at least find no indications for particles escaping from the emission region.

In the absence of evidence for particle escape, we need to consider if confinement over the cluster lifetime is reasonable or not. As already discussed in Sect. 3.3.2, this is not straightforward in the case of CR protons, which in a disordered magnetic field would diffuse much too quickly. A possible way to circumvent this problem would be a magnetic field topology in which field lines are preferentially in the plane orthogonal to the radial direction, which can substantially inhibit the radial diffusion.

## 4. Conclusion

We have presented a detailed analysis of HESS J1646-458, a very-high-energy γ-ray source positionally coincident with the young massive stellar cluster Westerlund 1. HESS J1646-458 is largely extended (~2° in diameter), and exhibits a complex, shell-like structure, with Westerlund 1 close to its centre. We found no indications for energy-dependent morphology. The energy spectrum of HESS J1646-458 extends to at least several tens of TeV, with a spectral index of ~-2.3, and a gradual steepening above ~10 TeV. Energy spectra extracted within 16 signal regions across the source region are very similar to each other, reinforcing the observation that the morphology of HESS J1646-458 does not vary with γ-ray energy.

In a hadronic scenario with CR protons producing the γ rays, the observed γ-ray spectrum implies proton energies in excess of several hundred TeV. However, our analysis of the H I and CO emission around Westerlund 1 indicates no clear correlation between hydrogen gas clouds and the γ-ray emission, as would be expected to some degree within such a scenario. Nevertheless, in particular due to uncertainties related to the distribution of target gas, a hadronic scenario remains viable in principle. On the other hand, the lack of significant energy-dependent morphology of the γ-ray emission represents a challenge for an interpretation within an IC-dominated, leptonic scenario.

Investigating the possible physical counterparts to HESS J1646-458, we found that – while the energetic pulsars PSR J1648-4611 and PSR J1650-4601 may be contributing to the emission in their immediate surroundings – no other known object besides Westerlund 1 can be made responsible for the bulk of the γ-ray emission. Particle acceleration due to the cluster may occur at various possible sites: at wind-wind or SN-wind interaction shocks within the cluster, at turbulences in the superbubble excavated by the collective cluster wind, or at the termination shock of the cluster wind. Models in which the acceleration takes place within the cluster generally need to overcome the problem of transporting the accelerated CRs into the larger region from which we observe the γ-ray emission without too severe energy losses, and explain the fact that the γ-ray emission does not peak towards the cluster position – in particular the latter argument rules out a leptonic scenario with continuous injection for this case. Attributing the CR acceleration to the superbubble as a whole seems disfavoured because HESS J1646-458 – although largely extended – is still significantly smaller than the expected size of the superbubble, which has, furthermore, so far eluded its detection at other wavelengths. Therefore, we deem most attractive the scenario of CR acceleration at the cluster wind termination shock, because it provides a natural explanation for the shell-like morphology of HESS J1646-458, and the wind is powerful enough to sustain

<footnote n="6">⁶ We have used the full-sky frequency map at 30 GHz available through the Planck Legacy Archive, http://pla.esac.esa.int/pla</footnote>

<footnote n="7">⁷ The electron spectral index derived in the NAIMA fit corresponds to the present-time population of electrons, whose energy spectrum is steepened with respect to the injection spectrum due to cooling.</footnote>
the γ-ray emission. Based on the available data, however, we are not able to firmly identify the acceleration mechanism at work.

Our results further support massive stellar clusters as CR accelerators, and motivate to investigate Westerlund 1 and other representatives of this class of objects more deeply in the future. In particular, we encourage a deep and broad coverage in X-ray observations of the region around Westerlund 1, which may enable the identification of the cluster wind termination shock. On the other hand, a more accurate measurement of the γ-ray emission in this region will be provided by the upcoming Cherenkov Telescope Array (Acharya et al. 2018), which is designed to be an order of magnitude more sensitive than the H.E.S.S. experiment. Exploiting the data from this and other observatories will be crucial in understanding the contribution of massive stellar clusters to the sea of CRs in the Milky Way.

**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'Energie atomique et aux Énergies alternatives (CEA), the UK Science and Technology Facilities Council (STFC), 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 GAMMAPY (https://gammapy.org; Deil et al. 2017, 2020), ASTROPY (https://www.astropy.org; Astropy Collaboration 2013, 2018), MATPLOTLIB (https:// matplotlib.org; Hunter 2007), and NAIMA (https://naima.readthedocs.io; Zabalza 2015) software packages.

## References

Abdalla, H., Abramowski, A., Aharonian, F., et al. 2018a, A&A, 612, A11

Abdalla, H., Aharonian, F., Ait Benkhali, F., et al. 2018b, A&A, 620, A66

Abdalla, H., Abramowski, A., Aharonian, F., et al. 2018c, A&A, 612, A1

Abdalla, H., Abramowski, A., Aharonian, F., et al. 2018d, A&A, 612, A2

Abdalla, H., Aharonian, F., Ait Benkhali, F., et al. 2019, A&A, 621, A116

Abdalla, H., Adam, R., Aharonian, F., et al. 2020, A&A, 635, A167

Abdo, A. A., Allen, B., Berley, D., et al. 2007, ApJ, 658, L33

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 723, 649

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33

Abeysekara, A. U., Archer, A., Aune, T., et al. 2018, ApJ, 861, 134

Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2021, Nat. Astron., 5, 465

Abramowski, A., Acero, F., Aharonian, F., et al. 2011, A&A, 525, A46

Abramowski, A., Acero, F., Aharonian, F., et al. 2012, A&A, 537, A114

Abramowski, A., Aharonian, F., Ait Benkhali, F., et al. 2014a, MNRAS, 439, 2828

Abramowski, A., Aharonian, F., Ait Benkhali, F., et al. 2014b, ApJ, 794, L1

Abramowski, A., Aharonian, F., Ait Benkhali, F., et al. 2014c, Phys. Rev. D, 90, 122007

Abramowski, A., Aharonian, F., Ait Benkhali, F., et al. 2015, Science, 347, 406

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

Ackermann, M., Ajello, M., Allafort, A., et al. 2011, Science, 334, 1103

Ackermann, M., Ajello, M., Atwood, W. B., et al. 2012, ApJ, 750, 3

Aghakhanloo, M., Murphy, J. W., Smith, N., et al. 2020, MNRAS, 492, 2497

Aghakhanloo, M., Murphy, J. W., Smith, N., et al. 2021, Res. Notes Am. Astron. Soc., 5, 14

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

Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007, A&A, 467, 1075

Aharonian, F., Yang, R., & de Oña Wilhelmi, E. 2019, Nat. Astron., 3, 561

Akrami, Y., Argüeso, F., Ashdown, M., et al. 2020, A&A, 641, A2

Amenomori, M., Bao, Y. W., Bi, X. J., et al. 2021, Phys. Rev. Lett., 126, 141101

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

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

Ballet, J., Burnett, T. H., Digel, S. W., & Lott, B. 2020, ArXiv e-prints [arXiv:2005.11208]

Beasor, E. R., Davies, B., Smith, N., Gehrz, R. D., & Figer, D. F. 2021, ApJ, 912, 16

Belczynski, K., & Taam, R. E. 2008, ApJ, 685, 400

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

Bhadra, S., Gupta, S., Nath, B. B., & Sharma, P. 2022, MNRAS, 510, 5579

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207

Braiding, C., Wong, G. F., Maxted, N. I., et al. 2018, PASA, 35, e029

Brandner, W., Clark, J. S., Stolte, A., et al. 2008, A&A, 478, 137

Breuhaus, M., Hinton, J. A., Joshi, V., Reville, B., & Schoorlemmer, H. 2022, A&A, 661, A72

Bykov, A. M. 2014, A&ARv, 22, 77

Bykov, A. M., Ellison, D. C., Gladilin, P. E., & Osipov, S. M. 2015, MNRAS, 453, 113

Bykov, A. M., Marcowith, A., Amato, E., et al. 2020, Space Sci. Rev., 216, 42

Cao, Z., Aharonian, F. A., An, Q., et al. 2021, Nature, 594, 33

Cesarsky, C. J., & Montmerle, T. 1983, Space Sci. Rev., 36, 173

Chandrasekhar, S. 1943, Rev. Mod. Phys., 15, 1

Chu, Y.-H. 2008, in Proc. Int. Astron. Union, eds. F. Bresolin, P. A. Crowther, & J. Puls, 250, 341

Clark, J. S., Negueruela, I., Crowther, P., & Goodwin, S. P. 2005, A&A, 434, 949

Clark, J. S., Muno, M. P., Negueruela, I., et al. 2008, A&A, 477, 147

Clark, J. S., Najarro, F., Negueruela, I., et al. 2019, A&A, 623, A83

Crowther, P. A., Hadfield, L. J., Clark, J. S., Negueruela, I., & Vacca, W. D. 2006, MNRAS, 372, 1407

Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792

Davies, B., & Beasor, E. R. 2019, MNRAS, 486, L10

de Neurois, M., & Rolland, L. 2009, Astropart. Phys., 32, 231

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

Deil, C., Wood, M., Hassan, T., et al. 2018, https://doi.org/10.5281/zenodo.1409831

Deil, C., Donath, A., Terrier, R., et al. 2020, https://doi.org/10.5281/zenodo.4701492

Eichler, D., & Usov, V. 1993, ApJ, 402, 271

Forman, W., Jones, C., Cominsky, L., et al. 1978, ApJS, 38, 357

Fujii, M. S., & Portegies Zwart, S. 2016, ApJ, 817, 4

Gupta, S., Nath, B. B., Sharma, P., & Eichler, D. 2018, MNRAS, 473, 1537

Heyer, M., & Dame, T. M. 2015, ARA&A, 53, 583

Hillas, A. M. 1984, ARA&A, 22, 425

Hinton, J. A., & Hofmann, W. 2009, ARA&A, 47, 523

Holler, M., Berge, D., van Eldik, C., et al. 2015, Proc. 34th Int. Cosmic Ray Conf. (ICRC2015), 847

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

Israel, G. L., Campana, S., Dall'Osso, S., et al. 2007, ApJ, 664, 448

Kantzas, D., Markoff, S., Lucchini, M., et al. 2022, MNRAS, 510, 5187

Kavanagh, P. J., Norci, L., & Meurs, E. J. A. 2011, New Astron., 16, 461

Kavanagh, P. J., Vink, J., Sasaki, M., et al. 2019, A&A, 621, A138

Kissmann, R. 2014, Astropart. Phys., 55, 37

Kissmann, R., Werner, M., Reimer, O., & Strong, A. W. 2015, Astropart. Phys., 70, 39

Kissmann, R., Niederwanger, F., Reimer, O., & Strong, A. W. 2017, AIP Conf. Proc., 1792, 070011

Knödlseder, J., Mayer, M., Deil, C., et al. 2016, A&A, 593, A1

Koo, B.-C., & McKee, C. F. 1992, ApJ, 388, 93

Kothes, R., & Dougherty, S. M. 2007, A&A, 468, 993

Kramer, M., Bell, J. F., Manchester, R. N., et al. 2003, MNRAS, 342, 1299

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

Longair, M. S. 1992, High Energy Astrophysics. Vol. 1: Particles, Photons and their Detection (Cambridge: Cambridge University Press)

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

Martí-Devesa, G., & Reimer, O. 2021, A&A, 654, A44

McClure-Griffiths, N. M., Dickey, J. M., Gaensler, B. M., et al. 2005, ApJS, 158, 178

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

Morlino, G., Blasi, P., Peretti, E., & Cristofari, P. 2021, MNRAS, 504, 6096

Muno, M. P., Law, C., Clark, J. S., et al. 2006a, ApJ, 650, 203

Muno, M. P., Clark, J. S., Crowther, P. A., et al. 2006b, ApJ, 636, L41

Negueruela, I., Alfaro, E. J., Dorda, R., et al. 2022, A&A, 664, A146
Ohm, S., van Eldik, C., & Egberts, K. 2009, Astropart. Phys., 31, 383

Ohm, S., Hinton, J. A., & White, R. 2013, MNRAS, 434, 2289

Parizot, E., Marcowith, A., van der Swaluw, E., Bykov, A. M., & Tatischeff, V. 2004, A&A, 424, 747

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

Piron, F., Djannati-Ataï, A., Punch, M., et al. 2001, A&A, 374, 895

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

Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431

Pshirkov, M. S. 2016, MNRAS, 457, L99

Rate, G., Crowther, P. A., & Parker, R. J. 2020, MNRAS, 495, 1209

Reimer, A., Pohl, M., & Reimer, O. 2006, ApJ, 644, 1118

Rohlfs, K., & Wilson, T. L. 2004, Tools of Radio Astronomy, 4th edn. (Berlin, Heidelberg: Springer-Verlag)

Sakai, M., Matsumoto, H., Haba, Y., Kanou, Y., & Miyamoto, Y. 2013, PASJ, 65, 64

Specovius, A. 2021, PhD Thesis, Friedrich-Alexander-Universität Erlangen-
Nürnberg

Sun, X.-N., Yang, R.-Z., Liang, Y.-F., et al. 2020a, A&A, 639, A80

Sun, X.-N., Yang, R.-Z., & Wang, X.-Y. 2020b, MNRAS, 494, 3405

Tavani, M., Sabatini, S., Pian, E., et al. 2009, ApJ, 698, L142

Vieu, T., Gabici, S., Tatischeff, V., & Ravikularaman, S. 2022, MNRAS, 512, 1275

Vink, J. S., & Gräfener, G. 2012, ApJ, 751, L34

Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377

Werner, M., Reimer, O., Reimer, A., & Egberts, K. 2013, A&A, 555, A102

Westerlund, B. 1961, PASP, 73, 51

Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010, ApJ, 716, 1191

Yadav, N., Mukherjee, D., Sharma, P., & Nath, B. B. 2017, MNRAS, 465, 1720

Yang, R.-Z., & Aharonian, F. 2017, A&A, 600, A107

Yang, R.-Z., & Wang, Y. 2020, A&A, 640, A60

Yang, R., de Oña Wilhelmi, E., & Aharonian, F. 2018, A&A, 611, A77

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

Zorn, J. 2019, PhD Thesis, Ruprecht-Karls-Universität Heidelberg

¹<fnref n="4" /> Sorbonne Université, Université Paris Diderot, Sorbonne Paris Cité, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Énergies, LPNHE, 4 Place Jussieu, 75252 Paris, France

¹<fnref n="5" /> Instytut Fizyki Jądrowej PAN, ul. Radzikowskiego 152, 31-342 Kraków, Poland

¹<fnref n="6" /> School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue, Braamfontein, Johannesburg 2050, South Africa

¹<fnref n="7" /> Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université de Paris, 92190 Meudon, France

¹<fnref n="8" /> Aix Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

¹<fnref n="9" /> School of Physical Sciences, University of Adelaide, Adelaide 5005, Australia

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

²<fnref n="1" /> Laboratoire Univers et Particules de Montpellier, Université Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, 34095 Montpellier Cedex 5, France

²<fnref n="2" /> Landessternwarte, Universität Heidelberg, Königstuhl 69117, Heidelberg, Germany

²<fnref n="3" /> Université Bordeaux, CNRS, LP2I Bordeaux, UMR 5797, 33170 Gradignan, France

²<fnref n="4" /> University of Oxford, Department of Physics, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK

²<fnref n="5" /> Institut für Astro- und Teilchenphysik, Leopold-Franzens-Universität Innsbruck, 6020 Innsbruck, Austria

²<fnref n="6" /> Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149, 22761 Hamburg, Germany

²<fnref n="7" /> Obserwatorium Astronomiczne, Uniwersytet Jagielloński, ul. Orla 171, 30-244 Kraków, Poland

²<fnref n="8" /> Institute of Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland

²<fnref n="9" /> Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18, 00-716 Warsaw, Poland

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

³<fnref n="1" /> Yerevan Physics Institute, 2 Alikhanian Brothers St., 375036 Yerevan, Armenia

³<fnref n="2" /> Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Strasse 24/25, 14476 Potsdam, Germany

³<fnref n="3" /> Department of Physics, Konan University, 8-9-1 Okamoto, Higashinada, Kobe, Hyogo 658-8501, Japan

³<fnref n="4" /> Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS), The University of Tokyo, 5-1-5 Kashiwa-no-Ha, Kashiwa, Chiba 277-8583, Japan

³<fnref n="5" /> RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan

³<fnref n="6" /> Department of Physics, University of the Free State, PO Box 339, Bloemfontein 9300, South Africa

³<fnref n="7" /> School of Physics, University of New South Wales, Sydney, NSW 2052, Australia

³<fnref n="8" /> Armagh Observatory and Planetarium, Armagh BT61 7HT, UK

³<fnref n="9" /> Western Sydney University, Locked Bag 1797, Penrith South DC, NSW 2751, Australia

<footnote n="1">¹ Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2, Ireland</footnote>

<footnote n="2">² Max-Planck-Institut für Kernphysik, PO Box 103980, 69029 Heidelberg, Germany</footnote>

<footnote n="3">³ Laboratoire Leprince-Ringuet, École Polytechnique, CNRS, Institut Polytechnique de Paris, 91128 Palaiseau, France</footnote>

<footnote n="4">⁴ University of Namibia, Department of Physics, Private Bag 13301, Windhoek 10005, Namibia</footnote>

<footnote n="5">⁵ Centre for Space Research, North-West University, Potchefstroom 2520, South Africa</footnote>

<footnote n="6">⁶ DESY, 15738 Zeuthen, Germany</footnote>

<footnote n="7">⁷ Université de Paris, CNRS, Astroparticule et Cosmologie, 75013 Paris, France</footnote>

<footnote n="8">⁸ Department of Physics and Electrical Engineering, Linnaeus University, 351 95 Växjö, Sweden</footnote>

<footnote n="9">⁹ Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany</footnote>

<footnote n="10">¹⁰ Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, 72076 Tübingen, Germany</footnote>

<footnote n="11">¹¹ Université Savoie Mont Blanc, CNRS, Laboratoire d'Annecy de Physique des Particules – IN2P3, 74000 Annecy, France</footnote>

<footnote n="12">¹² IRFU, CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="13">¹³ Astronomical Observatory, The University of Warsaw, Al. Ujazdowskie 4, 00-478 Warsaw, Poland</footnote>
## Appendix A: Contribution from Galactic diffuse γ-ray emission

We describe in this appendix our estimation for the contribution of Galactic diffuse emission to the γ-ray emission from HESS J1646−458, based on a prediction obtained with the PICARD CR propagation code (Kissmann 2014). The PICARD simulation is based on an analytical, continuous distribution of CR sources that follows the spiral arms of the Milky Way (for more details, see Kissmann et al. 2015). Kissmann et al. (2017) used these simulations to derive predictions for the flux of diffuse γ rays from three different processes: bremsstrahlung and IC emission from CR electrons, and the decay of neutral pions (and other short-lived particles) created in interactions of hadronic CRs with the interstellar medium. Template maps of these predictions were provided to us by R. Kissmann (priv. comm.). The predictions depend on various model assumptions (e.g. the CR source distribution, the distribution of gas and radiation fields, ...), meaning that in particular the predicted absolute flux levels are rather uncertain. Nevertheless, we used the templates as an estimate for the level of diffuse emission in the Westerlund 1 region, noting that contributions larger than predicted would quickly lead to a worse agreement between our background model (which predicts the level of hadronic background, but does not take into account diffuse γ-ray emission) and the observed data (cf. Figs. 1a and 2).

We show in Fig. A.1 the total predicted flux of diffuse γ rays (i.e. the sum of the flux for each of the aforementioned processes) in the energy range above 0.37 TeV. In this region and energy range, IC emission from CR electrons represents the dominant contribution to the Galactic diffuse emission. In Fig. A.2, we show the same flux maps as in Fig. 3, but with the diffuse γ-ray flux as predicted by the PICARD code in the respective energy range subtracted.

For a more quantitative comparison, we integrated the γ-ray flux measured with H.E.S.S. and that predicted as diffuse emission by PICARD in each of the 16 signal regions (a–p) defined in Fig. 1b, the resulting values are displayed in Table A.1. Above the lowest threshold energy of 0.37 TeV, the ratio between the predicted diffuse emission and the measured γ-ray flux varies between ~50% for region ‘a’ and ~12% for region ‘m’. For

higher energy thresholds we obtained smaller ratios in general, down to ~4% contamination above 4.9 TeV for region ‘j’. Summing up the fluxes over all signal regions, we found for energy thresholds of 0.37 TeV, 1 TeV, and 4.9 TeV a ratio of ~24%, ~17%, and ~8%, respectively. We stress again, however, that the absolute flux level of the diffuse emission is rather uncertain, and so are the estimates of its contribution to the total observed emission.

Finally, we note that, due to the omnipresence of the diffuse γ-ray emission along the Galactic plane, it is likely that at least part of it has been absorbed by the hadronic background model when we adjusted it to the observed data (cf. Sect. 2.1). This would reduce the measured γ-ray flux, and hence render the ratios $F_{\text{diff}}/F_{\text{meas}}$ listed in Table A.1 overestimates, even if $F_{\text{diff}}$ was precisely known.

<figure><img src="image_10.png" /><figcaption>Fig. A.1. Prediction of total diffuse γ-ray flux above 0.37 TeV from the PICARD code. The total emission includes contributions from bremsstrahlung, IC emission, and hadronic processes. The position of Westerlund 1 is marked by the black star symbol; the grey, dashed line shows the Galactic plane. Blue lines are contour lines of the flux map shown in Fig. 3a.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. A.2.** Flux maps after subtraction of Galactic diffuse γ-ray emission. The maps are the same as in Fig. 3, except that the Galactic diffuse γ-ray flux as predicted by the PICARD code has been subtracted. The colour axis scales and flux levels of contour lines are identical to those in Fig. 3.</figcaption></figure>

<table><caption>**Table A.1.** Comparison of measured γ-ray flux and predicted diffuse flux for signal regions, above different energy thresholds.</caption><thead><tr><th rowspan="2">Signal region</th><th colspan="3">E > 0.37 TeV</th><th colspan="3">E > 1 TeV</th><th colspan="3">E > 4.9 TeV</th></tr><tr><th>F<sub>meas</sub></th><th>F<sub>diff</sub></th><th>F<sub>diff</sub>/F<sub>meas</sub></th><th>F<sub>meas</sub></th><th>F<sub>diff</sub></th><th>F<sub>diff</sub>/F<sub>meas</sub></th><th>F<sub>meas</sub></th><th>F<sub>diff</sub></th><th>F<sub>diff</sub>/F<sub>meas</sub></th></tr></thead><tbody><tr><td>a</td><td>34.8</td><td>17.5</td><td>0.50</td><td>7.32</td><td>3.12</td><td>0.43</td><td>0.457</td><td>0.218</td><td>0.48</td></tr><tr><td>b</td><td>38.5</td><td>18.1</td><td>0.47</td><td>7.64</td><td>3.24</td><td>0.42</td><td>0.546</td><td>0.227</td><td>0.42</td></tr><tr><td>c</td><td>64.8</td><td>17.6</td><td>0.27</td><td>19.3</td><td>3.14</td><td>0.16</td><td>1.86</td><td>0.219</td><td>0.12</td></tr><tr><td>d</td><td>56.3</td><td>15.6</td><td>0.28</td><td>21.2</td><td>2.76</td><td>0.13</td><td>3.61</td><td>0.191</td><td>0.05</td></tr><tr><td>e</td><td>26.1</td><td>12.5</td><td>0.48</td><td>10.6</td><td>2.19</td><td>0.21</td><td>1.38</td><td>0.150</td><td>0.11</td></tr><tr><td>f</td><td>48.3</td><td>13.7</td><td>0.28</td><td>9.00</td><td>2.44</td><td>0.27</td><td>2.33</td><td>0.170</td><td>0.07</td></tr><tr><td>g</td><td>75.3</td><td>15.9</td><td>0.21</td><td>19.4</td><td>2.83</td><td>0.15</td><td>2.69</td><td>0.197</td><td>0.07</td></tr><tr><td>h</td><td>71.8</td><td>17.0</td><td>0.24</td><td>18.4</td><td>3.03</td><td>0.17</td><td>2.44</td><td>0.211</td><td>0.09</td></tr><tr><td>i</td><td>60.0</td><td>8.39</td><td>0.14</td><td>13.8</td><td>1.47</td><td>0.11</td><td>2.04</td><td>0.101</td><td>0.05</td></tr><tr><td>j</td><td>84.4</td><td>10.9</td><td>0.13</td><td>20.3</td><td>1.91</td><td>0.09</td><td>3.74</td><td>0.132</td><td>0.04</td></tr><tr><td>k</td><td>80.3</td><td>13.8</td><td>0.17</td><td>20.1</td><td>2.44</td><td>0.12</td><td>2.72</td><td>0.170</td><td>0.06</td></tr><tr><td>l</td><td>74.1</td><td>18.3</td><td>0.25</td><td>16.8</td><td>3.27</td><td>0.19</td><td>2.72</td><td>0.228</td><td>0.08</td></tr><tr><td>m</td><td>50.4</td><td>6.27</td><td>0.12</td><td>7.36</td><td>1.08</td><td>0.15</td><td>1.24</td><td>0.074</td><td>0.06</td></tr><tr><td>n</td><td>38.6</td><td>8.21</td><td>0.21</td><td>8.50</td><td>1.43</td><td>0.17</td><td>1.75</td><td>0.098</td><td>0.06</td></tr><tr><td>o</td><td>72.9</td><td>11.0</td><td>0.15</td><td>20.2</td><td>1.94</td><td>0.10</td><td>2.73</td><td>0.135</td><td>0.05</td></tr><tr><td>p</td><td>58.4</td><td>15.8</td><td>0.27</td><td>14.4</td><td>2.81</td><td>0.20</td><td>2.04</td><td>0.196</td><td>0.10</td></tr></tbody><tfoot><tr><td>All</td><td>935</td><td>221</td><td>0.24</td><td>234</td><td>39.1</td><td>0.17</td><td>34.3</td><td>2.72</td><td>0.08</td></tr></tfoot></table>

**Notes.** See Fig. 1b for the definition of the signal regions. $F_{\text{meas}}$ and $F_{\text{diff}}$ denote the measured and predicted diffuse γ-ray flux, respectively, in units of $10^{-11} \text{ cm}^{-2} \text{ s}^{-1}$.
**Appendix B: Radio maps for other velocity intervals**

As a comparison to Fig. 8, we show H I and CO maps for velocity intervals of $v = [-48.5, -38.5]$ km s⁻¹ and $v = [-37, -27]$ km s⁻¹ in Figs. B.1 and B.2, respectively.

<figure><img src="image_12.png" /><figcaption>**Fig. B.1.** Same as Fig. 8, but adopting an interval in velocity with respect to the local standard of rest of $v = [-48.5, -38.5]$ km s⁻¹.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. B.2.** Same as Fig. 8, but adopting an interval in velocity with respect to the local standard of rest of $v = [-37, -27]$ km s⁻¹.</figcaption></figure>