Understanding the TeV $\gamma$-ray emission surrounding the young massive star cluster Westerlund 1

Context: Young massive star clusters (YMCs) have come increasingly into the focus of discussions on the origin of galactic cosmic rays (CRs). The proposition of CR acceleration inside superbubbles (SBs) blown by the strong winds of these clusters avoids issues faced by the standard paradigm of acceleration at supernova remnant shocks. Aims: We provide an interpretation of the latest TeV $\gamma$-ray observations of the region around the YMC Westerlund 1 taken with the High Energy Stereoscopic System (H.E.S.S.) in terms of diffusive shock acceleration at the cluster wind termination shock, taking into account the spectrum and morphology of the emission. As Westerlund 1 is a prototypical example of a YMC, such a study is relevant to the general question about the role of YMCs for the Galactic CR population. Methods: We generate model $\gamma$-ray spectra, characterise particle propagation inside the SB based on the advection, diffusion, and cooling timescales, and constrain key parameters of the system. We consider hadronic emission from proton-proton interaction and subsequent pion decay and leptonic emission from inverse Compton scattering on all relevant photon fields, including the CMB, diffuse and dust-scattered starlight, and the photon field of Westerlund 1 itself. The effect of the magnetic field on cooling and propagation is discussed. Klein-Nishina effects are found to be important in determining the spectral evolution of the electron population. Results: A leptonic origin of the bulk of the observed $\gamma$-rays is preferable. The model is energetically plausible, consistent with the presence of a strong shock, and allows for the observed energy-independent morphology. The hadronic model faces two main issues: confinement of particles to the emission region and an unrealistic energy requirement.


