© The Authors 2025

1 Introduction

Star clusters (SCs) have long been proposed as possible sources contributing to the Galactic cosmic ray (CR) flux. The interest in such sources is mainly due to the fact that SCs might host the right conditions for the acceleration of CRs up to the knee and provide an explanation of the well-known ²²Ne/²⁰Ne excess in CRs (Gupta et al. 2020).

Particle acceleration could occur at the termination shock of the collective wind launched by the young stars located in the cluster core, if the SC is sufficiently compact (Morlino et al. 2021; Blasi & Morlino 2023), or occur due to supernova explosions (Vieu et al. 2022), for sufficiently old SCs. Second order Fermi acceleration in the bubble has also been invoked as a possible mechanism of particle acceleration (Vink 2024). The maximum energy that can be achieved at the termination shock depends rather crucially on details of CR transport in the bubble. For instance, as discussed by Morlino et al. (2021); Blasi & Morlino (2023), although a formal estimate of the maximum energy, E_max, attainable at the termination shock can easily approach ~Peta-electronVolt, the spherical symmetry of the problem introduces a curvature in the spectrum of accelerated particles that starts at energies much lower than E_max. The only exception to this conclusion is the case of Bohm diffusion, but it is unlikely to apply to the case of extrinsic turbulence.

From the observational point of view, gamma ray telescopes such as the High-Altitude Water Cherenkov Observatory (HAWC) (Abeysekara et al. 2021) and the Large High Altitude Air Shower Observatory (LHAASO) (Cao et al. 2021) have measured the flux of gamma rays from selected SCs, and the spectrum appears to drop appreciably at energies much lower than Peta-ElectronVolt despite the fact that photons with energy as high as 2.5 PeV have recently been detected by LHAASO from a region of ~6° in the direction of the Cygnus OB association (Lhaaso Collaboration 2024). Gamma ray emission has now been detected from a number of SCs, such as Westerlund 1 (Abramowski et al. 2012), Westerlund 2 (Yang et al. 2018), the Cygnus cocoon (Ackermann et al. 2011; Aharonian et al. 2019), NGC 3603 (Saha et al. 2020), BDS2003 (Albert et al. 2021), W40 (Sun et al. 2020), and 30 Doradus in the Large Magellanic Cloud (H. E. S. S. Collaboration 2015). Several very young SCs have also been detected in gamma rays in the ≲100 GeV energy range (Peron et al. 2024a).

Detailed calculations of particle acceleration and transport in the bubble excavated by the wind carried out by Blasi & Morlino (2023) (see also Menchiari (2023)) and the more phenomenological estimates presented in Cao et al. (2021) suggest that in order to account for the level of gamma ray emission in the Cygnus SC, a mean density of ~5–10 cm^–3 is needed in the bubble. Such gas is not the diffuse gas associated with the collective wind of the SC, which has a much lower density, typically ~10^–3 cm^–3. The mean density should rather be dominated by atomic and molecular gas clumped in dense regions with a larger local density. Such clumps can either be the leftover of the passage of the shock associated with the collective wind or formed as a result of the cooling of the gas plowed away toward the outer edge of the bubble. Accepting the presence of this gas and attributing the gamma ray emission to hadronic interactions with such gas, Blasi & Morlino (2024) showed that the spectra of H and He nuclei escaping the SC are different, which is in agreement with what was inferred from CR data collected by the CREAM collaboration (Ahn et al. 2010), by PAMELA (Adriani et al. 2011), AMS-02 (Aguilar et al. 2015b,a), DAMPE (An et al. 2019; Ale- manno et al. 2021), and CALET (Adriani et al. 2022, 2023). The disappointing by-product of this finding is that, in this case, heavier nuclei are severely depleted due to spallation inside the bubble, so SCs cannot contribute an appreciable flux of nuclei heavier than He.

