Open Access
Issue
A&A
Volume 671, March 2023
Article Number A4
Number of page(s) 8
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202245444
Published online 27 February 2023

© The Authors 2023

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

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 22Ne/20Ne 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−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 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 × 104M; 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) stars1, 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 × 1034 erg s−1 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.

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

Table 1.

Input parameters used in the default scenario.

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 , with the mass-loss rate, , and the wind speed, vw. Lw 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 Lw = 1037 erg s−1 (Seo et al. 2018). With its 24 WR stars, Westerlund 1 therefore has Lw ≳ 2.4 × 1038 erg s−1. 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 × 104M 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 Mstar > 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 × 1039 erg s−1. The same approach can be used to estimate the average power hitherto ejected by SNe, which yields 3.7 × 1038 − 1.5 × 1039 erg s−1, assuming 1051 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 × 1039 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 Lw = 1039 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−5M yr−1, which results in a cluster mass-loss of ∼ 5 × 10−4 M yr−1. 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−4 M yr−1. For the cluster wind velocity, we obtain

(1)

which means that vw takes values of ∼1300 − 5000 km s−1 for the ranges of Lw and given above. The lower and upper bound represent quite extreme edge cases. Typical vw are 2000−3000 km s−1 (e.g. Stevens & Hartwell 2003).

In addition to the characteristics of the wind, the external density, next, is needed to determine the extent of the SB. Neutral gas tracers reveal an average density of 10.5 cm−3 for H2, 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 took next = 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 nint = 0.02 − 0.10 cm−3 (Mac Low & McCray 1988). In the default scenario (Table 1), nint = 0.078 cm−3. We note that nint weakly depends on ξb (), 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, Rb, and the termination shock, Rts, which are

(2)

(3)

according to analytic SB theory (Weaver et al. 1977), where Lw was rescaled by ξb (see Vieu et al. 2022a). At 3.9 kpc, Rts = 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 Rts ∼ 20 − 60 pc and Rb ∼ 60 − 180 pc. Rts = 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, Bacc, 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 TeV and TeV. With these values, the Hillas limit (Hillas 1984) constrains the upstream Bacc to

(4)

for electrons and protons, respectively. In addition, efficient particle acceleration requires MA ≫ 1, where the Alfvénic Mach number, MA, is the ratio of the wind velocity to the Alfvén speed, . We relate the gas density ρ to the mass loss, 4πρ(R)R2vw = , where R is the radial distance from the cluster, to obtain

(5)

Upstream of the shock, we require

(6)

where MA = 10 was chosen as a reference value for a strong shock. Jointly with Eq. (4), Eq. (6) places a strong constraint on Bacc for the standard parameters that are given, especially for protons, as they require a higher Emax. The parameter space is especially narrow if the wind is slow. We select Bacc = 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 available2, which yielded Lbol > 5 × 1040 erg s−1. The sample only contains a subset of the brightest stars. We therefore conservatively adopted an estimate Lbol ∼ 1041 erg s−1. Comparing this value to our estimate for the wind power, we note that Lbol ≈ 100Lw presents a good approximation. The energy density of the cluster photons is Ucl = Lbol(4πr2c)−1 ∼ 42(Lbol/1041 erg s−1)(Rts/20 pc)−2 eV cm−3. In the following, Ucl is always evaluated at the termination shock and therefore serves as an upper bound to the energy density inside the SB, at R > Rts. 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, Teff = 10 000 − 50 000 K. The threshold electron energy above which Klein-Nishina (KN) suppression is important at these temperatures is EKN = 75(T/104 K)−1 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)