Introduction
The standard paradigm of cosmic ray (CR) acceleration in the Galaxy states that supernova remnants (SNR) are the dominant source class at ∼GeV-PeV. This scenario is faced with several long-standing problems (for a review see Gabici et al. 2019), for example the mis-match of models to the observed 22 Ne/ 20 Ne ratio and the fact that acceleration to PeV energies is only conceivable under certain conditions, namely high shock velocities in dense environments. In addition, the γ-ray spectra of several SNR have been found to cut off early in the TeV-band, disfavouring recent PeV acceleration (see Funk 2015). Superbubbles (SBs) forming around young star clusters and associations present an alternative scenario, as was recognised early on (Cesarsky & Montmerle 1983). Young star clusters typically contain hundreds of massive stars with powerful, supersonic winds (see Portegies Zwart et al. 2010). These winds blow a cavity in the remnant of the parent molecular cloud. If the cluster is sufficiently compact, a termination shock forms inside this cavity. Beyond the termination shock radius lies the hot, shocked interior of the SB, which is delimited at the outer edge by a contact discontinuity separating the wind from a thin shell of sweptup material (Weaver et al. 1977;Mac Low & McCray 1988). This picture applies after a brief initial phase when the shell becomes radiative. While observations have confirmed the early theory in broad strokes, there remains a systematic discrepancy in bubble size and temperature, which is known as the SB ene-mail: lucia.haerer@mpi-hd.mpg.de ergy crisis (see Oey 2009;Vieu 2021, Ch. 1.5). Losses due to non-negligible interstellar pressure, a porous and thermally conducting shell, dust emission, and escaping non-thermal particles are suggested to reduce the energy available to inflate and heat up the SB. The analytic theory can be corrected by introducing an empirical pre-factor, ξ b , scaling the energy input by the wind. Vieu et al. (2022a) estimate ξ b ≈ 22% from observations, a value that is in the range predicted by current simulations (e.g. Gupta et al. 2018b).
Several models explore the SB scenario for galactic CRs (e.g. Bykov 2001;Ferrand & Marcowith 2010;Morlino et al. 2021;Vieu et al. 2022a), discussing multiple sites with favourable conditions for particle acceleration: the cluster wind termination shock, the turbulent bubble interior, and the cluster itself. The termination shock can reach high Mach numbers and potentially accelerate particles to PeV energies (see recent works by Morlino et al. 2021;Vieu et al. 2022b) due to its size (∼10s of pc, Weaver et al. 1977) and the high velocities expected for cluster winds (∼2000-3000 km s −1 , Stevens & Hartwell 2003). In the turbulent bubble interior, particles could then be stochastically reaccelerated by a Fermi-II-type process (Klepach et al. 2000;Bykov 2001). Particles could also be accelerated by supernovae (SNe) exploding in the cluster and expanding inside the lowdensity wind, therefore reaching higher shock velocities, or at colliding winds inside the cluster (see, e.g. discussion in Vieu et al. 2022a).
In the last decade, evidence for particle acceleration in star cluster environments has accumulated through γ-ray observa-Article number, page 1 of 8 arXiv:2301.10496v1 [astro-ph.HE] 25 Jan 2023 A&A proofs: manuscript no. main tions (Ackermann et al. 2011;Abramowski et al. 2012Abramowski et al. , 2015Yang et al. 2018;Aharonian et al. 2019;Abeysekara et al. 2021;Mestre et al. 2021). In this work, we examine the TeV γ-ray observations of the region around Westerlund 1 taken with the High Energy Stereoscopic System (H.E.S.S.) and reported in Aharonian et al. (2022) (henceforth HC22). Westerlund 1 is a young (4-5 Myr, Beasor et al. 2021), massive (3-5 × 10 4 M , Brandner et al. 2008, see also Portegies Zwart et al. 2010and Lim et al. 2013, compact (half-mass radius ∼1 pc) star cluster, located at a distance of ∼4 kpc from Earth and ∼4.6 kpc from the Galactic Center (GC) (Kothes & Dougherty 2007;Davies & Beasor 2019;Negueruela et al. 2022). The cluster contains a large collection of young massive stars, including 24 Wolf-Rayet (WR) stars 1 , one Luminous Blue Variable (LBV), 10 Yellow Hypergiants (YHGs) and Red Supergiants (RSGs), and several bright OB supergiants (Clark et al. 2020). Extended TeV γ-ray emission from the vicinity of the cluster was first reported by the H.E.S.S. collaboration in Abramowski et al. (2012) and designated HESS J1646−458. A Fermi analysis by Ohm et al. (2013) revealed a GeV source, slightly off-set from the TeV emission and the cluster. Muno et al. (2006) detected diffuse, non-thermal X-rays extending at least 5' outwards from the cluster. The 2022 H.E.S.S. results (HC22) reveal an emission region ∼1.5 • in diameter which is centred just slightly off the cluster position (see Fig. 1). The emission has a ring-like, energy-independent morphology, and a γ-ray spectrum with spectral index ∼2.4 up to ∼80 TeV, which is constant across the source within the range of uncertainty. The total γ-ray luminosity is 9×10 34 erg s −1 between 0.37 and 100 TeV, adopting a distance of 3.9 kpc. HC22 conclude that particle acceleration in the Westerlund 1 SB or inside the cluster itself is the most promising interpretation of the results. This idea is supported by work based on the first H.E.S.S. publication on the source (e.g. Aharonian et al. 2019). The region harbours other source candidates, but they do not provide sufficient power or cannot account for the extent of the emission (HC22). Here, we investigate the scenario of particle acceleration at the cluster wind termination shock in the light of the 2022 H.E.S.S. results. In Sect. 2, we constrain properties of the SB and the cluster wind and outline basic requirements for shock acceleration. In Sect. 3 we characterise the morphology expected from proton and electron cooling and transport. Section 4 discusses model spectra. Section 5 then combines the findings to a full picture.

Characterisation of the Westerlund 1 region
In this section, we characterise the Westerlund 1 region and estimate key parameters. We discuss the substantial uncertainty afflicting many parameters and select fiducial values for the following analysis, which are summarised in Table 1.

The superbubble
The size, termination shock position, and interior density of the SB are largely determined by cluster wind characteristics, the density of the environment, which we term external density, and the age of the cluster. A key parameter is the cluster wind power, which is defined as L w = 0.5Ṁv 2 w , with the mass-loss rate,Ṁ, and the wind speed, v w . L w is set by the sum of the contributions from stellar winds and SNe, which we discuss in turn. The WR wind power can serve as a lower bound on the former, neglecting winds from other star classes. At a rather low mass of 2 µG magnetic field (emission region) B acc 2 µG magnetic field (acceleration region) n ext 100 cm −3 external density, outside SB L bol 10 41 erg s −1 cluster bolometric luminosity T eff 40,000 K cluster effective temperature ξ b 22% empirical L w scaling factor 10-15 M , a WR star is expected to have an average wind power of L w = 10 37 erg s −1 (Seo et al. 2018). With its 24 WR stars, Westerlund 1 therefore has L w 2.4 × 10 38 erg s −1 . For a more accurate estimate, we calculate the time-dependent wind power of a toy-model cluster. To accommodate uncertainty we vary the key parameters around a range of plausible values. We first populate the cluster according to an initial mass function (IMF), ξ(m) ∝ m −Γ , assuming a cluster mass of M = 3-5 × 10 4 M and lower and upper bounds for stellar masses of 0.4 M and 120 M . As Westerlund 1 has a top-heavy IMF (Lim et al. 2013), Γ is set to 1.8-2.0. Then, we let the cluster evolve by computing the mass-loss in each time-step and updating stellar masses accordingly. In addition, stars that surpass their lifetime according to Limongi & Chieffi (2006) are removed or enter a WR phase if M star > 20 M . Wind power, mass-loss, and WR phase lifetime are taken from Seo et al. (2018). We adjust the WR population in the model cluster to that of Westerlund 1 at the time the wind power is evaluated. The resulting cluster wind power at 4-5 Myr is 1.6-2.8 × 10 39 erg s −1 . The same approach can be used to estimate the average power hitherto ejected by SNe, which yields 3.7 × 10 38 -1.5 × 10 39 erg s −1 , assuming 10 51 erg per SN. Taking the WR power discussed in the beginning of this section as a lower limit for the wind power, the estimate for the total power is ∼0.5-4 × 10 39 erg s −1 . The γ-ray luminosity cited in Sect. 1 amounts to ∼0.002-0.02% of this estimate, a broadly plausible value for the combined efficiencies of particle acceleration and γ-ray emission. For the following analysis, we conservatively select L w = 10 39 erg s −1 as a fiducial value. For the mass-loss, we again obtain a lower estimate from Seo et al. (2018). A WR star of 10 M loses 2 × 10 −5 M yr −1 , which results in a cluster mass-loss ofṀ ∼ 5 × 10 −4 M yr −1 . As WR have mass-loss rates greatly exceeding those of main sequence stars, we assume they dominate the total mass-loss and takeṀ cl = 5-10 × 10 −4 M yr −1 . For the cluster wind velocity, we obtain which means that v w takes values of ∼1300-5000 km s −1 for the ranges of L w andṀ given above. Note that the lower and upper bound represent quite extreme edge cases. Typical v w are 2000-3000 km s −1 (e.g. Stevens & Hartwell 2003).
In addition to the characteristics of the wind, the external density, n ext , is needed to determine the extent of the SB. Neutral gas tracers reveal an average density of 10.5 cm −3 for H 2 , as traced by CO and 3.2 cm −3 for HI, from 21 cm observations, assuming a distance of 3.9 kpc (HC22). However, the distribution of material is not homogeneous and can be as high as 190 cm −3 , as measured in one CO cloud. In addition, a significant fraction of the material is likely to be ionised due to the young cluster. This material and also dust are not seen by the tracers discussed above. We therefore take n ext = 100 cm −3 as a fiducial value for the following analysis and 10 cm −3 as a lower bound. This range gives a density in the bubble interior of n int = 0.02-0.10 cm −3 (Mac Low & McCray 1988). In the default scenario (Table 1), n int = 0.078 cm −3 . Note that n int weakly depends on ξ b (n int ∝ ξ 6/35 b ), a scaling factor which encompasses losses due to the SB energy crisis (see Sect. 1).
The estimates for wind power, mass-loss, and external density let us calculate the position of the forward shock, R b , and the termination shock, R ts , which are according to analytic SB theory (Weaver et al. 1977), where L w was rescaled by ξ b (see Vieu et al. 2022a). At 3.9 kpc, R ts = 20 pc corresponds to an angular radius of ∼0.3 • , which coincides well with the inner boundary of the γ-ray ring (see HC22 and qualitatively Fig. 1). Allowing for the uncertainty in the parameters discussed earlier yields R ts ∼ 20-60 pc and R b ∼ 60-180 pc. R ts = 60 pc would place most of the H.E.S.S. emission upstream of the termination shock, which is implausible in a scenario where acceleration mainly occurs at the termination shock, leading us to disfavour the combination of high wind power and low external density.
Finally, note that the H.E.S.S. emission is slightly asymmetric and off-set from the cluster. Clumps and density inhomogeneities in the surrounding medium are a plausible origin of these features, especially as the elongation is perpendicular to the galactic plane (see Fig. 1). Asymmetry is a common feature of SBs, due to the inherent clumpiness of molecular clouds (see Chu 2008).

The magnetic field in the acceleration region
For particle acceleration at the termination shock to power the HESS J1646−458 emission, the shock must, in addition to being super-Alfvénic, have a magnetic field in the acceleration region, B acc , that allows for the observed maximum energy. The H.E.S.S. spectrum lacks a clear cut-off and extends up to ∼80 TeV, from which we infer the parent electron and proton energies to be E e max 100 TeV and E p max 480 TeV. With these values, Hillas' limit (Hillas 1984) constrains the upstream B acc to number, M A , is the ratio of the wind velocity to the Alfvén speed, v A = B acc / 4πρ(R ts ). We relate the gas density ρ to the massloss, 4πρ(R)R 2 v w =Ṁ, where R is the radial distance from the cluster, to obtain Upstream of the shock, we require (6) where M A = 10 was chosen as a reference value for a strong shock. Jointly with Eq. 4, Eq. 6 places a strong constraint on B acc for the standard parameters that are given, especially for protons as they require a higher E max . The parameter space is especially narrow if the wind is slow. We select B acc = 2 µG for the following analysis.

The cluster photon field
Westerlund 1 contains a large sample of luminous stars, as discussed in the introduction. To obtain a lower bound for the bolometric luminosity of the cluster, we sum those of the brightest stars, for which measurements are available 2 , which yields L bol > 5 × 10 40 erg s −1 . The sample only contains a subset of the brightest stars, we therefore conservatively adopt an estimate L bol ∼ 10 41 erg s −1 . Comparing this value to our estimate for the wind power, we note that L bol ≈ 100L w presents a good approximation. The energy density of the cluster photons is U cl = L bol (4πr 2 c) −1 ∼ 42(L bol /10 41 erg s −1 )(R ts /20 pc) −2 eV cm −3 . In the following, U cl is always evaluated at the termination shock and therefore serves as an upper bound to the energy density inside the SB, at R > R ts . We model the emission from the cluster as thermal emission. In our sample, the YHGs are as luminous as the WR stars and the OB super-giants combined, which makes it difficult to draw a clear conclusion on whether hot or cool stars dominate the emission. We therefore investigate a broad range of effective temperatures for the cluster, T eff = 10,000-50,000 K. The threshold electron energy above which Klein-Nishina (KN) suppression is important at these temperatures is E KN = 75(T/10 4 K) −1 GeV (e.g. Blumenthal & Gould 1970). Thus, despite its large energy density, the cluster photon field does not result in a cooling time of the order of a few kyr, as would be expected for TeV particles in the Thomson regime.

Morphology
Two notable features of the source are its ring-like shape and energy-independent morphology. In the previous section, we showed that the position of the termination shock correlates with the size of the γ-ray ring. The morphology and extent of the emission region is determined by the transport and cooling of particles in the SB interior, which we investigate in this section. We also discuss the requirements to obtain an energy independent morphology.

Cooling and transport -Protons
The cooling time for protons through proton-proton (PP) interactions in the interior of the bubble can be written as (e.g. Aharonian 2004) where we use the cross-section, σ pp , of Kafexhiu et al. (2014), normalised such that the scaling factor x = 1 at 10 TeV. The value of x varies by less than 50% in the TeV-PeV range. Clearly, t cool t sys , as n int ≤ 0.1 cm −3 (see Sect. 2). This result has two important consequences: firstly, the efficiency of γ-ray production via PP interactions inside SBs is low, of the order t sys /t cool , which increases the power requirement. Secondly, the size of the emission region scales with the age of the system. We discuss the consequences of advection and diffusion in turn. A common result of SB models is that the density is approximately constant inside the SB. Mass continuity requires vρR 2 = const., such that ρ(R) = const. ⇒ v =Ṙ ∝ R −2 . The flow velocity immediately downstream of the termination shock is v(R ts ) ≈ v w /4 in the strong shock limit. The time taken to advect from the shock to the inferred gamma-ray source radius R source ≈ 60 pc is obtained by solving the above differential equation: which is much shorter than the age of the system. In fact, protons are expected to fill the interior of the SB completely. This finding is inconsistent with the size of HESS J1646−458, unless the wind is exceptionally weak, such that R b = R source ≈ 60 pc. In addition, the Bohm limit places a lower bound on the diffusive transport, assuming isotropy. We parameterise the diffusion coefficient as where we introduce the Hall parameter, h(E), as the ratio of a particle's mean free path to its gyroradius, h = λ/r g ≥ 1, with r g = E(eB) −1 . The diffusion length is p = 4Dt sys , which, in the Bohm limit (h = 1), can be written as where B is the magnetic field in the emission region. From Eq. 10 and R ts 20 pc, it follows that protons could fill the SB interior by diffusion alone, even for high B such as 10 µG. Considering that diffusion and advection are both at play and that particles can diffuse across the contact discontinuity, which extends the emission region beyond R b , these findings effectively rule out that the observed morphology is dominated by protons. The effect of radial diffusion would have to be essentially reduced to zero, in addition to the presence of a weak wind causing the bubble to be small. Furthermore, such a scenario would raise the question why the emission region does not coincide with the expected location of the shell, despite its large density (∼10 3 n int , see Gupta et al. 2018a). The scenario that the observed emission is shell emission, would require R b ≈ 20 pc, which is inconsistent with SB theory and observations (see Sect. 2.1 and Chu 2008).

Cooling and transport -Electrons
Electrons mainly cool by interacting with ambient photon fields and the magnetic field. We consider direct and dust-scattered starlight (Popescu et al. 2017), the Cosmic Microwave Background (CMB, 2.7 K, 0.25 eV cm −3 ) and the cluster photon field as discussed in Sect. 2.3. The model by Popescu et al. (2017) is axisymmetric with respect to the GC and assumes a distance GC-Earth of 8 kpc. We calculate the cooling times using the GAMERA library 3 , taking into account inverse Compton (IC) scattering on the above mentioned photon fields, synchrotron radiation, and bremsstrahlung. The results are shown in Fig. 2, on the left hand-side for a set of parameters that we adopt as default case in the following (see Table 1) and on the right hand-side for limiting cases of maximal and minimal cooling (see caption of Fig. 2). Several conclusions can be drawn from the figure. First, the cooling break is below 100 GeV, even in the case of minimal cooling. TeV particles are therefore expected to display a cooled spectrum. Second, cooling at 10-100 TeV is dominated by IC scattering on the CMB and synchrotron radiation. Even though the cluster photon field has the highest energy density, it is strongly affected by the KN suppression in this energy range. The total cooling time at 10 TeV is 55.1 kyr for default parameters. The diffusion length is given by e = 11 D 1.66 × 10 26 cm 2 s −1 0.5 t cool 55.1 kyr 0.5 pc , where D is the diffusion coefficient at the Bohm limit with B = 2 µG. The magnetic field carried by the wind from the cluster is expected to be disordered on scales R ts . Thus particles with gyro-radius equal to the outer-scale of the turbulence will have h ≈ 1, but will have an energy scaling that reflects the turbulence 10 3 10 1 10 1 10 3 Energy [TeV] 10 0 10 1 10 2 10 3 10 5 E 1 E 1/2 E 1/3 ALL Fig. 2. Left: electron cooling times for the default case (see Table 1). The broad grey line is the sum of all the components shown in colour. For a description of the photon fields for IC scattering see the text. Right: range of plausible cooling times for synchrotron radiation (orange), bremsstrahlung (green), and the sum of all IC components (blue), resulting from a cluster photon field energy density of U cl = 3-180 eV cm −3 , a density inside the superbubble of n int = 0.02-0.1 cm −3 , an effective cluster temperature of T eff = 10,000-50,000 K, a magnetic field of B = 1-10µG, and an enhancement of the diffuse stellar and dust emission by a factor 1-3. The enhancement of the diffuse component is motivated by the proximity to the cluster and increased dust density in the region, compared to the standard ISM. The inset shows the behaviour in the TeV-band in more detail highlighting the loss-time scaling with energy. The default case (broad grey line) shows a behaviour that would require diffusion close to the Kraichnan regime (D ∼ E 1/2 ) to reproduce energy-independent morphology.
where r inj is the injection scale of the turbulence, which can be assumed to be of the order of the average distance between stars in the cluster. The turbulence scaling inside SB is an open question, which is difficult to tackle with numerical simulations. Existing results are broadly consistent with δ 0.5 (see Vieu 2021).  Fig. 3 in HC22), which we discuss in detail in Sect. 3.3. The cooling time for 10 TeV electrons is ≈ 70 kyr. In this time electrons are advected to a radius of ∼37 pc approximately half way across the emission region. Diffusive transport, while uncertain, can account for the source size for range of values for δ. The situation changes if the photon fields and the magnetic field deviate from our standard assumed values. With maximal IC cooling (see Fig. 2, right panel) the total cooling time at 10 TeV is reduced to 28.7 kyr, reducing the advection length by ∼6 pc. This can be compensated if diffusion is closer to the Kolmogorov regime. In the case of larger magnetic field strength, both the cooling time and diffusion coefficient will reduce. For B 5 µG, the transport should be diffusion dominated with δ 1/2. B = 10 µG requires δ ≈ 1/3.