This was the first symptom of a more general implication of the detection of gamma rays from Cygnus and other clusters: if of hadronic origin, this radiation forces us to conclude that CRs in SCs need to cross a large grammage inside the bubble. For the reference values of the mean density quoted by Cao et al. (2021), such grammage exceeds the total that CRs are expected to accumulate in the Galaxy while diffusing out. One could argue that this conclusion may differ by changing the model of CR injection. Hence, below we consider three models for the origin of the gamma ray emission in an SC. In the first model, CRs are all accelerated in the center of the SC and diffuse outward while interacting with dense clouds (this is the reference model considered by Cao et al. (2021) for the Cygnus region). In the second model, particle acceleration occurs at the termination shock, and particles advect and diffuse outward (Morlino et al. 2021). In the third model, CRs are assumed to be injected in an impulsive event in the center and diffuse outward. The event might be a supernova explosion in the core of the cluster.

We show that in all models listed above, the local grammage traversed by CRs with E ≳ 10 TeV while escaping the SC exceeds the one inferred from secondary/primary ratios (Evoli et al. 2019, 2020; Schroer et al. 2021). The implications for lower energy CRs are also discussed. In light of this finding, we are forced to conclude that either SCs only contribute a negligible fraction of the CRs observed at the Earth or CRs originate in SCs and also accumulate their grammage there rather than in the interstellar medium (ISM) on Galactic scales. This second scenario would require a major revision of the foundations of the model of CR transport in the Galaxy. Alternatively, the gamma ray emission that we measure might be due to leptonic processes, which makes the predictions independent of the gas density. These models (see, e.g., Härer et al. (2023) for such a model applied to Westerlund 1) require a lower strength of the magnetic field in the bubble, which in turn implies that SCs cannot reasonably play the role of PeVatrons. Finally, it might be that the SCs that have been detected so far in gamma rays are unusual, and thus their characteristics are irrelevant for the bulk of the SCs. Although this cannot be excluded, it would definitely be surprising.

If the constraints to the grammage that are inferred here are meaningful for the bulk of SCs, then one should keep them in mind even when supernova remnants (SNRs) are considered as the main source of CRs, since the majority of such sources are expected to be located in SCs. From the observational point of view, an important clue to the origin of the gamma radiation and the constraints on grammage might come from the detection of gamma rays from very young SCs, for which we can be reasonably sure of the absence of SN events taking place at the present time.

This article is organized as follows: in Sect. 2 we discuss the gamma ray emission in the three models of CR production illustrated above. In Sect. 3, we elaborate on the implications of the gamma ray measurements for CR grammage, comparing the results with standard models of CR transport in the Galaxy. In Sect. 4, we outline our conclusions.

2 Models of gamma ray production in star clusters

As we discuss below, the effect of the grammage is so macroscopic that a sophisticated treatment of proton-proton interactions is not necessary. Moreover, in the analysis below, we focus on gamma ray energies >1 TeV, where scaling relations are sufficiently accurate for our purposes here. We assumed that gamma rays with a given energy Eγ are produced by protons with energy E such that E_γ = ηE, with η ≃ 0.1. The cross section for inelastic pp interactions is considered to be constant with energy, σ_pp = 32 mb.

Here we consider three models of CR transport in the SC. In Model 1 we assumed, following (Lhaaso Collaboration 2024), that a source in the center of the SC continuously injects CRs with a spectrum Q(E) = A(E/m_p)^–α, with α > 2. We also follow (Lhaaso Collaboration 2024) in assuming that transport is purely diffusive with the diffusion coefficient D(E). In the energy region of interest here, E_γ ≳ 1 TeV, the diffusion equation can be used in its stationary form so that the CR density of energy E at a distance r from the center is $$n_{CR}(r,E)\approx \frac{Q(E)}{4\pi rD(E)}.$$(1)

We normalized the injection spectrum to a fraction ξ_CR of the luminosity L_SC of the SC (the one that involves plasma motion and can be the seed of particle acceleration) through $$Q(E)={\xi }_{CR}\frac{L_{SC}(\alpha -2)}{{\left(m_p{c}^2\right)}^2}{\left(\frac{E}{m_p}\right)}^{-\alpha },$$(2)

