A deep spectromorphological study of the $\gamma$-ray emission surrounding the young massive stellar cluster Westerlund 1

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) 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 $\gamma$-ray source HESS J1646-458 has been detected in the vicinity of Westerlund 1 in the past, its association could not be firmly identified. We aim to identify the physical processes responsible for the $\gamma$-ray emission around Westerlund 1 and thus to better understand the role of massive stellar clusters in the acceleration of Galactic CRs. Using 164 hours of data recorded with the High Energy Stereoscopic System (H.E.S.S.), we carried out a deep spectromorphological study of the $\gamma$-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. We detected large-scale ($\sim 2^\circ$ diameter) $\gamma$-ray emission with a complex morphology, exhibiting a shell-like structure and showing no significant variation with $\gamma$-ray energy. The combined energy spectrum of the emission extends to several tens of TeV, and is uniform across the entire source region. We did not find a clear correlation of the $\gamma$-ray emission with gas clouds as identified through H I and CO observations. We conclude that, of the known objects within the region, only Westerlund 1 can explain the bulk of the $\gamma$-ray emission. Several CR acceleration sites and mechanisms are conceivable, and discussed in detail. (abridged)


Introduction
Young massive stellar clusters are environments of copious star formation, and 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 excavated by the com-Through interaction with ambient gas and radiation fields, high-energy hadronic CRs produce high-energy γ rays, which provides strong motivation for the search for γ-ray emission from massive stellar clusters (this has been realised already 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 γ-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 lead to several detections: the Cygnus region harbours several sources of TeV-energy γ rays (Abdo et al. 2007;Abeysekara et al. 2018), and has recently been detected up to energies of hundreds of TeV (Abeysekara et al. 2021;Cao et al. 2021); -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); -the young stellar cluster Westerlund 1 (Abramowski et al. 2012), which will be introduced in more detail below; -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); -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 γ-ray emission can be linked directly to the stellar clusters in only some of the above cases, the precise CR acceleration sites are yet unidentified, 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 γ 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 and positrons 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 γ-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) (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 (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  (3,4) L w kinetic luminosity of cluster wind 10 39 erg s −1 (5) v w velocity of cluster wind 3000 km s −1 (6) References.
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(Aghakhanloo et al. , 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 ∼10 39 erg s −1 (Muno et al. 2006b) and for the wind velocity v w ∼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 (Muno et al. 2006b), which was later identified as likely being of thermal origin with XMM-Newton observations (Kavanagh et al. 2011). Additionally, Muno et al. (2006a) have identified an X-ray magnetar, CXOU J164710.2−455216, which supposedly was created in the explosion of a very massive (> 40 M ) progenitor star 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 (Ė > 2 × 10 35 erg s −1 ) pulsars, PSR J1648−4611 and PSR J1650−4601 (Manchester et al. 2005), as well as the lowmass 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 Article number, page 2 of 19 F. Aharonian et al.: A deep spectromorphological study of the γ-ray emission surrounding Westerlund 1 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 H I 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.

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. 2018d).
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 andOctober 11, 2017. The runs amount to a total observation time of 164.2 hours, which represents a significant increase with respect to the previous publication (Abramowski et al. 2012) (33.8 hours). 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. 2014c) 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(Deil et al. , 2020. 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 CRinduced 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 3dimensional 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 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. (2018b) was employed to generate an exclusion map.

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 Article number, page 3 of 19  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.
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 smooth-ing 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 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 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. 2014a;Amenomori et al. 2021), these measurements do not provide a good constraint for the level of diffuse emission in the region of HESS J1646−458. Therefore, in order to assess the possible contamination with diffuse emission of the γ-ray signal of HESS J1646−458, we have used a prediction of the diffuse γ-ray flux based on the Picard CR propagation code (Kissmann 2014;Kissmann et al. 2015Kissmann et al. , 2017. This analysis is described in more detail in Appendix A, where we show in Fig. A.2 the same flux maps as in Fig. 3, but with the pre-dicted flux due to diffuse emission subtracted. We conclude that, while the Galactic diffuse emission likely contributes at a considerable level -∼24% (∼17%/∼8%) above a threshold energy of 0.37 TeV (1 TeV/4.9 TeV), according to the Picard template -, it cannot explain the bulk of the γ-ray emission, and does not alter the source morphology in a significant way. For these reasons, and because of the rather large uncertainties associated with any estimate of the Galactic diffuse emission in a particular region of the sky, we have performed the subsequent analysis without explicitly taking it into account, noting that none of the conclusions drawn in this paper are affected by this.
In order to further characterise the morphology of the emission -and its apparent invariance with respect to energywe derived radial profiles of the observed excess. Noting that Article number, page 5 of 19 A&A proofs: manuscript no. westerlund1_paper Notes. 'Excess' denotes the number of observed excess events within the white, dashed circle in Fig. 1b, excluding a circular region with 0.6 • radius around HESS J1640−465. The χ 2 values and corresponding p-values are a measure of the compatibility of the shape of the radial excess profile of this energy band with the total radial excess profile (see text for details).
Westerlund 1 is not located precisely at the centre of the shelllike structure, the profiles were computed not with respect to the stellar cluster position, but to a slightly shifted coordinate (R.A. 16 h 46 m 36 s , Dec. −46 • 01 12 ), which corresponds to the barycentre of the γ-ray excess. The asymmetry of the observed emission with respect to the cluster position could for example be caused by inhomogeneities in the ISM surrounding the stellar cluster. Fig. 4 shows the profiles for different energy bands (upper panel), and for the five segments defined in Fig. 1b (lower panel). The excess profiles: (i) confirm the shell-like structure, with a peak at around 0.5 • (corresponding to a projected distance of ∼34 pc for a cluster distance of 3.9 kpc), followed by a relatively slow falloff; (ii) exhibit the same shape in all energy bands, that is, they show no indications for an energy-dependent morphology of the excess; (iii) are also largely compatible between the five segments, with only minor small-scale deviations discernible.
In order to reinforce the second point above, we carried out χ 2 tests in which we compared the profile for each energy band with one computed using all events outside this band (thus ensuring statistically independent sets). The results, listed in Table 2, show that each of the profiles is compatible with the total profile in terms of its shape within the statistical uncertainties.

Energy spectra
The complex morphology of HESS J1646−458 prohibits a simple extraction of an energy spectrum by means of modelling the emission with a single spatial model. We therefore considered the 16 signal regions indicated in Fig. 1b and extracted energy spectra for each of these regions. To obtain the spectra, we fitted (using a forward-folding likelihood fit; Piron et al. 2001) a power-law (PL) model, with E 0 = 1 TeV kept fixed in the fit, to the observed γ-ray excess in each signal region and subsequently derived flux points based on this fitted model. Table 3 lists the observed γ-ray excess as well as the fitted power-law model parameters. A comparison of the shapes of the energy spectra is provided in Fig. 5. Finally, Fig. 6 shows the fitted spectral index for each signal region as a function of its angular separation from Westerlund 1.  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 ap-Article number, page 6 of 19 F. Aharonian et al.: A deep spectromorphological study of the γ-ray emission surrounding Westerlund 1 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. φ 0 and Γ are the best-fit parameter values of the power-law fit for each region (cf. Eq. (1)). √ ∆TS denotes the square root of the difference in test statistic (TS = −2 ln(L)) between the best-fit power-law model and the null hypothesis (corresponding to φ 0 = 0). plied 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), for which we obtained (with E 0 = 1 TeV kept fixed in the fit) φ 0 = (1.00 ± 0.03) × 10 −11 TeV −1 cm −2 s −1 , Γ = 2.30 ± 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 γ ∼ 9 × 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 highenergy 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 pro-ton spectrum we obtained a normalisation (at E 0 = 1 TeV) of φ p 0 = (1.28 ± 0.17) × 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 Γ p = 2.33 ± 0.06, and a cut-off energy of E p c = (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 p c > 214 TeV. Extrapolating the proton spectrum down to an energy of 1 GeV, the required energy in primary protons is W p ∼ 6 × 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, ρ IR = 0.74 eV cm −3 ) and an optical field (T opt = 2400 K, ρ opt = 1.4 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 SC = 40 000 K (Crowther et al. 2006) and derived 3 an energy density of ρ SC ∼ 30 eV cm −3 at a distance of 34 pc from the cluster, which approximately corresponds to the distance at which the radial γ-ray excess profile peaks (cf. Fig. 4). In this case, the best-fit parameters are φ e 0 = (4.7 ± 0.5) × 10 35 (d/3.9 kpc) 2 eV −1 , Γ e = 2.97 ± 0.07, and E e c = (180 +200 −70 ) TeV, with a 95% confidence level lower limit on E e c of 87 TeV -the resulting γ-ray spectrum is shown 3 To derive the energy density, we used ρ SC = L / (4πr 2 c), where L is the total cluster luminosity and r the distance from the cluster. For a wind efficiency η =Ṁv w / (L/c) 1 (Vink & Gräfener 2012), L = 2 (v w /c) −1 L w ∼ 2 × 10 41 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 ∼ 6 × 10 40 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. by the red, dashed-dotted line in Fig. 7. Electrons down to energies of ∼0.4 TeV contribute to the γ-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 ∼ 7.2 × 10 48 (d/3.9 kpc) 2 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 × 10 35 (d/3.9 kpc) 2 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 µG, L e > 1.7 × 10 36 (d/3.9 kpc) 2 erg s −1 .
The chosen analysis method and tools in principle also enable a three-dimensional modelling of the γ-ray emission detected with H.E.S.S., that is, to decompose the emission into several components with separate energy spectra and morpho- logical 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.

Analysis of radio line data
Attributing the γ-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 (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 −1 , 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 ve- 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 HI = 1.823 × 10 18 cm −2 / (K km s −1 ) (  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. The bottom panel shows the 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. sity of n H I,Wd1 = 3.2 cm −3 . 4 Similarly, from the CO data we get 5 M CO,Wd1 = 4.3 × 10 5 M and n CO,Wd1 = 10.5 cm −3 , where n CO is the equivalent density for atomic hydrogen and can thus be directly compared to n H I . We stress, however, that in particular 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 gasa density of n CO,cloud ∼ 190 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 4 Abramowski et al. (2012) derived, for a similar region, a much smaller value of n H I = 0.24 cm −3 . We attribute this to the usage of an erroneous formula in that paper. 5 Due to the more indirect nature of the estimate, the CO-to-H 2 conversion factor, X CO , is less well constrained than X H I . Here we used X CO = 1.5 × 10 20 cm −2 / (K km s −1 ), which Ackermann et al. (2012) indicate as an appropriate value for the galactocentric radius of Westerlund 1, R ≈ 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 −2 / (K km s −1 ) recommended by Bolatto et al. (2013). We have neglected the possible contribution of 4 He, which could increase the mass estimate by ∼25%.  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 ≈ 3.9 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 de-Article number, page 9 of 19 A&A proofs: manuscript no. westerlund1_paper rived 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, 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).

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.

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

PSR J1648-4611 & PSR J1650-4601
The two energetic pulsars PSR J1648−4611 (with spin-down powerĖ = 2.1 × 10 35 erg s −1 ) and PSR J1650−4601 (Ė = 2.9 × 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.  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] 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. 2018b), 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.

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.

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,Ṗ = 9.2 × 10 −13 s s −1 ; Israel et al. 2007) imply a rotational spin-down power of onlyĖ = 3 × 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. 2018c). Even though the measured X-ray luminosity of L X ≈ 3 × 10 33 erg s −1 (Muno et al. 2006a) 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.

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. (2006b) 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 ergpossibly more in the case of CXOU J164710.2−455216, if its progenitor was indeed as heavy as 40 M (Muno et al. 2006a) -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 H cm −3 gives a required energy in protons of 1.2 × 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%.

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 ≥ 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 H I 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 ∼ 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 ∼ 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).

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.

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 ∼100 pc, the Hillas criterion (Hillas 1984) implies for an average turbulent fluid velocity u = 100 km s −1 a minimum magnetic field strength of ∼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 ρ 0 , the radius of the superbubble is given by R SB = 0.76 (L w / ρ 0 ) 1/5 τ 3/5 (Weaver et al. 1977;Koo & McKee 1992), or R SB ∼ 256 (ρ 0 / 1 m H cm −3 ) −1/5 pc for our assumptions (cf. Table 1). Adopting, for example, ρ 0 = 5 m H cm −3 , we obtain R SB ∼ 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 association 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 H I 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).

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 / ρ 0 ) 3/10 v −1/2 w τ 2/5 (Koo & McKee 1992), or R TS ∼ 51 (ρ 0 / 1 m H cm −3 ) −3/10 pc with our adopted parameter values. Inserting ρ 0 = 5 m H cm −3 yields a radius of R TS ∼ 32 pc. Notably, this value coincides well with the radial distance of ∼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 ∼ 10 39 erg s −1 and τ ∼ 4 Myr, the total available energy is ∼1.3 × 10 53 erg, which in principle suffices to explain the required energy in CR protons, W p ∼ 6 × 10 51 (n/1 cm −3 ) −1 ergalthough, 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 ∼ √ 6Dt ∼ 81 (E/100 TeV) 1/2 (B/10 µG) −1/2 (t/1 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 ∼50 pc for only 1 Myr already requires a rather large magnetic field strength of B ∼ 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 possible an interpretation within the leptonic scenario. 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 ∼ 1.7 × 10 36 erg s −1 (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 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 −1 . 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 −1 . 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 lowenergy 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).

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 Γ e ∼3.0 7 and Γ p ∼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.

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 positionin 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 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.
card 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 Kissmann (2022). 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 diff /F meas listed in Table A 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.
Article number, page 17 of 19 A&A proofs: manuscript no. westerlund1_paper    Article number, page 19 of 19