Energy independent morphology
HC22 show H.E.S.S. maps for E > 0.37, >1, and >4.9 TeV in their Fig. 3. Only a negligible change in source morphology is present between these bands. In this section, we discuss the implications of this energy independence for our model. For electrons, we predict the source size in each map based on the map's mean logarithmic γ-ray energy, where the mean is weighted with the log dN/dE flux. Binned data from the combined energy spectrum are used (see Fig. 7 in HC22). The averages are 4.9, 8.3, and 20 TeV for maps (a) E > 0.37, (b) >1, and (d) >4.9 TeV, respectively. Next, we determine the cooling time of the parent electrons in the ambient photon fields. Electron energies are obtained by calculating IC spectra from several narrow-band electron injection spectra and selecting the spectrum with the highest flux at the given γ-ray energy, which yields 16, 27, and 43 TeV for maps (a), (b), and (d), respectively. Figure 3 shows that the difference between map (a) and (d) is expected to be ∼5 pc, which corresponds to ∼4.4 , for advection dominated transport, which is smaller than the kernel used to smooth the maps and therefore consistent with the observation of energy independent morphology.
Energy independence is also achievable in the case of diffusion dominated transport. In the Thomson regime, t cool ∝ E −1 , for Bohm diffusion, as D ∝ E, the diffusion length is energy independent. In the transition to the KN regime, energy independent morphology can arise for other diffusion scenarios as well, for example Kolmogorov diffusion, where D ∝ E 1/3 , would have to be compensated by an E −1/3 scaling of the cooling time. The inset in Fig. 2 on the right hand-side illustrates this result. Energy independent morphology in the TeV arises in the Bohm case only if synchrotron cooling dominates and, in contrast, for Kolmogorov if IC cooling is strong and dominates. The default case requires diffusion close to Kraichnan scaling to reproduce energy independent morphology. However, due to the many uncertainties involved, no definite conclusions on the scaling should be drawn.
Finally, consider the hadronic scenario. As the energy dependence of the proton cooling time is weak (see Eq. 7) and diffusion in the radial direction has to be suppressed, an energy independent morphology is trivially expected. Note, however, that the source morphology is hard to reconcile with a hadronic scenario and it is therefore disfavoured (see Sect. 3.1).