where we assumed that the energy in the form of CRs is dominated by relativistic particles, namely E ≳ m_p. Within a factor of order unity, this is a good approximation for 2 ≲ α ≲ 3, as expected from diffusive shock acceleration (DSA). These source spectra are also compatible with those found in Lhaaso Collaboration (2024) that provide a good fit to the observed gamma ray emission for the Cygnus region. We note that, quantitatively, including particles with energy lower than m_p in the normalization procedure would strengthen the conclusions drawn below.

In the second model (Model 2), we follow Morlino et al. (2021); Blasi & Morlino (2023) in assuming that particle acceleration takes place at the termination shock of the collective wind. The structure of the bubble is properly described by Morlino et al. (2021): the location of the termination shock is at r = Rsh , while the bubble extends to r = Rb . However, as we demonstrate below, these quantities do not need to be written down explicitly. Particle acceleration is assumed to transform a fraction ξ_CR of the luminosity of the wind, $L_w=(1/2)\dot{M}{v}_w^2$, into CRs. Here, Ṁ is the rate of mass loss in the form of a collective wind, and v_w is the wind velocity. The structure of the bubble is such that the shocked wind (i.e., the plasma behind the termination shock) has a velocity v(r) = (v_w/4)(r/R_sh)^-2, assuming that the termination shock is strong (compression factor of four). The density downstream of the shock is constant. In the wind region (upstream of the termination shock), the velocity is constant to the value v_w , while the density is $${\rho }_w(r)=\frac{\dot{M}}{4\pi r^2{v}_w}.$$(3)

Assuming again that the density of accelerated particles at the termination shock has the shape of a power law, n_CR(E) = A(E/m_p)^–α, and that the energy in the form of accelerated particles carries a fraction ξ_CR of L_w, it is easy to derive $$n_{CR}(E)=2{\xi }_{CR}\frac{L_w(\alpha -2)}{4\pi R_{sh}^2{\left(m_p{c}^2\right)}^2}{\left(\frac{E}{m_p}\right)}^{-\alpha }.$$(4)

In the following, we assume that this density remains constant downstream of the termination shock. While this is a better approximation at E ≲ 1 TeV than at higher energies (where escape from the bubble becomes dominated by diffusion rather than advection), it also maximizes the gamma ray emission, making the results illustrated below even stronger.

In the third model (Model 3), the particles are assumed to be produced in the core of the SC, as in Model 1, but the source is impulsive. This situation might mimic the case of an SN explosion. It is not clear how likely it is to have had an SN explosion in Cygnus or Westerlund 1, although there is evidence for a pulsar in the core in the Cygnus region, with an estimated age of ~10⁵ yr. Since there is no evidence of a shell in the bubble, we assumed that the SN event occurred ≲10⁵ years ago. Once produced, CRs are assumed to diffuse outward and eventually escape.

If E_SN is the kinetic energy of the SN ejecta, then the density of particles in the bubble can be estimated as $$n_{CR}(E)\approx \frac{{\xi }_{CR}{E}_{SN}(\alpha -2)}{{\left(m_p{c}^2\right)}^2{[4\pi D(E)t]}^{3/2}}{\left(\frac{E}{m_p}\right)}^{-\alpha },\text{r}<\sqrt{4Dt}.$$(5)

We note that the condition that particles of energy E are still inside the SC at time t imposes an upper limit on the diffusion coefficient D(E): $$D(E)<7.7\times {10}^{27} {\text{cm}}^2/\text{s}{\left(\frac{\text{t}}{{10}^5\text{yr}}\right)}^{-1}{\left(\frac{{\text{R}}_{\text{b}}}{100\text{\hspace{0.17em}pc}}\right)}^2.$$(6)

We also note that the particles at time t uniformly fill a sphere of radius $\sqrt{4Dt}$, while their density drops∝ exp[–r²/4Dt] at larger distances. For simplicity, we assumed that the particles are all concentrated with the diffusion distance.

2.1 Gamma ray emission in Model 1: Diffusion model

In the following, we denote as n_gas the mean gas density in the bubble, which is to be interpreted as the volume averaged gas density, dominated by the clumpy structure of atomic and molecular gas (see discussion in Blasi & Morlino 2024). We assumed that this density is constant throughout the bubble since the available information on the spatial distribution of the gas is very poor (also see the discussion in Menchiari 2023).