(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, tcool ≫ tsys, as nint ≤ 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 tsys/tcool, 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ρR2 = const., such that ρ(R) = const. ⇒ v = R−2. The flow velocity immediately downstream of the termination shock is v(Rts)≈vw/4 in the strong-shock limit. The time taken to advect from the shock to the inferred γ-ray source radius Rsource ≈ 60 pc is obtained by solving the differential equation

(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 Rb = Rsource ≈ 60 pc. In addition, the Bohm limit places a lower bound on the diffusive transport, assuming isotropy. We parametrise the diffusion coefficient as

(9)

where we introduce the Hall parameter, h(E), as the ratio of the mean free path of a particle to its gyroradius, h = λ/rg ≥ 1, with rg = E(eB)−1. The diffusion length is , which in the Bohm limit (h = 1) can be written as

(10)

where B is the magnetic field in the emission region. From Eq. (10) and Rts ≳ 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 Rb, 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 (∼103nint, see Gupta et al. 2018a). The scenario that the observed emission is shell emission would require Rb ≈ 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−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 calculated the cooling times using the GAMERA library3, 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

(11)

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

(12)

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

thumbnail 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, tcool ∝ 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 ∝ E1/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/Ecutoff). The normalisation was obtained assuming that the power continuously injected into the particle spectra in the range 1 GeV−1 PeV is ηpLw and ηeLw 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.

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

(13)

Figure 5 shows the cooling and accelerations times for Bacc = 1 − 10 μG in the Bohm limit. We note that for h > 1, the acceleration time increases compared to Fig. 5, lowering Ecutoff. For the cooling, we considered the default case shown in Table 1 and in Fig. 2 (left) and varied Bacc. 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 Bacc, the IC model underpredicts the data at ≳10 TeV. We note that the magnetic field in the acceleration region Bacc need not be the same as the average value in the post-shock cooling region.

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

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

thumbnail 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 Ecutoff = 400 TeV − 1.2 PeV, although the required efficiency is an obstacle for the hadronic model. The standard setting, Lw = 1039 erg s−1 and next = 100 cm−3, 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, Ecutoff = 1.2 PeV). Harder injection spectra do not account well for the steepening at ≲1 TeV. The efficiency depends inversely on Lw and nint. 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 Rb ≲ 60 pc. From this constraint and the efficiency requirement set above, it follows that next ≳ 300 cm−3 (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, Lsh, can be obtained from the ratio of the residence time in the shell, tres, to the proton cooling time (see Eq. (7)). This yields

(14)

where we have assumed that 1% of Lw 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−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 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 Bacc ≲ 2 μG. A higher Bacc 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 Bacc ∼ 1 − 2 μG fulfils the Hillas limit and our requirement of a strong termination shock. If Bacc ≳ 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.


1

To be specific, 16 WN and eight WC stars. WN and WC are WR types for which nitrogen and carbon dominate the spectrum, respectively.

2

The sample includes 18 WR stars (16 WN and two WC, Crowther et al. 2006), six YHGs (Clark et al. 2005), and the 26 OB super-giants that have an absolute magnitude larger than −6 in Negueruela et al. (2010).

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

  1. Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2021, Nat. Astron., 5, 465 [NASA ADS] [CrossRef] [Google Scholar]
  2. Abramowski, A., Acero, F., Aharonian, F., et al. 2012, A&A, 548, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Abramowski, A., Aharonian, F., Ait Benkhali, F., et al. 2015, Science, 347, 406 [NASA ADS] [CrossRef] [Google Scholar]
  4. Ackermann, M., Ajello, M., Allafort, A., et al. 2011, Science, 334, 1103 [NASA ADS] [CrossRef] [Google Scholar]
  5. Actis, M., Agnetta, G., Aharonian, F., et al. 2011, Exp. Astron., 32, 193 [NASA ADS] [CrossRef] [Google Scholar]
  6. Aharonian, F. A. 2004, Very High Energy Cosmic Gamma Radiation (Singapore: World Scientific Publishing Co., Pte. Ltd.) [CrossRef] [Google Scholar]
  7. Aharonian, F., Yang, R., & de Oña Wilhelmi, E. 2019, Nat. Astron., 3, 561 [Google Scholar]
  8. Aharonian, F., Ashkar, H., Backes, M., et al. 2022, A&A, 666, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Beasor, E. R., Davies, B., Smith, N., Gehrz, R. D., & Figer, D. F. 2021, ApJ, 912, 16 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bhadra, S., Gupta, S., Nath, B. B., & Sharma, P. 2022, MNRAS, 510, 5579 [NASA ADS] [CrossRef] [Google Scholar]
  11. Blumenthal, G. R., & Gould, R. J. 1970, Rev. Mod. Phys., 42, 237 [Google Scholar]
  12. Brandner, W., Clark, J. S., Stolte, A., et al. 2008, A&A, 478, 137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Bykov, A. M. 2001, Space Sci. Rev., 99, 317 [NASA ADS] [CrossRef] [Google Scholar]
  14. Cesarsky, C. J., & Montmerle, T. 1983, Space Sci. Rev., 36, 173 [NASA ADS] [CrossRef] [Google Scholar]
  15. Chu, Y. H. 2008, IAU Symp., 250, 341 [NASA ADS] [Google Scholar]
  16. Clark, J. S., Negueruela, I., Crowther, P. A., & Goodwin, S. P. 2005, A&A, 434, 949 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Clark, J. S., Ritchie, B. W., & Negueruela, I. 2020, A&A, 635, A187 [EDP Sciences] [Google Scholar]
  18. Crowther, P. A., Hadfield, L. J., Clark, J. S., Negueruela, I., & Vacca, W. D. 2006, MNRAS, 372, 1407 [Google Scholar]
  19. Davies, B., & Beasor, E. R. 2019, MNRAS, 486, L10 [NASA ADS] [CrossRef] [Google Scholar]
  20. Drury, L. O. 1983, Rep. Progr. Phys., 46, 973 [Google Scholar]
  21. Ferrand, G., & Marcowith, A. 2010, A&A, 510, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Funk, S. 2015, Annu. Rev. Nucl. Part. Sci., 65, 245 [NASA ADS] [CrossRef] [Google Scholar]
  23. Gabici, S., Evoli, C., Gaggero, D., et al. 2019, Int. J. Mod. Phys. D, 28, 1930022 [CrossRef] [Google Scholar]
  24. Gupta, S., Nath, B. B., & Sharma, P. 2018a, MNRAS, 479, 5220 [NASA ADS] [CrossRef] [Google Scholar]
  25. Gupta, S., Nath, B. B., Sharma, P., & Eichler, D. 2018b, MNRAS, 473, 1537 [NASA ADS] [Google Scholar]
  26. Hillas, A. M. 1984, ARA&A, 22, 425 [Google Scholar]
  27. Hinton, J., & SWGO Collaboration 2022, 37th International Cosmic Ray Conference, 23 [Google Scholar]
  28. Kafexhiu, E., Aharonian, F., Taylor, A. M., & Vila, G. S. 2014, Phys. Rev. D, 90, 123014 [Google Scholar]
  29. Klepach, E. G., Ptuskin, V. S., & Zirakashvili, V. N. 2000, Astropart. Phys., 13, 161 [NASA ADS] [CrossRef] [Google Scholar]
  30. Kothes, R., & Dougherty, S. M. 2007, A&A, 468, 993 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Lim, B., Chun, M.-Y., Sung, H., et al. 2013, AJ, 145, 46 [NASA ADS] [CrossRef] [Google Scholar]
  32. Limongi, M., & Chieffi, A. 2006, ApJ, 647, 483 [NASA ADS] [CrossRef] [Google Scholar]
  33. Mac Low, M.-M., & McCray, R. 1988, ApJ, 324, 776 [NASA ADS] [CrossRef] [Google Scholar]
  34. Mestre, E., de Oña Wilhelmi, E., Torres, D. F., et al. 2021, MNRAS, 505, 2731 [CrossRef] [Google Scholar]
  35. Morlino, G., Blasi, P., Peretti, E., & Cristofari, P. 2021, MNRAS, 504, 6096 [NASA ADS] [CrossRef] [Google Scholar]
  36. Muno, M. P., Law, C., Clark, J. S., et al. 2006, ApJ, 650, 203 [NASA ADS] [CrossRef] [Google Scholar]
  37. Negueruela, I., Clark, J. S., & Ritchie, B. W. 2010, A&A, 516, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Negueruela, I., Alfaro, E. J., Dorda, R., et al. 2022, A&A, 664, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Oey, M. S. 2009, AIP Conf. Ser., 1156, 295 [NASA ADS] [CrossRef] [Google Scholar]
  40. Ohm, S., Hinton, J. A., & White, R. 2013, MNRAS, 434, 2289 [NASA ADS] [CrossRef] [Google Scholar]
  41. Popescu, C. C., Yang, R., Tuffs, R. J., et al. 2017, MNRAS, 470, 2539 [NASA ADS] [CrossRef] [Google Scholar]
  42. Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431 [NASA ADS] [CrossRef] [Google Scholar]
  43. Seo, J., Kang, H., & Ryu, D. 2018, J. Korean Astron. Soc., 51, 37 [NASA ADS] [Google Scholar]
  44. Stevens, I. R., & Hartwell, J. M. 2003, MNRAS, 339, 280 [NASA ADS] [CrossRef] [Google Scholar]
  45. Vieu, T. 2021, Ph.D. Thesis, Université de Paris, France [Google Scholar]
  46. Vieu, T., Gabici, S., Tatischeff, V., & Ravikularaman, S. 2022a, MNRAS, 512, 1275 [NASA ADS] [CrossRef] [Google Scholar]
  47. Vieu, T., Reville, B., & Aharonian, F. 2022b, MNRAS, 515, 2256 [NASA ADS] [CrossRef] [Google Scholar]
  48. Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377 [Google Scholar]
  49. Yang, R.-Z., de Oña Wilhelmi, E., & Aharonian, F. 2018, A&A, 611, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Tables

Table 1.

Input parameters used in the default scenario.

All Figures

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

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.