Spectrum
We construct both a hadronic (PP interaction) and a leptonic (IC emission) model of the TeV γ-ray spectrum, using the GAMERA library to accurately include KN corrections. The particle injection spectra are set to be power laws with an exponential cut-off, dN/(dEdt) ∝ E −α inj exp(−E/E cutoff ). The normalisation is obtained assuming that the power continuously injected into the particle spectra in the range 1 GeV-1 PeV is η p L w and η e L w for protons and electrons, respectively. In other words, the efficiencies η p and η e give the fraction of wind power converted into CR power at the termination shock. From the injection spectra, GAMERA is used to calculate the cooled spectra and consequentially the γ-ray emission. We use the GEANT 4.10.0 hadronic emission model and standard parameters for the leptonic emission. We consider the ambient radiation fields discussed in Sect. 3.2: the cluster photon field, diffuse starlight, dust-scattered starlight, and the CMB. The best-fitting models are summarised in Fig. 4.

Leptonic emission
In the electron spectrum, E cutoff is set to the energy where the cooling time equals the acceleration time, which assuming approximately equal upstream and downstream residence times is (Drury 1983) (13) Figure 5 shows the cooling and accelerations times for B acc = 1-10 µG in the Bohm limit. Note that for h > 1, the acceleration time increases compared to Fig. 5, lowering E cutoff . For the cooling, we consider the default case shown in Table 1 and in Fig. 2 (left) and vary B acc . A 2 µG field places the cut-off in the particle spectrum at ∼170 TeV. The resulting IC spectra (Fig. 4 Fig. 4. The best models for the spectrum of HESS J1646−458 (data taken from HC22). Standard parameters are assumed (see Table 1). The injected particle spectrum is a power law with an exponential cut-off, where α inj is the index and E cutoff the cut-off energy. The GAMERA library is used to calculate the cooled particle spectrum and γ-ray production. cool. times acc. times 1 µG 2 µG 5 µG 10 µG Fig. 5. Cooling and acceleration times (solid and dashed lines, respectively) for a magnetic field in the acceleration region B acc = 1-10 µG. The cut-off in the electron spectrum is to set the energy, where the two timescales are equal. Photon fields are chosen as in Fig. 2. lines) are consistent with the H.E.S.S. data and require efficiencies of only 0.09-0.28% for α inj = 2.1-2.3 and default parameters (see Table 1). For larger B acc , the IC model under-predicts the data at 10 TeV. We note that the magnetic field in the acceleration region B acc need not be the same as the average value in the post-shock cooling region. Figure 6 illustrates the dependence of the IC spectrum on the photon fields. The left-hand panel concerns the cluster and diffuse photon fields. A low effective temperature steepens the spectrum between ∼0.1-10 TeV and results in a higher flux, because the photons are less affected by KN suppression due to their lower energy. Such models are generally easier to reconcile with the slight steepening observed in the spectrum at 1 TeV. For the diffuse fields (right-hand panel), we enhance the values predicted by Popescu et al. (2017) by a factor 2-3. This accounts for reprocessed and lower energy cluster photons. An enhancement of the diffuse starlight and dust emission smooths the feature introduced by the cluster photon field at ∼10-100 GeV, flattening the spectrum. In summary, IC emission explains the spec- trum of HESS J1646−458 well, although the required balance between cooling and acceleration times is a strong constraint on the model.
The energetic electrons are also expected to produce synchrotron emission (see Fig. 7). HC22 deduced the brightness at 30 GHz in a region with 1 • radius around Westerlund 1 using data from the Planck satellite, which can serve as an upper bound on the synchrotron component. Figure 7 shows that our models for B = 1-10 µG are consistent with this bound.