As discussed above, in Model 1 (Lhaaso Collaboration 2024), we assume that a source in the SC core injects CRs at a constant rate and that these particles move outward only due to diffusion. In this simple picture there is no bulk motion of the plasma, so it mimics a situation where the SC is unable to launch a collective wind (although even in that case it is difficult to imagine that a collective motion may be completely absent). As argued by Vieu et al. (2024), this might be the case of the Cygnus cocoon.

The emissivity in the form of gamma rays can be written as $$J_{\gamma }\left(E_{\gamma },r\right)d{E}_{\gamma }=n_{CR}(E,r)dE{n}_{𝑔as}c{\sigma }_{pp}.$$(7)

Using the assumption that E_γ = ηE, one can write $$J_{\gamma }\left(E_{\gamma },r\right)E_{\gamma }^2={\xi }_{CR}\frac{L_{SC}(\alpha -2)}{4\pi rD(E)}{\eta }^{\alpha -1}{\left(\frac{E_{\gamma }}{m_p}\right)}^{-\alpha +2}{n}_{𝑔as}c{\sigma }_{pp}.$$(8)

It follows that the gamma ray flux from the entire bubble can be estimated as $$\begin{array}{l}{E}_{\gamma }^2{\text{Φ}}_{\gamma }\left(E_{\gamma }\right)={\displaystyle {\int }_0^{R_b}d}r\frac{r^2}{d^2}{\xi }_{CR}\frac{L_{SC}(\alpha -2)}{4\pi rD(E)}{\eta }^{\alpha -1}{\left(\frac{E_{\gamma }}{m_p}\right)}^{-\alpha +2}\times n_{𝑔as}C{\sigma }_{pp}\\ ={\xi }_{CR}\frac{L_{SC}(\alpha -2)}{4\pi d^{2}D(E)}{\eta }^{\alpha -1}{\left(\frac{E_{\gamma }}{m_p}\right)}^{-\alpha +2}\frac{R_b^2}{2}{n}_{𝑔as}c{\sigma }_{pp}.\end{array}$$(9)

The escape time from the bubble in a purely diffusive case is ${\tau }_b=R_b^2/6D(E)$, so we can introduce the grammage traversed by particles while escaping the bubble as $$X(E)=n_{𝑔as}{m}_{p}c{\tau }_b(E)=n_{𝑔as}{m}_{p}c\frac{R_b^2}{6D(E)}.$$(10)

This allows us to rewrite Eq. (9) as $$E_{\gamma }^2{\text{Φ}}_{\gamma }\left(E_{\gamma }\right)=\frac{3{\xi }_{CR}{L}_{SC}(\alpha -2){\eta }^{\alpha -1}}{4\pi d^2}{\left(\frac{E_{\gamma }}{m_p}\right)}^{-\alpha +2}\frac{X(E)}{X_{cr}},$$(11)

where X_cr = m_p/σ_pp is the critical grammage for pp interactions. As written here, this expression can also be used, within factors of order unity, when the gas is not uniformly distributed, for instance when the whole gas is concentrated in a thin shell at the edge of the bubble.

2.2 Gamma ray emission in Model 2: Acceleration at the termination shock

With the same assumptions on pp collisions adopted above, the gamma ray emissivity in the case of particle acceleration at the termination shock can be written as $$J_{\gamma }\left(E_{\gamma }\right)=2{\xi }_{CR}\frac{L_w(\alpha -2)}{4\pi R_{sh}^2{\left(m_p{c}^2\right)}^2}{\eta }^{\alpha -1}{\left(\frac{E_{\gamma }}{m_p}\right)}^{-\alpha }{n}_{𝑔as}c{\sigma }_{pp}.$$(12)

The integral over the downstream of the termination shock, the region that dominates the gamma ray emission in this model, leads to $$E_{\gamma }^2{\text{Φ}}_{\gamma }\left(E_{\gamma }\right)=2\frac{{\xi }_{CR}{L}_w(\alpha -2)}{4\pi R_{sh}^2}{\eta }^{\alpha -1}{\left(\frac{E_{\gamma }}{m_p}\right)}^{-\alpha +2}\frac{R_b^3}{3d^2}{n}_{𝑔as}c{\sigma }_{pp}.$$(13)

