EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 666 (October 2022) A&A, 666 (2022) A124 Full HTML
Open Access
Issue
A&A
Volume 666, October 2022
Article Number A124
Number of page(s) 18
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202244323
Published online 17 October 2022

© F. Aharonian et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model.

Open access funding provided by Max Planck Society.

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 excavated 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 CRs1 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 γ rays, which provides strong motivation for the search for γ-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 γ-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, 2020b). 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 γ 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 γ-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 γ 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 positrons2 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) 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 105 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 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 Lw ∼ 1039 erg s−1 (Muno et al. 2006a) and for the wind velocity vw ∼ 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.

Table 1.

Parameter values for Westerlund 1 assumed in this work.

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. 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⊙) 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 × 10 35 $ \dot{E} > 2 \times 10^{35} $ erg s−1) 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 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.

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

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

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

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

thumbnail 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 ( TS > 3 $ \sqrt{\mathrm{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 16h49m4.8s, Dec −46° 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)×10−9 cm−2 s−1 sr−1 for panels a and b, at F = (3/5.5/8)×10−9 cm−2 s−1 sr−1 for panel c, and at F = (1/1.5)×10−9 cm−2 s−1 sr−1 for panel d.

thumbnail Fig. 4.

Radial excess profiles. Shown are exposure-corrected excess counts per unit sky area. Upper panel: profiles for different energy bands. Lower panel: profiles for different segments as defined in Fig. 1b. The black curve, showing the total excess above threshold in all segments, is the same in both panels. All profiles are normalised to equal area, to allow an easy comparison. The profiles have been computed with respect to a centre point at RA 16h46m36s, Dec −46° 01′12″, slightly shifted from the Westerlund 1 position. In the calculation of the profiles, we discarded pixels within 0.6° of the position of HESS J1640−465.

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. 2014c; 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. 2015, 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 predicted 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 energy – we derived radial profiles of the observed excess. Noting that Westerlund 1 is not located precisely at the centre of the shell-like structure, the profiles were computed not with respect to the stellar cluster position, but to a slightly shifted coordinate (RA 16h46m36s, 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. Figure 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 fall-off; (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.

Table 2.

Morphological fit results.

2.2.2. 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,

d N d E = ϕ 0 · ( E E 0 ) − Γ , $$ \begin{aligned} \frac{\mathrm{d} N}{\mathrm{d} E} = \phi _0\cdot \left(\frac{E}{E_0}\right)^{-\Gamma }, \end{aligned} $$(1)

with E0 = 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.

Table 3.

Spectral analysis results for signal regions.

thumbnail 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 Γ = 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σ) flux point are shown.

thumbnail Fig. 6.

Spectral index Γ 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, ⟨Γ⟩=2.41 ± 0.02, respectively.

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.

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

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

d N d E = ϕ 0 · ( E E 0 ) − Γ · exp ( − E E c ) , $$ \begin{aligned} \frac{\mathrm{d} N}{\mathrm{d} E} = \phi _0\cdot \left(\frac{E}{E_0}\right)^{-\Gamma }\cdot \exp \left(-\frac{E}{E_{\rm c}}\right), \end{aligned} $$(2)

for which we obtained (with E0 = 1 TeV kept fixed in the fit) ϕ0 = (1.00 ± 0.03)×1011 TeV−1 cm−2 s−1, Γ = 2.30 ± 0.04, and E c = ( 44 − 11 + 17 ) $ E_{\mathrm{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 × 1034 (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 E0 = 1 TeV) of ϕ 0 p = ( 1.28 ± 0.17 ) × 10 38 ( d / 3.9 k p c ) 2 ( n / 1 c m − 3 ) − 1 e V − 1 $ \phi_0^p = (1.28 \pm 0.17) \times 10^{38}\,(d/{3.9}\,{\rm kpc})^2\,(n/1\,{\rm cm}^{-3})^{-1}\,{\rm 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 c p = ( 400 − 130 + 250 ) $ E_{\mathrm{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 $ E_{\rm c}^p > {214} $ TeV. Extrapolating the proton spectrum down to an energy of 1 GeV, the required energy in primary protons is Wp ∼ 6 × 1051 (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 (TIR = 26 K, ρIR = 0.74 eV cm−3) and an optical field (Topt = 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 TSC = 40 000 K (Crowther et al. 2006) and derived3 an energy density of ρSC ∼ 30 eV cm3 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 ϕ 0 e = ( 4.7 ± 0.5 ) × 10 35 ( d / 3.9 k p c ) 2 e V − 1 $ \phi_0^e = (4.7 \pm 0.5)\times 10^{35}\,(d/{3.9}\,{\rm kpc})^2\,{\rm eV}^{-1} $, Γe = 2.97 ± 0.07, and E c e = ( 180 − 70 + 200 ) $ E_{\mathrm{c}}^e = (180^{+200}_{-70}) $ TeV, with a 95% confidence level lower limit on E c e $ E_{\rm c}^e $ of 87 TeV – the resulting γ-ray spectrum is shown 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 We ∼ 7.2 × 1048 (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 Le >  4.1 × 1035 (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, Le >  1.7 × 1036 (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 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 γ-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, H2 (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 velocity intervals, v = [ − 48.5, −38.5] km s−1 (d ≈ 3.3 kpc) and v = [ − 37, −27] km s−1 (d ≈ 2.7 kpc), are shown in Appendix B.

thumbnail Fig. 8.

Maps showing H I emission (McClure-Griffiths et al. 2005) (top panel) and 12CO 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−1, 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.

thumbnail 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] 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.

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 XH I = 1.823 × 1018 cm−2/(K km s−1; Rohlfs & 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 MH I, Wd1 = 1.3 × 105 M⊙. This translates into an average density of nH I, Wd1 = 3.2 cm−3 4. Similarly, from the CO data we get5MCO, Wd1 = 4.3 × 105 M⊙ and nCO, Wd1 = 10.5 cm−3, where nCO is the equivalent density for atomic hydrogen and can thus be directly compared to nH 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 gas – a density of nCO, 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 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 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 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).

3.1.2. PSR J1648–4611 and PSR J1650–4601

The two energetic pulsars PSR J1648−4611 (with spin-down power E ˙ = 2.1 × 10 35 $ \dot{E} = 2.1 \times 10^{35} $ erg s−1) and PSR J1650−4601 ( E ˙ = 2.9 × 10 35 $ \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, Ṗ = 9.2 × 10−13 s s−1; Israel et al. 2007) imply a rotational spin-down power of only E ˙ = 3 × 10 31 $ \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 LX ≈ 3 × 1033 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 1051 erg – possibly more in the case of CXOU J164710.2−455216, if its progenitor was indeed as heavy as 40 M⊙ (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 mH cm−3 gives a required energy in protons of 1.2 × 1051 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 ≥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/RWd1)2/3 times larger (Longair 1992), where RWd1 ∼ 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).

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 ∼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 RSB = 0.76 (Lw/ρ0)1/5τ3/5 (Weaver et al. 1977; Koo & McKee 1992), or RSB ∼ 256 (ρ0/1 mH cm−3)−1/5 pc for our assumptions (cf. Table 1). Adopting, for example, ρ0 = 5 mH cm−3, we obtain RSB ∼ 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).

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 / ρ 0 ) 3 / 10 v w − 1 / 2 τ 2 / 5 $ R_{\mathrm{TS}} = 0.92\,(L_{\mathrm{w}}/\rho_0)^{3/10} v_{\mathrm{w}}^{-1/2} \tau^{2/5} $ (Koo & McKee 1992), or RTS ∼ 51 (ρ0/1 mHcm−3)−3/10 pc with our adopted parameter values. Inserting ρ0 = 5 mH cm−3 yields a radius of RTS ∼ 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 Lw ∼ 1039 erg s−1 and τ ∼ 4 Myr, the total available energy is ∼1.3 × 1053 erg, which in principle suffices to explain the required energy in CR protons, Wp ∼ 6 × 1051 (n/1 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 ∼ 6 D t ∼ 81 ( E / 100 TeV ) 1 / 2 ( B / 10 μ G ) − 1 / 2 ( t / 1 Myr ) 1 / 2 $ L \sim \sqrt{6Dt} \sim 81\,(E/{100}\,\mathrm{TeV})^{1/2}\,(B/10\,{\upmu}\mathrm{G})^{-1/2}\,(t/{1}\,\mathrm{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 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 Le ∼ 1.7 × 1036 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 GHz6, 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 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 Γe ∼ 3.07 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.

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

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.

2

Hereafter, we use the term ‘electrons’ to refer to both electrons and positrons.

3

To derive the energy density, we used ρSC = L/(4πr2c), where L is the total cluster luminosity and r the distance from the cluster. For a wind efficiency η = Ṁvw/(L/c) ≃ 1 (Vink & Gräfener 2012), L = 2 (vw/c)−1Lw ∼ 2 × 1041 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 × 1040 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.

4

Abramowski et al. (2012) derived, for a similar region, a much smaller value of nH 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-H2 conversion factor, XCO, is less well constrained than XH I. Here we used XCO = 1.5 × 1020 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)×1020 cm−2/(K km s−1) recommended by Bolatto et al. (2013). We have neglected the possible contribution of 4He, which could increase the mass estimate by ∼25%.

6

We have used the full-sky frequency map at 30 GHz available through the Planck Legacy Archive, http://pla.esac.esa.int/pla

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.

Acknowledgments

The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of H.E.S.S. is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the German Research Foundation (DFG), the Helmholtz Association, the Alexander von Humboldt Foundation, the French Ministry of Higher Education, Research and Innovation, the Centre National de la Recherche Scientifique (CNRS/IN2P3 and CNRS/INSU), the Commissariat à l’Énergie atomique et aux Énergies alternatives (CEA), the 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

  1. Abdalla, H., Abramowski, A., Aharonian, F., et al. 2018a, A&A, 612, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Abdalla, H., Aharonian, F., Ait Benkhali, F., et al. 2018b, A&A, 620, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Abdalla, H., Abramowski, A., Aharonian, F., et al. 2018c, A&A, 612, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Abdalla, H., Abramowski, A., Aharonian, F., et al. 2018d, A&A, 612, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Abdalla, H., Aharonian, F., Ait Benkhali, F., et al. 2019, A&A, 621, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Abdalla, H., Adam, R., Aharonian, F., et al. 2020, A&A, 635, A167 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Abdo, A. A., Allen, B., Berley, D., et al. 2007, ApJ, 658, L33 [NASA ADS] [CrossRef] [Google Scholar]
  8. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 723, 649 [NASA ADS] [CrossRef] [Google Scholar]
  9. Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33 [Google Scholar]
  10. Abeysekara, A. U., Archer, A., Aune, T., et al. 2018, ApJ, 861, 134 [NASA ADS] [CrossRef] [Google Scholar]
  11. Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2021, Nat. Astron., 5, 465 [NASA ADS] [CrossRef] [Google Scholar]
  12. Abramowski, A., Acero, F., Aharonian, F., et al. 2011, A&A, 525, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Abramowski, A., Acero, F., Aharonian, F., et al. 2012, A&A, 537, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Abramowski, A., Aharonian, F., Ait Benkhali, F., et al. 2014a, MNRAS, 439, 2828 [NASA ADS] [CrossRef] [Google Scholar]
  15. Abramowski, A., Aharonian, F., Ait Benkhali, F., et al. 2014b, ApJ, 794, L1 [NASA ADS] [CrossRef] [Google Scholar]
  16. Abramowski, A., Aharonian, F., Ait Benkhali, F., et al. 2014c, Phys. Rev. D, 90, 122007 [NASA ADS] [CrossRef] [Google Scholar]
  17. Abramowski, A., Aharonian, F., Ait Benkhali, F., et al. 2015, Science, 347, 406 [NASA ADS] [CrossRef] [Google Scholar]
  18. Acharya, B. S., Agudo, I., Al Samarai, I., et al. 2018, Science with the Cherenkov Telescope Array (Singapore: World Scientific Publishing) [Google Scholar]
  19. Ackermann, M., Ajello, M., Allafort, A., et al. 2011, Science, 334, 1103 [NASA ADS] [CrossRef] [Google Scholar]
  20. Ackermann, M., Ajello, M., Atwood, W. B., et al. 2012, ApJ, 750, 3 [NASA ADS] [CrossRef] [Google Scholar]
  21. Aghakhanloo, M., Murphy, J. W., Smith, N., et al. 2020, MNRAS, 492, 2497 [NASA ADS] [CrossRef] [Google Scholar]
  22. Aghakhanloo, M., Murphy, J. W., Smith, N., et al. 2021, Res. Notes Am. Astron. Soc., 5, 14 [NASA ADS] [Google Scholar]
  23. Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, A&A, 457, 899 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007, A&A, 467, 1075 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Aharonian, F., Yang, R., & de Oña Wilhelmi, E. 2019, Nat. Astron., 3, 561 [CrossRef] [Google Scholar]
  26. Akrami, Y., Argüeso, F., Ashdown, M., et al. 2020, A&A, 641, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Amenomori, M., Bao, Y. W., Bi, X. J., et al. 2021, Phys. Rev. Lett., 126, 141101 [NASA ADS] [CrossRef] [Google Scholar]
  28. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  30. Ballet, J., Burnett, T. H., Digel, S. W., & Lott, B. 2020, ArXiv e-prints [arXiv:2005.11208] [Google Scholar]
  31. Beasor, E. R., Davies, B., Smith, N., Gehrz, R. D., & Figer, D. F. 2021, ApJ, 912, 16 [NASA ADS] [CrossRef] [Google Scholar]
  32. Belczynski, K., & Taam, R. E. 2008, ApJ, 685, 400 [NASA ADS] [CrossRef] [Google Scholar]
  33. Bell, A. R. 2013, Astropart. Phys., 43, 56 [NASA ADS] [CrossRef] [Google Scholar]
  34. Bhadra, S., Gupta, S., Nath, B. B., & Sharma, P. 2022, MNRAS, 510, 5579 [NASA ADS] [CrossRef] [Google Scholar]
  35. Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207 [CrossRef] [Google Scholar]
  36. Braiding, C., Wong, G. F., Maxted, N. I., et al. 2018, PASA, 35, e029 [NASA ADS] [CrossRef] [Google Scholar]
  37. Brandner, W., Clark, J. S., Stolte, A., et al. 2008, A&A, 478, 137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Breuhaus, M., Hinton, J. A., Joshi, V., Reville, B., & Schoorlemmer, H. 2022, A&A, 661, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Bykov, A. M. 2014, A&ARv, 22, 77 [NASA ADS] [CrossRef] [Google Scholar]
  40. Bykov, A. M., Ellison, D. C., Gladilin, P. E., & Osipov, S. M. 2015, MNRAS, 453, 113 [NASA ADS] [CrossRef] [Google Scholar]
  41. Bykov, A. M., Marcowith, A., Amato, E., et al. 2020, Space Sci. Rev., 216, 42 [CrossRef] [Google Scholar]
  42. Cao, Z., Aharonian, F. A., An, Q., et al. 2021, Nature, 594, 33 [NASA ADS] [CrossRef] [Google Scholar]
  43. Cesarsky, C. J., & Montmerle, T. 1983, Space Sci. Rev., 36, 173 [NASA ADS] [CrossRef] [Google Scholar]
  44. Chandrasekhar, S. 1943, Rev. Mod. Phys., 15, 1 [CrossRef] [Google Scholar]
  45. Chu, Y.-H. 2008, in Proc. Int. Astron. Union, eds. F. Bresolin, P. A. Crowther, & J. Puls, 250, 341 [NASA ADS] [Google Scholar]
  46. Clark, J. S., Negueruela, I., Crowther, P., & Goodwin, S. P. 2005, A&A, 434, 949 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Clark, J. S., Muno, M. P., Negueruela, I., et al. 2008, A&A, 477, 147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Clark, J. S., Najarro, F., Negueruela, I., et al. 2019, A&A, 623, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Crowther, P. A., Hadfield, L. J., Clark, J. S., Negueruela, I., & Vacca, W. D. 2006, MNRAS, 372, 1407 [NASA ADS] [CrossRef] [Google Scholar]
  50. Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792 [NASA ADS] [CrossRef] [Google Scholar]
  51. Davies, B., & Beasor, E. R. 2019, MNRAS, 486, L10 [NASA ADS] [CrossRef] [Google Scholar]
  52. de Naurois, M., & Rolland, L. 2009, Astropart. Phys., 32, 231 [NASA ADS] [CrossRef] [Google Scholar]
  53. Deil, C., Zanin, R., Lefaucheur, J., et al. 2017, Proc. 35th Int. Cosmic Ray Conf. (ICRC2017), 766 [Google Scholar]
  54. Deil, C., Wood, M., Hassan, T., et al. 2018, https://doi.org/10.5281/zenodo.1409831 [Google Scholar]
  55. Deil, C., Donath, A., Terrier, R., et al. 2020, https://doi.org/10.5281/zenodo.4701492 [Google Scholar]
  56. Eichler, D., & Usov, V. 1993, ApJ, 402, 271 [NASA ADS] [CrossRef] [Google Scholar]
  57. Forman, W., Jones, C., Cominsky, L., et al. 1978, ApJS, 38, 357 [Google Scholar]
  58. Fujii, M. S., & Portegies Zwart, S. 2016, ApJ, 817, 4 [NASA ADS] [CrossRef] [Google Scholar]
  59. Gupta, S., Nath, B. B., Sharma, P., & Eichler, D. 2018, MNRAS, 473, 1537 [NASA ADS] [Google Scholar]
  60. Heyer, M., & Dame, T. M. 2015, ARA&A, 53, 583 [NASA ADS] [CrossRef] [Google Scholar]
  61. Hillas, A. M. 1984, ARA&A, 22, 425 [Google Scholar]
  62. Hinton, J. A., & Hofmann, W. 2009, ARA&A, 47, 523 [NASA ADS] [CrossRef] [Google Scholar]
  63. Holler, M., Berge, D., van Eldik, C., et al. 2015, Proc. 34th Int. Cosmic Ray Conf. (ICRC2015), 847 [Google Scholar]
  64. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  65. Israel, G. L., Campana, S., Dall’Osso, S., et al. 2007, ApJ, 664, 448 [NASA ADS] [CrossRef] [Google Scholar]
  66. Kantzas, D., Markoff, S., Lucchini, M., et al. 2022, MNRAS, 510, 5187 [NASA ADS] [CrossRef] [Google Scholar]
  67. Kavanagh, P. J., Norci, L., & Meurs, E. J. A. 2011, New Astron., 16, 461 [NASA ADS] [CrossRef] [Google Scholar]
  68. Kavanagh, P. J., Vink, J., Sasaki, M., et al. 2019, A&A, 621, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Kissmann, R. 2014, Astropart. Phys., 55, 37 [NASA ADS] [CrossRef] [Google Scholar]
  70. Kissmann, R., Werner, M., Reimer, O., & Strong, A. W. 2015, Astropart. Phys., 70, 39 [NASA ADS] [CrossRef] [Google Scholar]
  71. Kissmann, R., Niederwanger, F., Reimer, O., & Strong, A. W. 2017, AIP Conf. Proc., 1792, 070011 [NASA ADS] [CrossRef] [Google Scholar]
  72. Knödlseder, J., Mayer, M., Deil, C., et al. 2016, A&A, 593, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Koo, B.-C., & McKee, C. F. 1992, ApJ, 388, 93 [NASA ADS] [CrossRef] [Google Scholar]
  74. Kothes, R., & Dougherty, S. M. 2007, A&A, 468, 993 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Kramer, M., Bell, J. F., Manchester, R. N., et al. 2003, MNRAS, 342, 1299 [NASA ADS] [CrossRef] [Google Scholar]
  76. Li, T., & Ma, Y. 1983, ApJ, 272, 317 [CrossRef] [Google Scholar]
  77. Longair, M. S. 1992, High Energy Astrophysics. Vol. 1: Particles, Photons and their Detection (Cambridge: Cambridge University Press) [Google Scholar]
  78. Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993 [NASA ADS] [CrossRef] [Google Scholar]
  79. Martí-Devesa, G., & Reimer, O. 2021, A&A, 654, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. McClure-Griffiths, N. M., Dickey, J. M., Gaensler, B. M., et al. 2005, ApJS, 158, 178 [NASA ADS] [CrossRef] [Google Scholar]
  81. Mohrmann, L., Specovius, A., Tiziani, D., et al. 2019, A&A, 632, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Morlino, G., Blasi, P., Peretti, E., & Cristofari, P. 2021, MNRAS, 504, 6096 [NASA ADS] [CrossRef] [Google Scholar]
  83. Muno, M. P., Law, C., Clark, J. S., et al. 2006a, ApJ, 650, 203 [NASA ADS] [CrossRef] [Google Scholar]
  84. Muno, M. P., Clark, J. S., Crowther, P. A., et al. 2006b, ApJ, 636, L41 [NASA ADS] [CrossRef] [Google Scholar]
  85. Negueruela, I., Alfaro, E. J., Dorda, R., et al. 2022, A&A, 664, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Ohm, S., van Eldik, C., & Egberts, K. 2009, Astropart. Phys., 31, 383 [NASA ADS] [CrossRef] [Google Scholar]
  87. Ohm, S., Hinton, J. A., & White, R. 2013, MNRAS, 434, 2289 [NASA ADS] [CrossRef] [Google Scholar]
  88. Parizot, E., Marcowith, A., van der Swaluw, E., Bykov, A. M., & Tatischeff, V. 2004, A&A, 424, 747 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Parsons, R. D., & Hinton, J. A. 2014, Astropart. Phys., 56, 26 [NASA ADS] [CrossRef] [Google Scholar]
  90. Piron, F., Djannati-Ataï, A., Punch, M., et al. 2001, A&A, 374, 895 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Popescu, C. C., Yang, R., Tuffs, R. J., et al. 2017, MNRAS, 470, 2539 [NASA ADS] [CrossRef] [Google Scholar]
  92. Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431 [NASA ADS] [CrossRef] [Google Scholar]
  93. Pshirkov, M. S. 2016, MNRAS, 457, L99 [NASA ADS] [CrossRef] [Google Scholar]
  94. Rate, G., Crowther, P. A., & Parker, R. J. 2020, MNRAS, 495, 1209 [NASA ADS] [CrossRef] [Google Scholar]
  95. Reimer, A., Pohl, M., & Reimer, O. 2006, ApJ, 644, 1118 [NASA ADS] [CrossRef] [Google Scholar]
  96. Rohlfs, K., & Wilson, T. L. 2004, Tools of Radio Astronomy, 4th edn. (Berlin, Heidelberg: Springer-Verlag) [CrossRef] [Google Scholar]
  97. Sakai, M., Matsumoto, H., Haba, Y., Kanou, Y., & Miyamoto, Y. 2013, PASJ, 65, 64 [NASA ADS] [CrossRef] [Google Scholar]
  98. Specovius, A. 2021, PhD Thesis, Friedrich-Alexander-Universität Erlangen-Nürnberg [Google Scholar]
  99. Sun, X.-N., Yang, R.-Z., Liang, Y.-F., et al. 2020a, A&A, 639, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  100. Sun, X.-N., Yang, R.-Z., & Wang, X.-Y. 2020b, MNRAS, 494, 3405 [NASA ADS] [CrossRef] [Google Scholar]
  101. Tavani, M., Sabatini, S., Pian, E., et al. 2009, ApJ, 698, L142 [NASA ADS] [CrossRef] [Google Scholar]
  102. Vieu, T., Gabici, S., Tatischeff, V., & Ravikularaman, S. 2022, MNRAS, 512, 1275 [NASA ADS] [CrossRef] [Google Scholar]
  103. Vink, J. S., & Gräfener, G. 2012, ApJ, 751, L34 [NASA ADS] [CrossRef] [Google Scholar]
  104. Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377 [NASA ADS] [CrossRef] [Google Scholar]
  105. Werner, M., Reimer, O., Reimer, A., & Egberts, K. 2013, A&A, 555, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Westerlund, B. 1961, PASP, 73, 51 [NASA ADS] [CrossRef] [Google Scholar]
  107. Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010, ApJ, 716, 1191 [NASA ADS] [CrossRef] [Google Scholar]
  108. Yadav, N., Mukherjee, D., Sharma, P., & Nath, B. B. 2017, MNRAS, 465, 1720 [NASA ADS] [CrossRef] [Google Scholar]
  109. Yang, R.-Z., & Aharonian, F. 2017, A&A, 600, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  110. Yang, R.-Z., & Wang, Y. 2020, A&A, 640, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  111. Yang, R., de Oña Wilhelmi, E., & Aharonian, F. 2018, A&A, 611, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  112. Zabalza, V. 2015, Proc. 34th Int. Cosmic Ray Conf. (ICRC2015), 922 [Google Scholar]
  113. Zorn, J. 2019, PhD Thesis, Ruprecht-Karls-Universität Heidelberg [Google Scholar]

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 Fdiff/Fmeas listed in Table A.1 overestimates, even if Fdiff was precisely known.

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

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

Table A.1.

Comparison of measured γ-ray flux and predicted diffuse flux for signal regions, above different energy thresholds.

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−1 and v = [ − 37, −27] km s−1 in Figs. B.1 and B.2, respectively.

thumbnail 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−1.

thumbnail 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−1.

All Tables

Table 1.

Parameter values for Westerlund 1 assumed in this work.

In the text
Table 2.

Morphological fit results.

In the text
Table 3.

Spectral analysis results for signal regions.

In the text
Table A.1.

Comparison of measured γ-ray flux and predicted diffuse flux for signal regions, above different energy thresholds.

In the text

All Figures

thumbnail 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).
In the text
thumbnail 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.
In the text
thumbnail 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 ( TS > 3 $ \sqrt{\mathrm{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 16h49m4.8s, Dec −46° 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)×10−9 cm−2 s−1 sr−1 for panels a and b, at F = (3/5.5/8)×10−9 cm−2 s−1 sr−1 for panel c, and at F = (1/1.5)×10−9 cm−2 s−1 sr−1 for panel d.
In the text
thumbnail Fig. 4.

Radial excess profiles. Shown are exposure-corrected excess counts per unit sky area. Upper panel: profiles for different energy bands. Lower panel: profiles for different segments as defined in Fig. 1b. The black curve, showing the total excess above threshold in all segments, is the same in both panels. All profiles are normalised to equal area, to allow an easy comparison. The profiles have been computed with respect to a centre point at RA 16h46m36s, Dec −46° 01′12″, slightly shifted from the Westerlund 1 position. In the calculation of the profiles, we discarded pixels within 0.6° of the position of HESS J1640−465.
In the text
thumbnail 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 Γ = 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σ) flux point are shown.
In the text
thumbnail Fig. 6.

Spectral index Γ 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, ⟨Γ⟩=2.41 ± 0.02, respectively.
In the text
thumbnail 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.
In the text
thumbnail Fig. 8.

Maps showing H I emission (McClure-Griffiths et al. 2005) (top panel) and 12CO 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−1, 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.
In the text
thumbnail 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] 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.
In the text
thumbnail 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.
In the text
thumbnail 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.
In the text
thumbnail 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−1.
In the text
thumbnail 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−1.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.