Hadronic emission
Putting aside the above mentioned limitations on morphology and advection, we consider a single-zone model for hadronic emission. The spectral shape can be reproduced with γ-ray emission from PP interaction using α inj = 2.2-2.4 and E cutoff = 400 TeV-1.2 PeV, although the required efficiency is an obstacle for the hadronic model. The standard setting, L w = 10 39 erg s −1 and n ext = 100 cm −3 , requires η p = 26-92%. These values are optimistic, especially because we do not consider particle escape and the best-fit model corresponds to the upper bound of this range (α inj = 2.4, E cutoff = 1.2 PeV). Harder injection spectra have trouble accounting for the steepening at 1 TeV. The efficiency depends inversely on L w and n int . An increase of their product by a factor of 4 brings η p down to 6.5-23%, which is more plausible for shock acceleration. We noted in Sect. 3 that the size of the source is considerably over-predicted in the hadronic scenario unless R b 60 pc. From this constraint and the efficiency requirement set above, it follows that n ext 300 cm −3 (cf. Eq. 2 and 3), which is much higher than typical values (see Sect. 2.1). An alternative explanation is the presence of a spectral break not far below 1 TeV.
The above consideration concerns the case in which the γ-ray emission originates in the bubble interior. An estimate for the γray luminosity expected from the SB shell, L sh , can be obtained from the ratio of the residence time in the shell, t res , to the proton 10 7 10 8 10 9 10 10 10 11 10 12 10 13 10 14 frequency [Hz] 10 8 10 6 10 4 10 2 10 0 Energy [eV] 10 0 10 1 10 2 10 3 Jydeg 2 Planck limit (HC22) B-field 1 G 2 G 5 G 10 G Fig. 7. Predictions for the synchrotron brightness for an emission region magnetic field B = 1-10 µG. The 2 µG model (black) uses the same set of parameters as the blue dotted model in Fig. 4 and was used to set the normalisation for all models. For reference, the source covers ∼1 deg 2 .
cooling time (see Eq. 7). This yields where we have assumed that 1% of L w is converted to protons above 10 TeV. This flux is more than an order of magnitude below that seen by H.E.S.S. (see Sect. 1). Since the shell emission is expected to be located at radii 60 pc and therefore is spread over a large solid angle, a non-detection in H.E.S.S. is not surprising. Future observations may reveal the shell component, which would greatly enhance our understanding of particle acceleration in the region. Note, however, that the numerical values used in Eq. 14 are uncertain.