To maximize the gamma ray flux predicted by the model and make our conclusions more solid, we assumed that particle escape from the bubble is dominated by advection. As discussed by Blasi & Morlino (2024), escape becomes dominated by diffusion at energies ≳1 TeV, but this leads to a smaller effective emission region for gamma rays and correspondingly smaller gamma rays ray flux; hence, our assumption is conservative in terms of deriving constraints on the grammage inside the bubble. Recalling that the plasma velocity drops downstream of the termination shock, we can write $$\frac{dr}{dt}=\frac{v_w}{4}{\left(\frac{r}{R_{sh}}\right)}^{-2}\to {\tau }_{adv}\approx \frac{4}{3}\frac{R_{sh}}{v_w}{\left(\frac{R_b}{R_{sh}}\right)}^3,$$(14)

where for simplicity we assumed that R_b >> R_sh. Replacing this expression in Eq. (13), we easily obtain $$E_{\gamma }^2{\text{Φ}}_{\gamma }\left(E_{\gamma }\right)\approx \frac{2{\xi }_{CR}{L}_w(\alpha -2){\eta }^{\alpha -1}}{4\pi d^2}{\left(\frac{E_{\gamma }}{m_p}\right)}^{-\alpha +2}\frac{X(E)}{X_{cr}},$$(15)

where the grammage is defined as X(E) = n_gasm_pcτ_adv. A quick comparison between Eqs. (11) and (15) immediately suggests that Models 1 and 2 are going to lead to very similar constraints in terms of grammage.

Fig. 1 Lower limits to the grammage inferred for the Cygnus region (left) and Westerlund 1 (right) if the observed gamma rays are of hadronic origin in three models of gamma ray production as described in the text. These limits are shown together with an up-to-date computation of the Galactic CR grammage inferred from secondary/primary ratios (Schroer et al. 2021). The yellow shadowed region shows the maximum source grammage still compatible with secondary to primary ratios (≃0.4 g cm–2).

2.3 Gamma ray emission in Model 3: Impulsive injection of cosmic rays

The detection of gamma rays of energy E_γ implies that the parent CR particles, with energy E = E_γ/η, are still confined inside the SC. This trivially implies that the gamma ray spectrum integrated over the volume is insensitive to the diffusion coefficient. Hence the flux of gamma rays can be written using Eq. (5) integrating on $r\le \sqrt{4D(E)t}$: $$E_{\gamma }^2{\text{Φ}}_{\gamma }\left(E_{\gamma }\right)=\frac{{\xi }_{CR}{E}_{SN}(\alpha -2){\eta }^{\alpha -1}}{3{\pi }^{3/2}{d}^2}{\left(\frac{E_{\gamma }}{m_p}\right)}^{-\alpha +2}{n}_{𝑔as}{\sigma }_{pp}c,$$(16)

with the usual meaning of the symbols. It is worth noting how this result, contrary to the ones obtained for the other two models, depends only on the density of the target gas and not on the grammage. On the other hand, as discussed above (see Eq. (6)), the condition that particles are still trapped in the bubble after a time t allows us to infer a lower limit to the grammage traversed by particles before they are released into the ISM.

3 Constraints on the grammage inside the Cygnus cocoon and Westerlund 1

