© The Authors 2023

1. 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 mismatch of models to the observed ²²Ne/²⁰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 swept-up 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 energy crisis (see Oey et al. 2009; Vieu 2021, Chap. 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 (∼10 s of pc, Weaver et al. 1977) and the high velocities expected for cluster winds (∼2000 − 3000 km s⁻¹, 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 low-density wind, therefore reaching higher shock velocities, or at colliding winds inside the cluster (see, e.g. the discussion in Vieu et al. 2022a).

In the past decade, evidence for particle acceleration in star cluster environments has accumulated through γ-ray observations (Ackermann et al. 2011; Abramowski et al. 2012, 2015; Yang 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⁴ M_⊙; Brandner et al. 2008, see also Portegies Zwart et al. 2010; 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 centre (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¹, one luminous blue variable (LBV), ten 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 that 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³⁴ erg s⁻¹ between 0.37 and 100 TeV, adopting a distance of 3.9 kpc. HC22 concluded 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.

Fig. 1. Map of HESS J1646−458 taken from HC22, overlaid with the termination shock position predicted for the parameters indicated in Eq. (2). The star indicates the position of Westerlund 1. The dashed grey line traces the Galactic plane.

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

2.1. 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_{\mathrm{w}}=0.5\dot{M}{v}_{\mathrm{w}}^2$, 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 10−15 M_⊙, a WR star is expected to have an average wind power of L_w = 10³⁷ erg s⁻¹ (Seo et al. 2018). With its 24 WR stars, Westerlund 1 therefore has L_w ≳ 2.4 × 10³⁸ erg s⁻¹. For a more accurate estimate, we calculated the time-dependent wind power of a toy-model cluster. To accommodate uncertainty, we varied the key parameters around a range of plausible values. We first populated the cluster according to an initial mass function (IMF), ξ(m)∝m^−Γ, assuming a cluster mass of M = 3 − 5 × 10⁴ 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), Γ was 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) were removed or entered a WR phase if M_star > 20 M_⊙. Wind power, mass loss, and WR phase lifetime were taken from Seo et al. (2018). We adjusted 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³⁹ erg s⁻¹. The same approach can be used to estimate the average power hitherto ejected by SNe, which yields 3.7 × 10³⁸ − 1.5 × 10³⁹ erg s⁻¹, assuming 10⁵¹ erg per SN. When the WR power discussed in the beginning of this section is taken as a lower limit for the wind power, the estimate for the total power is ∼0.5 − 4 × 10³⁹ erg s⁻¹. 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³⁹ erg s⁻¹ 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⁻⁵ M_⊙ yr⁻¹, which results in a cluster mass-loss of Ṁ ∼ 5 × 10⁻⁴ M_⊙ yr⁻¹. As the mass-loss rates of WR stars greatly exceed those of main-sequence stars, we assumed that they dominate the total mass loss and took Ṁ_cl = 5 − 10 × 10⁻⁴ M_⊙ yr⁻¹. For the cluster wind velocity, we obtain

$$\begin{array}{c}\hfill v_{\mathrm{w}}\sim 2500{\left(\frac{L_{\mathrm{w}}}{{10}^{39}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{0.5}{\left(\frac{\dot{M}}{5\times {10}^{-4}{M}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{-0.5}\mathrm{km}{\mathrm{s}}^{-1},\end{array}$$(1)

which means that v_w takes values of ∼1300 − 5000 km s⁻¹ for the ranges of L_w and Ṁ given above. The lower and upper bound represent quite extreme edge cases. Typical v_w are 2000−3000 km s⁻¹ (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⁻³ for H₂, as traced by CO and 3.2 cm⁻³ 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⁻³, 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 took n_ext = 100 cm⁻³ as a fiducial value for the following analysis and 10 cm⁻³ as a lower bound. This range gives a density in the bubble interior of n_int = 0.02 − 0.10 cm⁻³ (Mac Low & McCray 1988). In the default scenario (Table 1), n_int = 0.078 cm⁻³. We note that n_int weakly depends on ξ_b ($n_{\mathrm{int}}\propto {\xi }_{\mathrm{b}}^{6/35}$), which is a scaling factor that encompasses losses due to the SB energy crisis (see Sect. 1).

The estimates for wind power, mass loss, and external density allow us to calculate the position of the forward shock, R_b, and the termination shock, R_ts, which are

$$\begin{array}{cc}\hfill R_{\mathrm{ts}}& =20.4{\left(\frac{{\xi }_{\mathrm{b}}}{0.22}\right)}^{-1/5}{\left(\frac{\dot{M}}{5\times {10}^{-4}{M}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{3/10}{\left(\frac{n_{\mathrm{ext}}}{100{\mathrm{cm}}^{-3}}\right)}^{-3/10}\hfill \\ \hfill & \times {\left(\frac{v_{\mathrm{w}}}{2500\mathrm{km}{\mathrm{s}}^{-1}}\right)}^{1/10}{\left(\frac{t_{\mathrm{sys}}}{4\mathrm{Myr}}\right)}^{2/5}\mathrm{pc},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill R_{\mathrm{b}}& =74.2{\left(\frac{{\xi }_{\mathrm{b}}{L}_{\mathrm{w}}}{0.22\times {10}^{39}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{1/5}{\left(\frac{n_{\mathrm{ext}}}{100{\mathrm{cm}}^{-3}}\right)}^{-1/5}\hfill \\ \hfill & \times {\left(\frac{t_{\mathrm{sys}}}{4\mathrm{Myr}}\right)}^{3/5}\mathrm{pc},\hfill \end{array}$$(3)

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). When we allow for the uncertainty in the parameters discussed earlier, this 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. We therefore disfavour the combination of high wind power and low external density.

Finally, we 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 because of the inherent clumpiness of molecular clouds (see Chu 2008).

2.2. 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_{\max }^{\mathrm{e}}\gtrsim 100$ TeV and $E_{\max }^{\mathrm{p}}\gtrsim 480$ TeV. With these values, the Hillas limit (Hillas 1984) constrains the upstream B_acc to

$$\begin{array}{c}\hfill B_{\mathrm{acc}}>\frac{E_{\max }c}{e{v}_{\mathrm{w}}{R}_{\mathrm{ts}}}=0.7-3.1{\left(\frac{v_{\mathrm{w}}}{2500\mathrm{km}{\mathrm{s}}^{-1}}\right)}^{-1}{\left(\frac{R_{\mathrm{ts}}}{20\mathrm{pc}}\right)}^{-1}\mathrm{\mu }\mathrm{G}\end{array}$$(4)

for electrons and protons, respectively. In addition, efficient particle acceleration requires M_A ≫ 1, where the Alfvénic Mach number, M_A, is the ratio of the wind velocity to the Alfvén speed, $v_{\mathrm{A}}=B_{\mathrm{acc}}/\sqrt{4\pi \rho (R_{\mathrm{ts}})}$. We relate the gas density ρ to the mass loss, 4πρ(R)R²v_w = Ṁ, where R is the radial distance from the cluster, to obtain

$$\begin{array}{c}\hfill M_{\mathrm{A}}=\frac{v_{\mathrm{w}}}{v_{\mathrm{A}}}=\frac{\sqrt{\dot{M}{v}_{\mathrm{w}}}}{B_{\mathrm{acc}}{R}_{\mathrm{ts}}}·\end{array}$$(5)

Upstream of the shock, we require

$$\begin{array}{cc}& B_{\mathrm{acc}}<4.5{\left(\frac{R_{\mathrm{ts}}}{20\mathrm{pc}}\right)}^{-1}{\left(\frac{\dot{M}}{5\times {10}^{-4}{M}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{0.5}\hfill \\ \hfill & \times {\left(\frac{M_{\mathrm{A}}}{10}\right)}^{-1}{\left(\frac{v_{\mathrm{w}}}{2500\mathrm{km}{\mathrm{s}}^{-1}}\right)}^{0.5}\mathrm{\mu }\mathrm{G},\hfill \end{array}$$(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.

2.3. 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 summed those of the brightest stars for which measurements are available², which yielded L_bol > 5 × 10⁴⁰ erg s⁻¹. The sample only contains a subset of the brightest stars. We therefore conservatively adopted an estimate L_bol ∼ 10⁴¹ erg s⁻¹. 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²c)⁻¹ ∼ 42(L_bol/10⁴¹ erg s⁻¹)(R_ts/20 pc)⁻² eV cm⁻³. 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 modelled 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 investigated 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⁴ K)⁻¹ GeV (e.g. Blumenthal & Gould 1970). Thus, despite its high energy density, the cluster photon field does not result in a cooling time of the order of a few thousand years, as would be expected for TeV particles in the Thomson regime.

3. 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 for obtaining an energy-independent morphology.

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

$$\begin{array}{c}\hfill t_{\mathrm{cool}}={(0.5n_{\mathrm{int}}{\sigma }_{\mathrm{pp}}c)}^{-1}\approx 5.2\times {10}^{8}x{\left(\frac{n_{\mathrm{int}}}{0.1{\mathrm{cm}}^{-3}}\right)}^{-1}\mathrm{yr},\end{array}$$(7)

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⁻³ (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² = const., such that ρ(R) = const. ⇒ v = Ṙ ∝ R⁻². 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 γ-ray source radius R_source ≈ 60 pc is obtained by solving the differential equation

$$\begin{array}{cc}\hfill t_{\mathrm{adv}}=& 271{\left(\frac{v_{\mathrm{w}}}{2500\mathrm{km}{\mathrm{s}}^{-1}}\right)}^{-1}{\left(\frac{R_{\mathrm{ts}}}{20\mathrm{pc}}\right)}^{-2}\hfill \\ \hfill & \times [{\left(\frac{R_{\mathrm{source}}}{60\mathrm{pc}}\right)}^3-{\left(\frac{R_{\mathrm{ts}}}{20\mathrm{pc}}\right)}^3]\mathrm{kyr},\hfill \end{array}$$(8)

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 parametrise the diffusion coefficient as

$$\begin{array}{c}\hfill D=\frac{1}{3}h(E)r_{\mathrm{g}}c,\end{array}$$(9)

where we introduce the Hall parameter, h(E), as the ratio of the mean free path of a particle to its gyroradius, h = λ/r_g ≥ 1, with r_g = E(eB)⁻¹. The diffusion length is ${\ell }_{\mathrm{p}}=\sqrt{4D{t}_{\mathrm{sys}}}$, which in the Bohm limit (h = 1) can be written as

$$\begin{array}{c}\hfill {\ell }_{\mathrm{p}}=94{\left(\frac{E}{10\mathrm{TeV}}\right)}^{0.5}{\left(\frac{t_{\mathrm{sys}}}{4\mathrm{Myr}}\right)}^{0.5}{\left(\frac{B}{2\mathrm{\mu }\mathrm{G}}\right)}^{-0.5}\mathrm{pc},\end{array}$$(10)

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 essentially have to be reduced to zero, in addition to the presence of a weak wind causing the bubble to be small. Furthermore, this scenario would raise the question why the emission region does not coincide with the expected location of the shell, despite its high density (∼10³ 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).

3.2. Cooling and transport: Electrons

Electrons mainly cool by interacting with ambient photon fields and the magnetic field. We considered direct and dust-scattered starlight (Popescu et al. 2017), the cosmic microwave background (CMB; 2.7 K, 0.25 eV cm⁻³), 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 calculated the cooling times using the GAMERA library³, 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 side for a set of parameters that we adopt as default case in the following (see Table 1) and on the right side for limiting cases of maximum and minimum 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 minimum 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

$$\begin{array}{c}\hfill {\ell }_{\mathrm{e}}=11{\left(\frac{D}{1.66\times {10}^{26}{\mathrm{cm}}^2{\mathrm{s}}^{-1}}\right)}^{0.5}{\left(\frac{t_{\mathrm{cool}}}{55.1\mathrm{kyr}}\right)}^{0.5}\mathrm{pc},\end{array}$$(11)

Fig. 2. Electron cooling times. Left: 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 Ucl = 3 − 180 eV cm−3, a density inside the superbubble of nint = 0.02 − 0.1 cm−3, an effective cluster temperature of Teff = 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 ∼ E1/2) to reproduce an energy-independent morphology.

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 a gyro-radius equal to the outer scale of the turbulence will have h ≈ 1, but will have an energy scaling that reflects the turbulence spectrum. Generally, h scales as

$$\begin{array}{cc}\hfill h(E)& ={\left(\frac{r_{\mathrm{inj}}}{1\mathrm{pc}}\right)}^{1-\delta }{\left(\frac{E}{10\mathrm{TeV}}\right)}^{\delta -1}{\left(\frac{B}{2\mathrm{\mu }\mathrm{G}}\right)}^{1-\delta }\hfill \\ \hfill & \times \{\begin{array}{cc}14\hfill & \mathrm{for}\mathrm{Kraichnan},\delta =1/2\hfill \\ 33\hfill & \mathrm{for}\mathrm{Kolmogorov},\delta =1/3\hfill \end{array},\hfill \end{array}$$(12)

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 that is difficult to tackle with numerical simulations. Existing results are broadly consistent with δ ≲ 0.5 (see Vieu 2021).

Figure 3 compares diffusion and advection times to the electron cooling time at different energies, chosen to predict the size of the source in the H.E.S.S. maps with varying energy threshold (see 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 a range of values for δ. The situation changes when the photon fields and the magnetic field deviate from our standard assumed values. With maximum 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 for when diffusion is closer to the Kolmogorov regime. In the case of a 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.

Fig. 3. Diffusion (red) and advection (blue) timescales as a function of transported distance, assuming particles start at the termination shock at Rts = 20.4 pc. Timescales for three different diffusion coefficients are shown, following Eqs. (9) and (12). The grey band indicates the radius of the source, where the large range is due to the elongated shape. Cooling times for electrons of four different energies are indicated (black). In addition to a fiducial value (10 TeV), the average parent electron energies we determine for the H.E.S.S. maps (see HC22, Fig. 3) are shown. The technique used is detailed in the text.

3.3. Energy-independent morphology

HC22 showed 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 mean logarithmic γ-ray energy of the map, where the mean is weighted with the log dN/dE flux. Binned data from the combined energy spectrum were 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 determined the cooling time of the parent electrons in the ambient photon fields. Electron energies were 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 can also be achieved in the case of diffusion-dominated transport. In the Thomson regime, t_cool ∝ E⁻¹ 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 for by an E^−1/3 scaling of the cooling time. The inset in Fig. 2 on the right 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 diffusion if IC cooling is strong and dominates. The default case requires diffusion close to Kraichnan scaling to reproduce energy-independent morphology. However, because of the many uncertainties involved, no definite conclusions on the scaling should be drawn.

Finally, we considered 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. We note, however, that the source morphology is hard to reconcile with a hadronic scenario and is therefore disfavoured (see Sect. 3.1).

4. Spectrum

We constructed 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 were set to be power laws with an exponential cut-off, dN/(dEdt)∝E^−α_injexp(−E/E_cutoff). The normalisation was obtained assuming that the power continuously injected into the particle spectra in the range 1 GeV−1 PeV is η_pL_w and η_eL_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 was used to calculate the cooled spectra and consequentially the γ-ray emission. We used the GEANT 4.10.0 hadronic emission model and standard parameters for the leptonic emission. We considered 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.

Fig. 4. 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 Ecutoff the cut-off energy. The GAMERA library is used to calculate the cooled particle spectrum and γ-ray production.

4.1. Leptonic emission

In the electron spectrum, E_cutoff was set to the energy at which the cooling time equals the acceleration time, which assuming approximately equal upstream and downstream residence times is (Drury 1983)

$$\begin{array}{cc}\hfill t_{\mathrm{acc}}& \approx \frac{8D}{v_{\mathrm{w}}^2}=\frac{8h}{3}\frac{r_{\mathrm{g}}c}{v_{\mathrm{w}}^2}\hfill \\ \hfill & =6.7h\left(\frac{E}{100\mathrm{TeV}}\right){\left(\frac{B_{\mathrm{acc}}}{2\mathrm{\mu }\mathrm{G}}\right)}^{-1}{\left(\frac{v_{\mathrm{w}}}{2500{\mathrm{km}\mathrm{s}}^{-1}}\right)}^{-2}\mathrm{kyr}.\hfill \end{array}$$(13)

Figure 5 shows the cooling and accelerations times for B_acc = 1 − 10 μG in the Bohm limit. We note that for h > 1, the acceleration time increases compared to Fig. 5, lowering E_cutoff. For the cooling, we considered the default case shown in Table 1 and in Fig. 2 (left) and varied B_acc. A 2 μG field places the cut-off in the particle spectrum at ∼170 TeV. The resulting IC spectra (Fig. 4, blue 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 underpredicts 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.

Fig. 5. Cooling and acceleration times (solid and dashed lines, respectively) for a magnetic field in the acceleration region Bacc = 1 − 10 μG. The cut-off in the electron spectrum is to set the energy at which the two timescales are equal. Photon fields are chosen as in Fig. 2.

Figure 6 illustrates the dependence of the IC spectrum on the photon fields. The left 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. These models are generally easier to reconcile with the slight steepening observed in the spectrum at ≲1 TeV. For the diffuse fields (right panel), we enhanced 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 spectrum of HESS J1646−458 well, although the required balance between cooling and acceleration times is a strong constraint on the model.

Fig. 6. Effect of the variation of the cluster effective temperature (Teff, left) and an enhancement of the diffuse starlight and dust-scattered starlight radiation fields (right). Solid lines indicate the total IC spectrum. Dashed lines show only the cluster component on the left and the diffuse component on the right. Black lines show the best-fitting result using default parameters (Table 1), which is also displayed in Fig. 4 with the dotted blue line. The normalisation of all models is set equal to the default case.

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.

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 dotted blue model in Fig. 4 and was used to set the normalisation for all models. For reference, the source covers ∼1 deg2.

4.2. Hadronic emission

In addition to the limitations on morphology and advection described above, we considered 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³⁹ erg s⁻¹ and n_ext = 100 cm⁻³, requires η_p = 26 − 92%. These values are optimistic, especially because we did 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 do not account well 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 reduces η_p to 6.5 − 23%, which is more plausible for shock acceleration. We noted in Sect. 3 that the size of the source is considerably overpredicted 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⁻³ (cf. Eqs. (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.

This 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 cooling time (see Eq. (7)). This yields

$$\begin{array}{cc}\hfill L_{\gamma }^{\mathrm{sh}}& \approx {\eta }_{\mathrm{p}}{L}_{\mathrm{w}}\frac{t_{\mathrm{res}}}{t_{\mathrm{cool}}}={\eta }_{\mathrm{p}}{L}_{\mathrm{w}}\frac{R_{\mathrm{sh}}^2}{4D{t}_{\mathrm{cool}}}\hfill \\ \hfill & =3\times {10}^{33}\frac{{\eta }_{\mathrm{p}}(>10\mathrm{TeV})}{0.01}\frac{L_{\mathrm{w}}}{{10}^{39}\mathrm{erg}{\mathrm{s}}^{-1}}{\left(\frac{R_{\mathrm{sh}}}{10\mathrm{pc}}\right)}^2\hfill \\ \hfill & \times {\left(\frac{D_{\mathrm{sh}}(10\mathrm{TeV})}{{10}^{29}{\mathrm{cm}}^2{\mathrm{s}}^{-1}}\right)}^{-1}\frac{n_{\mathrm{sh}}}{200{\mathrm{cm}}^{-3}}\mathrm{erg}{\mathrm{s}}^{-1},\hfill \end{array}$$(14)

where we have assumed that 1% of L_w is converted into protons above 10 TeV. This flux is lower by more than an order of magnitude than 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. However, the numerical values used in Eq. (14) are uncertain.

5. Discussion and conclusion

We investigated the scenario of particle acceleration at a strong termination shock in a SB surrounding Westerlund 1 with regard to the H.E.S.S. observations reported in 2022 (HC22). In Sects. 3 and 4 we discussed how the TeV γ-ray morphology and spectrum of the source can be modelled with either IC emission or the decay of neutral pions produced via 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 an 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⁻³, 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 unlikely, especially one in which the particles reach PeV energies.

The leptonic model provides an overall consistent picture that agrees well with the available data. We investigated 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. Owing 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, in contrast 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 the 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. This 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. Although 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 the Cherenkov Telescope Array (Actis et al. 2011) and the Southern Wide-field Gamma-ray Observatory (Hinton & SWGO Collaboration 2022) 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 for understanding the spectral behaviour.

Acknowledgments

We thank Giada Peron, Richard Tuffs, and Felix Aharonian for helpful discussions. We also thank the referee for their comments, which helped improving the manuscript.

References

All Tables

All Figures

	Fig. 1. Map of HESS J1646−458 taken from HC22, overlaid with the termination shock position predicted for the parameters indicated in Eq. (2). The star indicates the position of Westerlund 1. The dashed grey line traces the Galactic plane.
In the text

Fig. 2. Electron cooling times. Left: 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 Ucl = 3 − 180 eV cm−3, a density inside the superbubble of nint = 0.02 − 0.1 cm−3, an effective cluster temperature of Teff = 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 ∼ E1/2) to reproduce an energy-independent morphology.

In the text

Fig. 3. Diffusion (red) and advection (blue) timescales as a function of transported distance, assuming particles start at the termination shock at Rts = 20.4 pc. Timescales for three different diffusion coefficients are shown, following Eqs. (9) and (12). The grey band indicates the radius of the source, where the large range is due to the elongated shape. Cooling times for electrons of four different energies are indicated (black). In addition to a fiducial value (10 TeV), the average parent electron energies we determine for the H.E.S.S. maps (see HC22, Fig. 3) are shown. The technique used is detailed in the text.

In the text

	Fig. 4. 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 Ecutoff the cut-off energy. The GAMERA library is used to calculate the cooled particle spectrum and γ-ray production.
In the text

	Fig. 5. Cooling and acceleration times (solid and dashed lines, respectively) for a magnetic field in the acceleration region Bacc = 1 − 10 μG. The cut-off in the electron spectrum is to set the energy at which the two timescales are equal. Photon fields are chosen as in Fig. 2.
In the text

Fig. 6. Effect of the variation of the cluster effective temperature (Teff, left) and an enhancement of the diffuse starlight and dust-scattered starlight radiation fields (right). Solid lines indicate the total IC spectrum. Dashed lines show only the cluster component on the left and the diffuse component on the right. Black lines show the best-fitting result using default parameters (Table 1), which is also displayed in Fig. 4 with the dotted blue line. The normalisation of all models is set equal to the default case.

In the text

	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 dotted blue model in Fig. 4 and was used to set the normalisation for all models. For reference, the source covers ∼1 deg2.
In the text