Discussion and conclusion
We investigated the scenario of particle acceleration at a strong termination shock in a superbubble (SB) surrounding Westerlund 1, with regard to the H.E.S.S. observations reported in 2022 (HC22). In Sect. 3 and 4 we discussed how the TeV γ-ray morphology and spectrum of the source can be modelled with either inverse Compton (IC) emission or the decay of neutral pions produced via proton-proton (PP) interactions.
The hadronic interpretation faces two main difficulties. Protons are expected to be uncooled for the given cluster age of 4 Myr. Transport by both advection and diffusion over this timescale significantly overpredicts the size of the emission region. The only conceivable, though unrealistic, solution is a SB radius of 60 pc constraining the advection length and an essentially complete suppression of radial diffusion. The second issue is the energetics: 26-92% of the cluster wind power have to be converted into γ-ray luminosity in the standard scenario. Alternative scenarios require unrealistic densities of the external medium of 300 cm −3 , considering that morphology has to be accounted for as well. The efficiency issue is especially relevant for models with steep injection spectra, as they require more power at TeV energies than models with power law indices around 2. Figure 4 shows that cut-off energies >1 PeV require such steep injection spectra. A potential solution to the efficiency issue is a break not far below the H.E.S.S. band. However, considering the general picture and that HC22 did not find an alignment between the γ-ray emission and the neutral gas distribution, we deem a purely hadronic scenario, especially one where the particles reach PeV energies, unlikely.
The leptonic model provides an overall consistent picture that is in good agreement with the available data. We investigate IC scattering on diffuse starlight and dust-scattered starlight, the CMB, and a T = 10,000-50,000 K thermal cluster photon field. Cooling shortens the advection and diffusion lengths considerably compared to the hadronic case. As a result, the model is consistent with both the radius of HESS J1646−458 and the energy independent morphology. We note that the mean γ-ray energy of the H.E.S.S. maps (Fig. 3 in HC22) only varies between 4.9 and 20 TeV. The predicted energy dependence of the transport in this range is weak and not detectable given the resolution of the maps. Due to continuous injection by the star cluster and the Klein-Nishina (KN) suppression, which causes the cooling times of high energy electrons to rise, a sufficient supply of high energy electrons is present in the bubble, contrary to what was reported in Bhadra et al. (2022). In fact, we find that the efficiency required for the conversion of cluster wind power to IC luminosity is <0.3%. The spectral shape is quite well predicted for an acceleration region magnetic field of B acc 2 µG. A higher B acc lowers the cut-off energy and hence underpredicts the flux at 10 TeV. This constraint is the largest caveat of the leptonic model, but does not threaten its validity as B acc ∼ 1-2 µG fulfils Hillas' limit and our requirement of a strong termination shock. If B acc 2 µG, a hard hadronic component could compensate for the early cut-off in the electron spectrum. Such a joint model is less challenging than a purely hadronic model in terms of efficiency, but requires electrons and protons to be constrained within the same emission region which presents an additional assumption. The spectral shape of HESS J1646−458 is well predicted by such models. Though the leptonic model is preferred globally, hadronic emission is expected at some level, and identifying its presence is a necessary step in determining the contribution of stellar clusters to the Galactic CR population.
A key open question is the position of the cut-off in the γ-ray spectrum, which southern hemisphere γ-ray observatories such as CTA and SWGO will be able to address. Radio data could constrain the magnetic field in the emission region, improving upon the upper bound set by HC22 (see Fig. 7). The GeV-band can also provide valuable insight, as many of the models shown in Fig. 4 separate in this range. The IC emission from the cluster photon field also peaks in the GeV (see Fig. 6) and is key in understanding the spectral behaviour.