The LHAASO experiment has recently measured the gamma ray flux with energy above ~1 TeV from the Cygnus cocoon at the level of $E_{\gamma }^2{\text{Φ}}_{\gamma }\left(E_{\gamma }\right)\simeq 6\times {10}^{-11}\text{erg\hspace{0.17em}\hspace{0.17em}}{\text{cm}}^{-2} {\text{s}}^{-1}$ over a region that extends ~150 pc from the center of Cygnus SC. The emission is claimed to have a hadronic origin, and it is difficult to imagine otherwise given the extended morphology. The slope of the spectrum is ~2.7 (at E ≳ 1 TeV) and probably reflects the energy-dependent diffusion coefficient, together with the shape of the spectrum near the maximum energy. This point, however, is of little interest here in that we focus on the lowest energy bin, at 1 TeV, and investigate the dependence of our constraints on the value of α, the slope of the injection spectrum. The HESS telescope has observed gamma rays from the direction of the SC Westerlund 1 (Aharonian et al. 2022), with $E_{\gamma }^2{\text{Φ}}_{\gamma }\left(E_{\gamma }\right)\simeq 1.6\times {10}^{-11}\text{erg\hspace{0.17em}}{\text{cm}}^{-2} {\text{s}}^{-1}$ at E_γ = 1 TeV. In the following, we adopt a distance of 1.4 kpc for Cygnus and 3.9 kpc for Westerlund 1.

In Model 1, we normalized the luminosity to L_SC = 10³⁹L³⁹ erg/s. We note that photons with energy Eγ are produced by protons with energy E = 10 TeV; hence, we are able to constrain the grammage traversed in the SC by protons with 10 TeV energy. The gamma ray flux computed using Eq. (11) can be written as $E_{\gamma }^2{\text{Φ}}_{\gamma }\approx 1.5\times {10}^{-9}{\xi }_{CR}{L}_{39}{d}_{kpc}^{-2}X(E)\text{erg\hspace{0.17em}}{\text{cm}}^{-2} {\text{s}}^{-1}(4.8\times {10}^{-10}{\xi }_{CR}{L}_{39}{d}_{kpc}^{-2}X(E)\text{erg\hspace{0.17em}}{\text{cm}}^{-2} {\text{s}}^{-1})$ for α = 2.2 (α = 2.4). Comparing this prediction with the flux as measured by LHAASO for the Cygnus cocoon, we obtain a grammage $X(E=10TeV)\approx 0.08{\xi }_{CR}^{-1}{L}_{39}^{-1} \text{g} {\text{cm}}^{-2}\left(0.25{\xi }_{CR}^{-1}{L}_{39}^{-1} \text{g} {\text{cm}}^{-2}\right)\text{for}\alpha =2.2(\alpha =2.4)$

Finally, when requiring a reasonable efficiency of particle acceleration, say ξ_CR ≲ 0.1, we immediately inferred that in order to account for the gamma ray flux observed by LHAASO, the grammage required for Model 1 has to be $X(E=10TeV)\gtrsim 0.8L_{39}^{-1} \text{g} {\text{cm}}^{-2}$ for α = 2.2 and $X(E=10TeV)\gtrsim 2.5L_{39}^{-1} \text{g} {\text{cm}}^{-2}$ for α = 2.4. A similar procedure leads to a lower limit to the grammage in the case of Westerlund 1. Both sets of limits are plotted in Fig. 1 for Cygnus (left panel) and Westerlund 1 (right panel). The figure also shows the grammage traversed by CRs during their journey through the Galaxy, as inferred from sec- ondary/primary ratios (see for instance Evoli et al. 2019, 2020; Schroer et al. 2021). The yellow shaded region shows the largest grammage (the so-called source grammage) that would still be compatible with the measured ratios. Interestingly, this grammage is not too far from the one expected to be accumulated while CRs are downstream of the shock of a SNR (Aloisio et al. 2015). Any grammage in excess of ~0.4g cm^–2 would force us to revisit the very foundations of the picture of CR transport in the Galaxy.

In fact, this seems to be the case for both Cygnus and Wester- lund 1 if the gamma ray emission is of hadronic origin. Model 1 leads us to expect a grammage in both cases appreciably larger than the Galactic grammage at 10 TeV. Even more concerning, in Model 1, this grammage would extend to lower energies as a power law, so the conclusions above would also apply to lower energies, where the grammage is even better constrained.

Finally, we note that the limit above has been specialized to an SC luminosity of 10³⁹ erg/s, which is probably somewhat larger than the actual luminosity of the two SCs considered here. For lower values of the luminosity, the limits would become correspondingly stronger.

In Model 2, we normalized the wind luminosity to L_w = 10³⁹L_w,39 erg/s. The gamma ray flux computed using Eq. (15) can be written as $E_{\gamma }^2{\text{Φ}}_{\gamma }\approx {10}^{-9}{\xi }_{CR}{d}_{kpc}^{-2}{L}_{w,39}X(E)\text{erg\hspace{0.17em}}{\text{cm}}^{-2} {\text{s}}^{-1}(3.2\times {10}^{-10}{\xi }_{CR}{d}_{kpc}^{-2}{L}_{w,39}X(E){\text{erg\hspace{0.17em}cm}}^{-2} {\text{s}}^{-1})$ for α = 2.2 (α = 2.4). Following the same procedure outlined for Model 1 and requiring efficiency ξ_CR < 0.1, we inferred the corresponding lower limits on the grammage, as shown in Fig. 1 (green symbols for Model 2). The problem discussed above with the grammage appears to be even more severe for Model 2.

We note that the limits derived for Model 2 should be considered as absolute lower limits to the grammage since we neglected the diffusive escape of CRs from the cocoon at E ≳ 1 TeV (Blasi & Morlino 2024), which would reduce the high energy CR density and increase the required grammage. These limits are in line with the results presented by Blasi & Morlino (2023, 2024), where a gas density n_gas ~ 5–10 cm^–3 was used, which is in line with the values inferred by the Lhaaso Collaboration (2024). Using the expressions above, this would correspond to a CR acceleration efficiency at the termination shock ξ_CR ~ 0.5%, in agreement with the values found by Blasi & Morlino (2023). This, however, implies a grammage in the cocoon of ~50 g cm^-2 and is thus well in excess of the Galactic grammage. Not surprisingly, such a grammage was found to be responsible for the complete spallation of heavy nuclei and a slight difference between the escaping spectra of H and He nuclei (Blasi & Morlino 2024). SCs with characteristics similar to those of the Cygnus cocoon and Westerlund 1 cannot account for the bulk of Galactic CRs unless CR transport is deeply revised. Alternatively, the gamma ray emission may not be due to hadronic interactions, as suggested by Härer et al. (2023) for Westerlund 1.

In Model 3, we tested the possibility that CRs have been injected in the bubble during an impulsive event at a time t in the past and that the gamma ray emission is being produced by such particles today. The condition that CRs of energy E are still confined inside the bubble leads to an upper limit to the diffusion coefficient, $D(E)<R_b^2/4t$. Comparing between Eq. (16) and the observed gamma ray flux from an SC allowed us to impose a lower limit to the gas density in the bubble by requiring ξ_CR ≤ 0.1. For Cygnus, this limit reads n_gas > N_min = 62 cm^–3 for α = 2.2 and n_gas > n_min = 131 cm^–3 for α = 2.4. The combination of these two conditions implies that the grammage must satisfy $$X(E)>\frac{2}{3}{n}_{min}{m}_{p}ct=1 \text{g} {\text{cm}}^{-2}\left(\frac{{\text{n}}_{\min }}{10 {\text{cm}}^{-3}}\right)\left(\frac{\text{t}}{{10}^5\text{yrs}}\right).$$(17)

One can see that not only the values of density in this model appear extreme and hard to reconcile with the formation of a bubble, but also the inferred grammage exceeds the Galactic grammage, more so than in Models 1 and 2 discussed above. In Fig. 1, we only show the limit on the grammage for Model 3 for the cases in which the limit is lower than 10 g cm^-2. Similar considerations apply to the case of Westerlund 1.

4 Discussion

Star clusters have long been thought to play an important role in the framework of the origin of Galactic CRs in at least two ways: first, by providing enough CRs to justify the excess of ²²Ne observed in the cosmic radiation and, second, by possibly providing the right conditions for acceleration up to the knee, something that SNRs seem to struggle to do. The recent detection of selected SCs in high energy gamma rays has fueled renewed interest in this problem, especially in the perspective of SCs playing the role of PeVatrons. As discussed by Morlino et al. (2021); Blasi & Morlino (2023), it seems unlikely that SCs would efficiently accelerate particles up to the knee energy at the termination shock generated by the collective wind. It cannot yet be excluded that higher energies may be reached in SN explosions occurring inside an SC. This latter scenario would apply to SCs that are old enough to host SN explosions.

Here, we have discussed the implications of the recently observed gamma ray emission from selected SCs in terms of their contribution to the bulk of Galactic CRs. In particular, we have shown that the gamma ray emission observed by LHAASO from the direction of the Cygnus cocoon (Cao et al. 2021; Lhaaso Collaboration 2024) and by HESS from the direction of West- erlund 1 (Aharonian et al. 2022) lead to an estimate of the grammage accumulated by particles while leaving the bubble region that exceeds the Galactic grammage at ~10 TeV. The way this finding extends to lower energies depends on the model adopted to describe the gamma ray emission. Even if this grammage were energy independent, the accumulated grammage at low energies would require substantial revision of the transport of Galactic CRs in order to accommodate the excess grammage.

If one takes for granted the average gas density in the Cygnus cocoon, n_gas = 5–10 cm^–3, as inferred in much of the recent literature (Aharonian et al. 2019; Cao et al. 2021; Menchiari 2023; Menchiari et al. 2024a), the inferred grammage would be X ~ 25–50 g cm^–2, which is larger than the total grammage traversed by CRs at ~ 10 GeV by about one order of magnitude. Not surprisingly, the calculations of Blasi & Morlino (2024), applied to the Cygnus cocoon, imply that nuclei heavier than He are basically destroyed inside the cocoon, and even for He the spallation reactions are sufficiently severe to harden the He spectrum with respect to H nuclei.

It is clear that the discovery of high energy gamma ray emission from SCs and most notably from Cygnus and Westerlund 1 represents a milestone in our field of investigation, and it provides a possible clue to the origin of high energy CRs. But the same detection also implies that either SCs cannot contribute but a small fraction of Galactic CRs or our understanding of how CRs accumulate grammage in the Galaxy needs to be deeply revised. Clearly these conclusions are based on assumptions that may be flawed. For instance, one could speculate that the SCs detected in gamma rays are not representative of the whole class of SCs, implying that perhaps most SCs are not bright in gamma rays. The other possibility is that the gamma radiation observed so far is of leptonic origin, so it becomes uncorrelated with the grammage traversed in the source.

We note that if most SCs are hadronic gamma ray emitters, as discussed by Menchiari et al. (2024b), they would contribute significantly to the diffuse Galactic gamma ray background. On the other hand, while there is substantial correlation between unidentified Fermi-LAT sources and young SCs, no such strong correlation arises in the TeV range or with SCs that are older than a few million years (Peron et al. 2024b), suggesting that Cygnus and Westerlund might be exceptional. Clearly both of these issues require more attention to assess the role of SCs as potential CR sources.

Models of the origin of Galactic CRs in which grammage is accumulated near sources were actually developed by Cowsik & Wilson (1975) and recently revised by Cowsik & Huth (2022) as a possible way to address the issue of a rising positron fraction. In general, these models struggle to accommodate the features observed in the spectra of primary nuclei and the decay of ¹⁰Be. Nevertheless, it is worth keeping in mind the possibility that at least part of the grammage inferred from measurements at the Earth may be accumulated in or around the sources rather than en route to the Earth.

Acknowledgements

This work has been partially funded by the European Union – Next Generation EU, through PRIN-MUR 2022TJW4EJ and by the European Union – NextGenerationEU under the MUR National Innovation Ecosystem grant ECS00000041 – VITALITY/ASTRA – CUP D13C21000430001. The author is grateful to E. Amato, C. Evoli, S. Menchiari, G. Morlino, G. Peron and E. Sobacchi for useful conversations and to B. Schroer for providing the grammage from Schroer et al. (2021).

References

All Figures

Fig. 1 Lower limits to the grammage inferred for the Cygnus region (left) and Westerlund 1 (right) if the observed gamma rays are of hadronic origin in three models of gamma ray production as described in the text. These limits are shown together with an up-to-date computation of the Galactic CR grammage inferred from secondary/primary ratios (Schroer et al. 2021). The yellow shadowed region shows the maximum source grammage still compatible with secondary to primary ratios (≃0.4 g cm–2).

In the text