© The Authors 2025

1 Introduction

Stellar clusters represent the fundamental building blocks of galaxies and are among the most intensively studied celestial objects in the Cosmos. Among all clusters, the youngest and most massive ones harbor hundreds of OB-type stars, whose violent winds create the perfect environment for the production of cosmic rays (CRs). Indeed, since the 1980s, it has been speculated that particle acceleration may occur in these environments, and several different scenarios have been proposed, such as acceleration by stellar wind-wind interaction (Cesarsky & Montmerle 1983; Klepach et al. 2000; Reimer et al. 2006), acceleration by supersonic turbulence (Bykov et al. 2020), and acceleration at the collective cluster wind termination shock (TS) (Morlino et al. 2021). Moreover, the large majority of massive stars (~96%) are believed to be born in clusters (Parker & Goodwin 2007), even if a non-negligible fraction (20–30%) is expelled from the parent cluster due to dynamical effects (Gies 1987; Parker & Goodwin 2007). Hence, most supernova explosions are expected to happen within a cluster. In such a scenario, particle acceleration can also occur thanks to a combination of winds and supernova explosions (Vieu et al. 2022).

Eventually, the detection, in the last decade, of diffuse γ-ray emission in coincidence with a dozen young massive star clusters (YMSCs), prominent among Cygnus OB2 (Bartoli et al. 2014; Abeysekara et al. 2021; Astiasarain et al. 2023; Lhaaso Collaboration 2024), Westerlund 1 (Abramowski et al. 2012; Aharonian et al. 2022), and Westerlund 2 (H. E. S. S. Collaboration 2011; Yang et al. 2018), has further strengthened the hypothesis that these objects are efficient CR factories. In fact, the observed emissions can be easily explained if ~1–10% of the power supplied by the stellar winds is used to accelerate CRs (Aharonian et al. 2019; Peron et al. 2024).

In general, the size of the detected γ-ray emission is of the order of 1° – 3°, but in some cases it can be extremely large, as in the case of the 100 deg² emission observed by the Large High Altitude Air Shower Observatory (LHAASO) around Cygnus OB2 (Lhaaso Collaboration 2024). Noticeably, in the well-established cases, the size of the emitting region is close to that expected for the wind-blown bubbles thought to surround such objects. The similarity can be easily explained if the accelerated particles remain mostly confined within the turbulent bubble, so that the γ-ray emission traces the projected radius of the latter.

The detection of large extended sources in γ-ray astronomy is often a challenging task due to the limited resolution of the telescopes, which prevents a clear disentanglement of multiple sources in crowded regions of the Galactic Plane. In addition, background subtraction is often problematic due to the relatively small field of view of most current instruments. These observational limitations make the detection of individual YMSCs very difficult.

At present, itis reasonable to expect that most of the emission coming from YMSCs is detected as a non-resolved contribution to the large-scale Galactic emission. In addition, recent analysis of the Galactic diffuse γ-ray emission has highlighted the presence of an excess above a few gigaelectronvolt up to petaelectronvolt energies which is possibly associated to emission from an unresolved population of sources (Zhang et al. 2023; Cao et al. 2023).

In this work, we argue that the contribution from YMSCs represents a non-negligible fraction of the diffuse Galactic emission, and could easily explain the observed γ-ray excess at very-high energies. Since the stellar cluster population is only known within a few kiloparsecs from the Sun (Cantat-Gaudin et al. 2020), in order to estimate such a contribution, we generated multiple synthetic populations of YMSCs based on known properties of observed local clusters. For each synthetic cluster, we simulated a mock population of stars, from which we estimated the parameters of the collective cluster wind. We computed the CR content in each YMSC using the model proposed by Morlino et al. (2021), and we then calculated the associated γ- ray emission due to hadronic interactions only. We include in our model three possible scenarios of plasma turbulence in the windblown bubble, which led to different predictions for the γ-ray spectrum.

While we are well aware that supernova explosions are likely to provide a sizeable contribution to the cluster kinetic energy, comparable to, or even larger than, that of stellar winds, we did not to include CR acceleration by these sources in the current calculation. This choice deserves some comment. In the first place, aimed at estimating a lower limit to the contribution of YMSCs to the diffuse γ-ray emission, we decided not to include Supernova Remnants (SNRs) because they could only increase our estimate and increase its uncertainty, due to the lack of a well-established description of CR production in SNRs expanding inside a star cluster. Our choice is justified a posteriori by our finding that even without inclusion of SNRs, and hence in the most conservative case, YMSCs can provide a sizeable contribution to the diffuse γ-ray emission.

The paper is structured as follows: in Sect. 2, we describe the method for the simulation of the mock star population in each cluster and the recipes used to model the stellar wind physics; in Sect. 3, we illustrate the adopted approach for the simulation of the synthetic YMSC population; in Sect. 4, we summarize the model of CR acceleration at the wind TS developed by Morlino et al. (2021) and the procedure used to compute the γ- ray emission; in Sect. 5, we compare the diffuse emission from the Galactic population with available data in the literature, and we then discuss the obtained results. Finally, we present our conclusions in Sect. 6.

2 Simulation of a mock stellar population and estimation of the cluster wind parameters

2.1 Generating the YMSCs stellar populations

For each YMSC, the mock stellar population is created based on two fundamental parameters: the cluster age (t_sc) and the number of stars (N_★) at the time the cluster was formed, with N_★ directly linked to the cluster mass (see Appendix A). We build our populations of stars in two steps: first, a population of N_★ stars is extracted by randomly sampling the initial stellar mass function (IMF); afterward, based on the cluster age, all the stars that are expected to have exploded as supernovae are removed from the population.

We here consider the IMF provided by Kroupa (2001). When sampling the IMF, we fix the minimum and maximum stellar masses that can be generated to M_★,min = 0.08 M_⊙ and M_★,max = 150 M_⊙. M_★,min is related to the minimum theoretical mass to support nuclear burning (Carroll & Ostlie 1996), while M_★,max is the maximum observed stellar mass in clusters (Zinnecker & Yorke 2007)¹.

In order to remove all the stars that have exploded as supernovae from the population of stars in a cluster of given age, we adopt the following criterion: a star of a specific mass M_★ will leave the main sequence, and subsequently explode as a supernova, at a turn-off time (t_to) given by the following (Buzzoni 2002): $${\log }_{10}\left(\frac{t_{\text{to}}}{1\text{yr}}\right)=0.825{\log }_{10}^2\left(\frac{M_{\star }}{120{\text{M}}_{\odot }}\right)+\mathrm{6.43.}$$(1)

We do not consider any post-main-sequence evolution (with the exception of stars in the Wolf Rayet phase, that will be investigated separately below): namely, all the stars with t_to < t_sc are removed from the cluster. We remind that, as anticipated in Section 1, we neglect the energy released by the supernova explosions in the calculations of the CR power.

As an additional case of investigation, we build a mock stellar population including the presence of Wolf-Rayet (WR) stars. To this purpose, we assume that all stars with M_★ > 25 M_⊙ undergo a WR phase for ∼0.3 Myr after t_to (Rosslowe & Crowther 2015), with the choice of the minimum WR mass motivated by the results of stellar evolution codes, showing that the minimum mass to develop a WR is in the range 22–37 M_⊙ (Eldridge & Vink 2006). Lower mass WR stars might in fact be born as a result of binary evolution. Those stars may provide a sizeable contribution to the wind power, but mostly at late times. In fact, this channel is important for stars in the 10–20 M_⊙ range which leave the main sequence at ages ≳10 Myr, hence this effect is bound to have a minor impact for the clusters considered in this work. In any case, neglect of this additional contribution can only lead to underestimate the wind kinetic luminosity, and hence make our estimate of the YMSC contribution to the diffuse γ-ray background a more conservative lower limit.

2.2 Modelling stellar winds

The physics of stellar winds is modelled using a purely empirical approach. The mass loss rate of a given main sequence star (Ṁ_★) is calculated using the relation provided by Nieuwenhuijzen & de Jager (1990): $${\dot{M}}_{\star }={10}^{-14.02}{\left(\frac{L_{\star }}{{\text{L}}_{\odot }}\right)}^{1.24}{\left(\frac{M_{\star }}{{\text{M}}_{\odot }}\right)}^{0.16}{\left(\frac{R_{\star }}{{\text{R}}_{\odot }}\right)}^{0.81}\frac{{\text{M}}_{\odot }}{\text{yr}}.$$(2)

The kinetic luminosity of the stellar winds (L_★,w) is defined as: $$L_{\star ,\text{w}}=\frac{1}{2}{\dot{M}}_{\star }{v}_{\text{w}}^2=\frac{1}{2}{\dot{M}}_{\star }\left\{C{\left(T_{\text{eff}}\right)}^2\left[\frac{2G{M}_{\star }\left(1-L_{\star }/L_{\text{Edd}}\right)}{R_{\star }}\right]\right\},$$(3)

where R_★ is the stellar radius, L_★ is the stellar bolometric luminosity, L_Edd is the Eddington luminosity and the term in braces represents the wind speed squared (Kudritzki & Puls 2000). The coefficient C(T_eff) depends on the stellar effective temperature T_eff and is inferred from observations (Kudritzki & Puls 2000). The stellar effective temperature is calculated using Boltzmann’s law $$T_{\text{eff}}={\left[\frac{L_{\star }}{4\pi R_{\star }^2{\sigma }_{\text{b}}}\right]}^{1/4},$$(4)

where σ_b is the Boltzmann constant. The stellar radii and luminosities are assigned based on empirical relations. For stellar radii, we use the relation proposed by Demircan & Kahraman (1991): $R_{\star }=0.85{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{0.67}{\text{R}}_{\odot }$. For stellar luminosities we use the broken power-law relation in Equation (B.4) of Menchiari et al. (2024).

When considering WR stars, we compute the mass loss rate according to Nugis & Lamers (2000): $${\dot{M}}_{\star ,\text{WR}}={10}^{-11.0}{\left(\frac{L_{\star ,\text{WR}}}{L_{\odot }}\right)}^{1.29}{\left(\frac{Y_{\text{WR}}}{Y_{\odot }}\right)}^{1.73}{\left(\frac{Z_{\text{WR}}}{Z_{\odot }}\right)}^{0.47}\frac{{\text{M}}_{\odot }}{\text{yr}},$$(5)

where Y_WR and Z_WR are the helium fraction and metallicity of WR stars respectively (both normalized to solar values). L_★,WR is the bolometric luminosity of the WR, which can be calculated using an ad hoc empirical mass-luminosity relation (see Equation (3) by Schaerer & Maeder 1992). We note that since L★,_WR depends on the WR mass, it can drastically change in time due to severe mass losses associated to the wind. We account for this dependence by considering the time evolution of the stellar mass in the evaluation of the mass loss rate, that is, Ṁ = Ṁ [L_★,WR(M_★(t))]. Finally, the kinetic luminosity of the WR wind is calculated using Equation (3), but considering a constant wind speed of 2000 km s⁻¹. In general, the wind speed in WR stars is not constant and depends on the type of WR. For example, WC and WO stars have a terminal wind speed in the range ~1000–3000 km⁻¹ s⁻¹ (Niedzielski & Skorzynski 2002; Crowther 2007), while WN stars have wind velocities spanning ~1000–2000 km⁻¹ s⁻¹ (Crowther 2007). Once the wind luminosity and mass loss rate of every i–th star of the cluster are known, the collective cluster wind mass loss rate (Ṁ) and luminosity (L_w) are obtained using mass and momentum conservation.

3 Generating synthetic population of Galactic YMSCs

In order to simulate the Galactic YMSCs, we start from the cluster formation rate, $${\xi }_{\text{sc}}\left(M_{\text{sc}},t,r,\theta \right)\equiv \frac{d{N}_{\text{sc}}}{d{M}_{\text{sc}}dtrdrd\theta },$$(6)

which gives the number of clusters born at time t per unit mass, time and surface. As commonly done in the literature, we here assume that the dependence of ξ_sc on space, time and mass, can be factorized, so that ξ_sc(M_sc, t, r, θ) = ψ(t) f (M_sc) ρ(r, θ), where ψ(t), f (M_sc) and ρ(r, θ) are, in order, the total cluster formation rate, the cluster initial mass function (CIMF) and the cluster spatial distribution.

3.1 Mass distribution of star clusters

The CIMF can be inferred from the from observation of the local SC population. In this regard, we rely on the work by Piskunov et al. (2018, based on the Milky Way Star Cluster Survey, where the CIMF is modelled as a broken power law: for clusters with M_sc > 10³ M_⊙ , the CIMF is $f\left(M_{\text{sc}}\right)=1.2{\left(M_{\text{sc}}/M_{\odot }\right)}^{-1.54}{\text{M}}_{\odot }^{-1}$. We neglect star clusters with mass smaller than 10³ M_⊙ because the integrated wind kinetic luminosity of those clusters is < 1% of the entire YMSC population. The maximum cluster mass is derived from observations (Piskunov et al. 2018), and is M_sc,max = 6.3 × 10⁴ M_⊙, corresponding to the mass of Westerlund 1. Notice that in generating the mock SC distribution, the number of stars is used in place of the cluster mass, see Appendix A.

Fig. 1 Wind power as a function of time for a YMSC with 103 M⊙ (upper plot) and 104 M⊙ (lower plot) is shown as a solid line. The shaded region shows how the wind luminosity can statistically vary for different random sampling of the initial stellar mass function. The limits of the shaded region are calculated as the 10% and 90% percentile out of 100 possible samplings of the initial stellar mass function.

3.2 Age distribution of star clusters

Similarly to the CIMF, the star clusters age distribution can be retrieved from the local cluster population. As we are interested in YMSCs, we consider only star clusters with age between t_sc,min = 0 and t_sc,max = 10 Myr, because the wind luminosity drops substantially at later times. This can be seen in Figure 1, where we show the wind power as a function of time for two YMSCs with masses of 10³ and 10⁴ M_⊙. Since the wind power depends on the number of massive stars, we also show the range in which it can vary, calculated as the 10% and 90% percentile out of 100 possible samplings of the stellar IMF. After accounting for the effect of cluster evolution, the age distribution of stellar clusters is found to be approximately uniform in the last few tens of megayears (Piskunov et al. 2018). This implies a constant cluster formation rate in the time interval we are interested in. This rate can be estimated from the star formation rate of local young stellar clusters within 1 kpc as measured by Bonatto & Bica (2011). Assuming the CIMF from Piskunov et al. (2018) (underfilling case) leads to an average local cluster formation rate of ψ = 1.8 Myr⁻¹ kpc⁻² (Menchiari 2023).

3.3 Spatial distribution of star clusters

The spatial distribution of young clusters is strictly linked to the value of the cluster formation rate across the Milky Way. This is well traced by the position of the giant molecular clouds (GMCs), which in turn follows the spiral pattern of the Galaxy. We spatially allocate the synthetic stellar clusters following a two-step procedure:

1. We generate stellar clusters with a galactocentric radial distribution following that of GMCs (under the assumption of an isotropic angular distribution) and assuming an exponential altitude distribution similar to that of the observed gas profile.

2. Based on its radial and angular position, we associate each synthetic cluster to a specific Galactic structure, i.e. spiral arm, galactic bar, etc. We detail this in the following.

For the radial distribution of GMCs and the modelling of the Milky Way spiral structure, we rely on the work by Hou & Han (2014), who also supply a complete catalogue of GMCs where, for each cloud, the galactic coordinates and kinematic velocity are reported. For a large portion of the GMCs, the kinematic ambiguity is resolved. We compute the kinematic distance of every cloud with disentangled kinematic ambiguity using the code developed by Wenger et al. (2018). When doing so, we adopt the classical rotation curve based method, using the state-of-the-art curve provided by Reid et al. (2019). A small fraction of GMCs in the catalog have also distances calculated using more reliable methods such as parallax. We hence adopt those estimates when available. Once the position of the GMCs is known, we can compute their radial distribution, which we express in terms of the surface mass density of molecular gas $\left({\mathrm{\Sigma }}_{H_2}^{GMC}\right)$. To calculate the surface mass density, we average the total GMCs mass in 18 rings with constant width of 1 kpc and radii spanning in the interval 0–18 kpc, and centred on the Galactic Centre. At this point we can formally define the radial distribution of YMSCs, normalized at the Sun position, as: $$\rho (r)=\frac{{\mathrm{\Sigma }}_{H_2}^{GMC}(r)}{{\mathrm{\Sigma }}_{H_2}^{GMC}(r=8.5\text{kpc})}.$$(7)

Notice that the normalization of ξ_sc at the Sun position is given by the local observed cluster formation rate ψ.

Once the radial distribution is known, we allocate stellar clusters following the Milky Way observed morphology. To reproduce the Milky Way structure we use the four logarithmic arms model provided by Hou & Han (2014), defined as: $$\ln \left(\frac{r}{R_i}\right)=\left(\frac{\theta -{\theta }_i}{1\mathrm{rad}}\right)\tan {\Psi }_i,$$(8)

where R_i, θ_i, and Ψ_i are parameters inferred from the fit to the position of HII regions, masers and GMCs. We report their values for each arms and for the local spur in Table 1. On top of the spiral structure, we also take into account the structure of the innermost region. Here we consider the presence of the Galactic bar and the Near 3 kpc and Far 3 kpc arms. The first is modelled as an ellipse with semi-major axis of 3.3 kpc and an aspect ratio of 4:10 (Churchwell et al. 2009). The Near/Far-3-kpc Arms are modelled using an elliptical annulus, with a semi-major axis of 4.1 kpc, and an aspect ratio of 0.54 (a semi-minor axis of 2.2 kpc) (Green et al. 2011).

To allocate the synthetic YMSCs in the Galaxy we use the following procedure: first, for every j-th cluster, we start by randomly generating its radial (r_j) and angular (θ_j) coordinates. Radial distances are extracted considering the probability distribution ρ(r) given by Equation (7), while the angular coordinate is chosen assuming a uniform distribution. Afterward, we randomly choose the structure (i.e. spiral arm, local spur or Near/Far-3-kpc Arms and Galactic Bar) to associate the cluster with, among the ones that satisfy the following criteria:

• Spiral Arm i = 1 to 4, if r_j > R_i;

• Local spur if 7.59 < r_j/kpc < 9.17 and 50.6° < θ_j < 110°;

• Spiral Arm i = 1, 3, 4 or Near/Far-3-kpc Arms and Galactic Bar, if R_i < r_j < 4.29 kpc;

• Near/Far-3-kpc Arms and Galactic Bar, if r_j < 3.27 kpc.

When more structures satisfy the criteria, they have the same probability to be chosen. Finally, when a cluster is associated with the Near/Far-3-kpc Arms and the Galactic Bar structures, the ultimate choice between which of the two structures should be matched is determined by a criterion of minimum distance.

Note that the angular coordinate θ_j is used only to check whether the YMSC should be associated with the Local Spur. Once the cluster is placed in a given arm, the coordinate θ_j is recalculated by inverting Equation (8): $${\theta }_j={\left(\frac{\tan {\Psi }_i}{1\mathrm{rad}}\right)}^{-1}\ln \left(\frac{r_j}{R_i}\right)+{\theta }_i.$$(9)

When a cluster is associated with the Galactic Bar we check whether its position is actually located within the Bar. If not, θ_j is randomly extracted until its position falls within the Bar. Finally, when a cluster is associated with the Near/Far-3-kpc Arm, we change its coordinates to the nearest point belonging to the ellipse defining these structures.

Once all clusters are associated to a structure, we proceed to perturb their positions according to a Gaussian distribution, following the fundamental concept that both Spiral Arms and the Near and Far structure possess an intrinsic thickness. The width of the spiral arm increases with galactocentric distance, so that the probability 𝒫 of finding a source displaced from the centre of the arm by Δr_j, is (Faucher-Giguère & Kaspi 2006): $$\mathcal{P}\left(\Delta r_j\right)=\frac{1}{2\pi \left(0.07r_j\right)}\exp \left[-\frac{\Delta r_j^2}{2{\left(0.07r_j\right)}^2}\right].$$(10)

As for the clusters belonging to the Near and Far arms, we extract from a Gaussian probability distribution the Galactic x and y coordinates with a spread of σ_x = σ_y = 0.1667 kpc, such that the probability at 3σ returns a scattering compatible with the observed radial thickness of 0.5 kpc (Green et al. 2011). Notice that the thickness of a spiral arm in YMSCs may be slightly different from that inferred for pulsars, however the result on the diffuse emission depends very weakly on this parameter.

At last, we generate the vertical coordinate (ɀ) following the observed gas distribution, which is an exponential profile with a characteristic spread of 100 pc (Strong et al. 2000): $$\rho (ɀ)=\exp \left(-\frac{ɀ}{100\text{pc}}\right).$$(11)

With the cluster formation rate ξ_sc defined it is possible to finally build the synthetic SC population. We find a total number of 747 clusters with average values for the wind luminosity and the mass loss rate equal to ⟨L_w⟩ ≈ 3 × 10³⁶ erg s⁻¹ and ⟨Ṁ⟩ ≈ 10⁻⁶ M_⊙ yr⁻¹, respectively. Figure 2 shows a possible realization of the Milky Way YMSCs population.

Fig. 2 Single realization of the Milky Way population of YMSCs. Top panel: face-on view of the Milky Way. The size and colour of each YMSC refer to its mass and age respectively, according to the scales indicated in the panel. Notice that the size distribution is continuous and the associations between mass and size reported in the panel are only for reference. The grey regions represents all the Galactic structures (spiral arms, Local Spur and Galactic bar) employed in the simulation of the Milky Way cluster population. The inset shows a zoomed view of the Milky Way center. The location of the Central Molecular Zone (CMZ) is also reported for reference. The star marks the position of the Solar System. Middle and bottom panels: edge on view of the Milky Way towards ROI1 (upper panel) and ROI2 (lower panel). Grey circles represent the projected sizes of the wind-blown bubbles surrounding YMSCs. The striped regions correspond to the 𝒜 ∪ ℬ masks employed in our analysis for this specific realization of the Milky Way (see Appendix B).

4 Cosmic ray acceleration at the collective wind termination shock of YMSCs

When the average distance between stars in a cluster is less than the size of individual stellar wind TS, the winds from massive stars are expected to combine to create a collective cluster wind. This scenario is favoured in young clusters as a direct consequence of mass segregation at birth (Lada & Lada 2003), namely the fact that the most massive stars (hence those producing the most powerful winds) forms at the very core of the SC. In our mock SC population, the star distribution is not simulated. We simply assume that all YMSCs are compact enough to generate a collective TS. As the wind material is shocked and heated, it expands adiabatically in the interstellar medium (ISM) generating large bubble-like structures similar to those observed close to single massive stars (Weaver et al. 1977). A model for particle acceleration in these systems was presented by Morlino et al. (2021) based on the idea that particles are accelerated at the collective wind TS via the diffusive shock acceleration mechanism, and subsequently escape from the acceleration site experiencing a combination of advection and diffusion within the hot bubble, until reaching its border. From there, CRs are free to escape in the unperturbed ISM.

Given the highly turbulent nature of the shocked wind, particle diffusion is likely suppressed, resulting in a long confinement time within the bubble where the majority of γ-ray emission is produced through hadronic interactions. The CR distribution in the bubble is (Morlino et al. 2021): $$f(r,p)=f_{\text{ts}}(p)\mathcal{G}\left(D_2,R_{\text{b}}\right),$$(12)

where the distribution function at the TS is $$f_{\text{ts}}(p)=\mathcal{K}\left({ϵ}_{\text{cr}},L_{\text{w}},\dot{M}\right){\left(\frac{p}{m_{\text{p}}c}\right)}^s\left[1+a_1{\left(\frac{p}{p_{\text{max}}}\right)}^{a_2}\right]e^{-a_3{\left(\frac{p}{p_{\max }}\right)}^{a_4}}.$$(13)

The parameters a_1,…,4 are determined by the energy dependence of the diffusion coefficient, which is ultimately related to the magnetohydrodynamic turbulence spectrum. In this work, we consider three different diffusion regimes: Kolmogorov, Kraich- nan and Bohm-like diffusion. The coefficients a_1,…,4 relevant for each of these three diffusion regimes are reported in Table 2. This coefficients are derived in Menchiari et al. (2024) and are valid for any YMSCs. The parameter s is the spectral injection slope of accelerated particles. We fix for all clusters s = −4.2. This value results from the modelling of the Cygnus OB2 emission performed by Menchiari et al. (2024), but a similar spectrum is also deduced from γ-ray measurements of other YMSCs (see e.g. Aharonian et al. 2019, 2022; Yang & Aharonian 2017; Liu et al. 2022; Mitchell et al. 2024) for which an energy spectrum of 2.3–2.4 is found. 𝒦 is a normalization factor that depends on the cluster wind luminosity, mass loss rate and CR acceleration efficiency (ϵ_cr) (see Equation (12) in Menchiari et al. 2024). The function 𝒢 (see Equations (10b) and (11) in Menchiari et al. 2024) depends on the diffusion coefficient in the bubble and the size of the bubble, which reads (Weaver et al. 1977): $R_{\text{b}}=0.798L_{\text{w}}^{1/5}{\rho }_0^{-1/5}{t}_{\text{sc}}^{3/5}$, where ρ₀ is the average density close to the YMSC that we assume to be 10m_p cm⁻³ (m_p is the mass of the proton). When calculating the diffusion coefficient in the system, we assume, for the Kolmogorov and Kraichnan cases, that the turbulence is injected at a characteristic length scale of ~ 1 pc, roughly corresponding to the size of the cluster core. In the Bohm case, we assume, instead, that the injection of turbulence occurs at all scales. The parameter p_max is the maximum momentum determined based on the Hillas (confinement) criterion: D_w(p_max)/v_w = R_ts, where D_w is the diffusion coefficient in the cold cluster wind, v_w = (2L_w/Ṁ)^1/2 is the collective cluster wind speed and $R_{\text{ts}}=0.791{\dot{M}}^{1/2}{v}_{\text{w}}^{1/2}{L}_{\text{w}}^{-1/5}{\rho }_0^{-3/10}{t}_{\text{sc}}^{2/5}$ is the size of the TS (Weaver et al. 1977).

The γ-ray flux from a single YMSC due to hadronic interactions is calculated as follows: $${\varphi }_{\gamma }\left(E_{\gamma }\right)=\frac{c{n}_0}{d^2}{{\displaystyle \iint }}_{R_{\text{ts}}}^{R_{\text{b}}}{r}^{2}f\left(r,E_p\right)\frac{d\sigma \left(E_p,E_{\gamma }\right)}{d{E}_p}drd{E}_p,$$(14)

where d is the distance from Earth and n₀ = 10 cm⁻³ is the average target density, which is expected to be close to the external ISM density, as a result of the expected fragmentation of the swept-up shell (Blasi & Morlino 2024). The differential cross section for γ-ray production, dσ/dE_p, is taken from Kafexhiu et al. (2014) and based on SYBILL.

Since, the γ-ray emissivity integrated along the line of sight is expected to depend on the projected angular radius only weakly (see Menchiari et al. 2024), in order to speed up the calculation, we use, as a spatial model for the sources, a uniform disk of size equal to that of the projected wind bubble.

Finally, we note that the γ-ray emission depends almost linearly on the bubble volume due to the swept-up mass. Hence, a correct estimate of the bubble size is important. The adiabatic model by Weaver et al. (1977) used here overestimates the bubble size because it does not account for the ISM external pressure nor for the radiative losses. While the former has a very minor impact, the latter may be relevant at times ≳10 Myr (Silich & Tenorio-Tagle 2013). Such a timescale may be even reduced due to the external shell fragmentation which increases the contact surface between the hot shocked wind and the cold ISM (Lancaster et al. 2021). On the other hand, when also SN explosion are taken into account, the effect of cooling will be compensated by the injection of additional power. Including all those effects is a not trivial task because it requires a full time dependent calculation able to account for the full 3D hydro-dynamical evolution of the bubble. However, at the end of Section 5, we evaluate the impact of a possible bubble size reduction by arbitrarily reducing its radius.

5 The contribution of YMSC to diffuse γ-ray emission

We now focus on two regions of the sky for which recent observations of the diffuse γ-ray emission are available in the literature, and compute the corresponding emission of our synthetic YMSCs. The first region (named ROI1), is defined for |b| ≤ 5° and 15° ≤ l ≤ 125°, while the second region (named ROI2) is defined for |b| ≤ 5°, and 125° ≤ l ≤ 235°. For both ROI1 and ROI2, flux measurements are provided by Fermi-LAT (Zhang et al. 2023) and LHAASO (Cao et al. 2023). In ROI1 additional measurements of the diffuse flux are provided by the ARGO experiment² (Bartoli et al. 2015).

5.1 Masking the contribution of resolved sources

Notice that all the data we use are obtained after removing known sources, to avoid an artificial increase of the Galactic diffuse emission. However, the procedure used for this purpose are different for each of the considered datasets. LHAASO measurements are obtained after masking sources from both the TeVCat (Wakely & Horan 2008) and the 1LHAASO catalog (Cao et al. 2024), while in the ARGO analysis only TeVCat sources have been masked³. The Fermi-LAT data results from applying the same masks used in the LHAASO analysis after removing the emission of the 4FGL sources (Abdollahi et al. 2020). This is done by employing a binned likelihood analysis to model and then subtract the contribution from point sources and extended sources of the 4FGL catalog located outside the masked regions.

To be consistent with this modus operandi, we also mask the simulated emission in the two ROIs following a similar method to that implemented in the LHAASO analysis. To do so, we use three different masks, referred to as 𝒜, ℬ, and 𝒜 ∪ ℬ (see Appendix B for additional details):

• 𝒜 removes the inner Galactic Plane (l ≤ 70°, |b| ≤ 1.5°) and the Local Arm (disk centred on [l = 73.5°, b = 0°]). This mask is designed to emulate the removal of Galactic sources, such as pulsar wind nebulae or SNRs, which are not considered in our realization of the Milky Way and that are predominantly located along the Galactic Plane, within the first and fourth Galactic Quadrants.

• ℬ masks all those positions of the sky where the significance of the emission from YMSCs at 100 TeV, evaluated as described in Appendix B, is larger than 5 times the diffuse emission. This mask is intended to remove all clusters whose emission should be bright enough to be detectable by LHAASO. We here assume that all the clusters detected by LHAASO should be also detected by Fermi-LAT. Hence, we do not mask nor remove the emission from YMSCs that should have been detected by Fermi-LAT but not by LHAASO. We comment a posteriori the implications of this assumption.

• 𝒜 ∪ ℬ is the union of 𝒜 and ℬ. The flux obtained using this mask is the most suitable for comparison with data, as 𝒜 ∪ ℬ most closely matches the masks used in the works by Cao et al. (2023) and Zhang et al. (2023).

We note that in ROI2 the mask 𝒜 is not defined, so in this region one will have ℬ ≡ 𝒜 ∪ ℬ (see Appendix B). The bottom panels in Figure 2 reports an example of the mask 𝒜 ∪ ℬ in both considered regions for a random realization of the Milky Way YMSCs population.

Fig. 3 Contribution to the diffuse γ-ray emission from a synthetic population of YMSCs in ROI1 (top row) and ROI2 (bottom row) for the three different cases of the diffusion coefficient considered: Kolmogorov (left column), Kraichnan (central column) and Bohm-like (right column). The spectrum is calculated after applying the mask 𝒜 ∪ ℬ. In each plot, the dashed line represents the median flux per energy bin after considering 100 different realisations of the Milky Way population of YMSCs, while the associated shaded region encloses the 25–75 percentile flux. The striped region instead encompasses the 25–75 percentile flux when WR stars are not included in the stellar mock populations. The violet solid lines in the right column and their associated shaded regions show the median and the 25–75 percentile diffuse γ-ray flux after the addition of emission from the CR sea.

5.2 Comparison with the observed galactic diffuse γ-ray emission

To account for the uncertainty associated with the stochastic nature of the generation of the synthetic population of YMSCs, we simulate 100 different realizations of the Milky Way and calculate, for each of them, the contribution to the diffuse γ-ray flux, adopting different masks and assumptions on the particle diffusion. Figure 3 shows the YMSCs contribution to the diffuse γ-ray emission in ROI1 (top row) and ROI2 (bottom row) for the three types of diffusion coefficient considered. The emission from YMSCs is calculated after applying 𝒜 ∪ ℬ in both ROI1 and ROI2. In all plots the dashed lines represent the median flux per energy bin over the different (100) realisations of the Milky Way. The shaded areas cover instead the range between the 25 and 75 percentile of the distribution, including (colour filled) or neglecting (striped) the contribution from WR stars.

In both ROIs, the γ-ray flux due to unresolved YMSCs is negligible when considering only the wind power from main sequence stars $\left({\varphi }_{\text{ymsc}}^{\text{MS}}\right)$. On the contrary, when WR stars are added to the stellar populations $\left({\varphi }_{\text{ymsc}}^{\text{WR}}\right)$ the contribution from YMSCs can account in some cases for a substantial fraction of the measured flux. To better visualize the relative importance of the emission from YMSCs, we also show the total expected diffuse γ-ray flux $\left({\varphi }_{\text{tot}}^{\text{WR}}\right)$ when the contribution from the bulk of diffuse CRs (i.e. the CR sea) is considered: ${\varphi }_{\text{tot}}^{\text{WR}}=F_{\gamma }+{\varphi }_{\text{ymsc}}^{\text{WR}}$, where F_γ is the flux from the CR sea (see Appendix C).

More precisely, in ROI1, the contribution of YMSCs in the Kolmogorov case can explain the diffuse γ-ray emission from ∼10 GeV to ∼1 TeV. In the Kraichnan scenario, however, the contribution of YMSCs is higher and the total γ-ray emission exceeds that observed from ∼10 GeV to ∼10 TeV by ∼20%. This implies that in the Kraichnan case, the contribution from YMSCs is able to totally account for the excess observed in the diffuse galactic emission if the product between ϵ_cr and n₀ is decreased by a factor ∼0.8. On the other hand, in the Bohm case, the contribution from YMSCs to the diffuse emission is negligible. This is because, when Bohm diffusion is adopted, many more clusters have γ-ray emission above the LHAASO threshold and are hence removed by applying the mask ℬ (see Appendix B).

In ROI2 the contribution from YMSCs is on average negligible, regardless of the presence of WR stars, with the exception of the Kraichnan case, where YMSCs could significantly contribute at ∼100 GeV. However, no solid conclusion can be drawn due to the large statistical uncertainties in this region. In fact, due to the low content of YMSCs in the outer Galaxy, some realisations have a much larger γ-ray flux. However, it is likely that the knowledge of the population of YMSCs in the outer region of the Galaxy is complete or quasi-complete. Consequently, one could consider to rely directly on observed clusters from Gaia (see, e.g. Celli et al. 2024; Mitchell et al. 2024), rather than using our approach.

These results must be considered with a few caveats in mind. First, the diffuse emission below ∼100 GeV is likely to be slightly overestimated, as we did not subtract the emission from those YMSCs that would have been detected by Fermi-LAT but not by LHAASO. In any case, we do not expect this choice to lead to a large error because the most powerful clusters are subtracted anyway. Second, as mentioned in Section 4, because of radiation losses, the size of the bubble could be smaller, reducing the γ-ray emission of each cluster. On average a radius reduction between 20% and 50% may be expected. If we assume that every bubble is systematically reduced by 20% we found that the overall flux is reduced by roughly the same factor. On the other hand, for a reduction of 50%, the flux is reduced by a factor ∼8, and the emission from YMSCs is drastically suppressed. However, the latter case is somehow unrealistic for two reasons: first, radiation losses are in general expected to appear after ∼10 Myr (Silich & Tenorio-Tagle 2013), and second, one should also consider the re-heating of the gas induced by supernova explosion in the cluster, which can compensate for the energy losses. Finally, we note that adding the contribution of SNe the total predicted flux would likely overshoot the data for the current values of CR acceleration efficiency and gas density, especially in the Kraich- nan scenario. We expect that consistency with the data could still be recovered by adjusting the product between ϵ_cr and n₀ , while still keeping both parameters in a reasonable range. However, fitting the observed spectrum is beyond the scope of this work, as it would also require the inclusion of additional populations of unresolved galactic γ-ray sources.

To summarize, we assess that non-resolved YMSCs do provide a considerable contribution to the diffuse galactic γ-ray emission in the inner region of the Galactic plane if the particle diffusion coefficient follows a Kolmogorov-like or Kraichnan- like spectral energy dependence. Such a conclusion is strengthened by the fact that our estimate neglects the emission coming from particles accelerated at SNRs.

6 Conclusions

In recent years, the improving sensitivity of γ-ray telescopes has led to the identification of several YMSCs as γ-ray sources. Nonetheless, nowadays, the detection of individual YMSCs still remains a challenging task due to their low surface brightness. This implies that most of the γ-ray flux from YMSCs is likely detected as an unresolved contribution to the diffuse Galactic γ-ray emission. Indeed, recent analyses of the diffuse γ-ray emission have shown a prominent excess, likely related to unresolved sources, ranging from a few tens of gigaelectronvots up to hundreds of teraelectronvolt. In this work, we estimated a lower limit to the contribution of YMSCs to the diffuse emission by simulating the entire Galactic populations of YMSCs, starting from observed properties of local clusters. We modelled the γ-ray spectrum of each synthetic YMSC considering a purely hadronic scenario and assuming that particle acceleration occurs at the collective cluster wind termination shock. We found that the contribution of the cluster population to the diffuse emission is sizeable above few tens of gigaelectronvots in the inner region of the Milky Way when the particle transport in the wind-blown bubble is governed by Kolmogorov or Kraichnan type diffusion. In the outer region of the Milky Way we could not draw any solid conclusion due to the large statistical uncertainties related to the realisation of the YMSC population. A more reliable estimate in this region should be made by directly using existing YMSC surveys.

The contribution of Wolf-Rayet stars turns out to dominate the collective wind power and hence the CR production, so that uncertainties in the modelling of these stars are the main caveat to our conclusions. Another source of uncertainty comes from the neglect of CRs accelerated by SNRs within the cluster. These are likely to further enhance the predicted γ-ray flux. However, at the moment, some important aspects concerning the evolution of SNRs inside a bubble and the associated particle acceleration are still unclear and require a dedicated effort to be fully understood. All the effects we neglected, if properly taken into account, can only lead to estimate a larger contribution of YMSCs to the diffuse γ-ray background, possibly reinforcing our conclusions: YMSCs provide a non-negligible contribution to the diffuse γ- ray emission of the Galaxy, and their contribution must be taken into account in any future realistic model of the Galactic diffuse emission.

Acknowledgements

We thank Pr. Ruo-Yu Liu for providing the detector exposure on behalf of the LHAASO collaboration. This work was partially funded by the European Union – NextGenerationEU RRF M4C2 1.1 under grant PRIN-MUR 2022TJW4EJ SM and GM are partially supported by the INAF Theory Grant 2022 “Star Clusters As Cosmic Ray Factories” and by the INAF Mini Grant 2023 "Probing Young Massive Stellar Cluster as Cosmic Ray Factories".

Appendix A Generation of the mock stellar population

When generating the mock stellar population in a cluster, it is operationally more straightforward to use the number of stars rather than the cluster mass. These two quantities are related by: $$N_{\star }=\Lambda M_{\text{sc}},$$(A.1)

with the parameter λ defined as: $$\Lambda =\frac{{\displaystyle {\int }_{M_{\star ,\min }}^{M_{\star }\max }{f}_{\star }}\left(M_{\star }\right)d{M}_{\star }}{{\displaystyle {\int }_{M_{\star ,\min }}^{M_{\star ,\max }}{M}_{\star }}{f}_{\star }\left(M_{\star }\right)d{M}_{\star }},$$(A.2)

where M_⋆ is the stellar mass, f_⋆ is the initial stellar mass function, and M_⋆,min and M_⋆,max are, respectively, the minimum and maximum stellar masses that can be generated in a cluster.

Because the number of stars in a cluster is the starting ingredient to generate the mock stellar population, for convenience, we operationally build up the cluster population in the Galaxy not as a function of the cluster mass, but rather as a function of N_⋆. This implies that the distribution function used to $f^{\prime }\left(N_{\star }\right)\equiv \frac{d{N}_{\text{sc}}}{d{N}_{\star }}$ and is related to f (M_sc) as a rewriting of the initial cluster mass function in terms of N_⋆ using Equation A.1 $$f^{\prime }\left(N_{\star }\right)=f\left(M_{\text{sc}}\right)\frac{d{M}_{\text{sc}}}{d{N}_{\star }}=1.2{\left(\Lambda M_{\odot }\right)}^{0.54}{N}_{\star }^{-1.54}.$$(A.3):

By doing so, the maximum and minimum cluster mass directly translate to a maximum and minimum number of stars, that is N_⋆,min = λM_sc,min and N_⋆,max = λM_sc,max.

Appendix B Definition of masks

The measurements of the Galactic diffuse emission in both ROIs are provided after masking known very-high-energy γ-ray sources. Therefore, in order to compare the emission from simulated YMSCs with the available data, our simulated sources should be masked following a similar criterion. The masks used by Zhang et al. (2023) and Cao et al. (2023) (M1_L for ROI1 and M2_L for ROI2 in the following) are built by removing the sources reported in the TeVcat plus the ones recently discovered by the KM2A detector of LHAASO (Cao et al. 2024). These masks cannot be applied to our simulated maps by simply using their spatial coordinates, because our sources are in different positions from the real ones and different from one realization to the other. The method we implemented to mimic the masks used by Cao et al. (2023) is described in the following.

We used three different kinds of masks. The first mask, A, is designed to remove regions of the Galactic Plane where, on average, most of the galactic sources are located. Following the shape of MI_L, one can easily see that, for l < 70°, a large fraction of the Galactic Plane is masked. We hence define 𝒜 so as to remove the zone of the sky defined between 15° < l < 70° and |b| < 1.5°. We also remove the region corresponding to the projected size of the Local Spiral arm, as an over-density of sources is expected there. Such a region corresponds to a circular area centred on (l = 73.5°, b = 0°) with radius 6.5°. Notice that in ROI2 a sparse distribution of the galactic sources is expected, as the outer region of the Milky Way is less crowded. Consequently, there is not a specific zone of the plane to mask. Hence, we define this specific type of mask only for ROI1.

The second type of masks, B, is defined so as to remove the emission from those YMSCs that would be detectable by the LHAASO experiment. To do so, we mask the positions of the sky where the significance of the surface brightness at 100 TeV (S₁₀₀) exceeds by 5 σ the brightness of the measured galactic diffuse emission, considered here as the nominal background. We use the emission at 100 TeV as it correspond to the energy where LHAASO has the best sensitivity (Cao et al. 2019). Hereafter, the subscript "100" indicates that the quantity is evaluated at 100 TeV. At each position of the sky belonging to the ROIs with coordinates (l₀, b₀), the significance of the emission is calculated as: $${\mathcal{S}}_{100}\left(l_0,b_0\right)={\displaystyle \iint \mathrm{\Theta }}\left(l_0,b_0,r_{100}^{\text{psf}}\right)\frac{N_{100}^{\text{sc}}(l,b)}{\sqrt{N_{100}^{\text{gde}}(l,b)}}dldb.$$(B.1)

𝒜 given position of the sky is then masked if 𝓢₁₀₀ > 5. In Equation B.1, $\Theta \left(l_0,b_0,r_{100}^{\text{psf}}\right)$ is a 2D top hat function centred on (l₀, b₀) with radius equal to the LHAASO point spread function $\left(r_{100}^{\text{psf}}={0.2}^{°}\right)$ (Cao et al. 2019), while $N_{100}^{\text{sc}}(l,b)$ and $N_{100}^{\text{gde}}(l,b)$ are the expected counts from the YMSCs and the galactic diffuse γ-ray emission, respectively. The latter quantities are calculated as: $$N_{100}^{\text{sc}}(l,b)={\varphi }_{100}^{\text{sc}}(l,b){ϵ}_{100}(l,b)\mathrm{\Delta }{E}_{100}$$(B.2a) $$N_{100}^{\text{gde}}(l,b)={\varphi }_{100}^{\text{bkg}}{ϵ}_{100}(l,b)\mathrm{\Delta }{E}_{100},$$(B.2b)

where ${\varphi }_{100}^{\text{sc}}(l,b)$ is the flux at the position (l, b) only due to YMSCs, while ${\varphi }_{100}^{\text{gde\hspace{0.17em}}}$ is the measured value of the diffuse galactic flux measured by LHAASO (Cao et al. 2023). Note that by doing so, we are assuming that the observed Galactic diffuse γ-ray emission is constant over the entire ROI. The value ∆E₁₀₀ is the energy bin width corresponding to the measured flux point at 100 TeV (Cao et al. 2023). Finally, ϵ₁₀₀(l, b) is the detector exposure at the position (l, b), which is calculated as: $${ϵ}_{100}(l,b)=A_{\text{eff},100}{t}_{\text{exp}}(l,b),$$(B.3)

where A_eff,100 is the KM2A effective area for γ-like events (Cao et al. 2019) and t_exp(l, b) is the exposure time of a given location in the sky. For the latter we consider that the detector array has operated also during the construction phase⁴, so that $$t_{\exp }(l,b)=\left[\left(\frac{1}{2}302+\frac{3}{4}219+423\right)\text{\hspace{0.17em}days\hspace{0.17em}}\right]\times t_{\text{hpd}}(l,b),$$(B.4)

where t_hpd(l, b) is the amount of time per day a given sky position is found at zenith angles (ζ) below 50°, which can be evaluated as: $$t_{\text{hpd}}(l,b)={\displaystyle {\int }_0^{24\text{\hspace{0.17em}h}}\mathscr{H}}\left({50}^{°}-\zeta (l,b,t)\right)\cos [\zeta (l,b,t)]dt,$$(B.5)

where H (50° – ζ(l, b, t)) is the Heaviside function, and $$\cos [\zeta (l,b,t)]=\sin \left({\delta }_{\text{L}}\right)\sin [\delta (l,b)]+\cos \left({\delta }_{\text{L}}\right)\cos [\delta (l,b)]\cos \left[{15}^{°}\left(\frac{t}{\text{\hspace{0.17em}h}}\right)\right].$$(B.6)

In Equation (B.6) the angles δ_L and δ(l, b) are the equatorial declination of the LHAASO site and the given position of the sky, respectively. Finally, the last type of mask, 𝒜 ∪ ℬ, is obtained by merging 𝒜 with ℬ. In ROI1 the 𝒜 ∪ ℬ mask is the one that best reproduces MI_L. Figure B.1 shows how the contribution to the diffuse gamma-ray flux in ROI1 changes after applying the 𝒜 (left column), ℬ (central column) and 𝒜 ∪ ℬ (right column) masks.

Figure B.2 shows an example of how 𝒜, ℬ, and 𝒜 ∪ ℬ compare with Ml_L. This specific configuration refers to the Galaxy realization showed in Figure 2 with the assumption of Kraichnan diffusion and inclusion of the WR contribution. Using all the 100 different realisations of the Milky Way, we can quantitatively estimate the difference between our masks and the original ones by evaluating the median ratio of the non-masked areas. For the 𝒜 ∪ ℬ mask we found this ratio to be 1.22, 1.19 and 1.1 for the Kolmogorov, Kraichnan and Bohm cases, respectively. The same ratios calculated for the ROI2 are 1.24, 1.24 and 1.22. Those numbers are estimated including WR stars in the stellar population. We found no significant difference when WRs are excluded. In summary, our sky regions is between 10% and 20% larger than the one considered by Cao et al. (2023).

The same method employed to generate the mask ℬ can also be used to estimate the number of YMSCs that should have been detected by LHAASO. This can be seen in Figure B.3 were the distribution of the number of detected clusters over the 100 different realisations is shown. We note that this number of detections must be considered as an order of magnitude estimate, as the detectability of a cluster should be evaluated after considering the instrument sensitivity over the entire energy range. Interestingly, the number of detections for the Kolmogorov and Kraichnan cases is consistent with the number of sources in the LHAASO catalog that can be associated to YMSCs (Mitchell et al. 2024). In contrast, the predicted median number of YMSCs in the Bohm case (~ 100 if WRs are included in the stellar populations) is of the same order of the total number of sources in the LHAASO catalog, suggesting that this scenario is likely to be too extreme.

Appendix C Diffuse emission from the CR sea

The γ-ray spectrum produced by the CR sea is calculated as in the work by Peron & Aharonian (2022), namely by assuming that the CR spectrum around the Galaxy is the same as the one locally measured at Earth. We consider the spectrum of CR protons, J(E_p), measured by AMS-02, Dampe, Kascade, and IceTop and fitted by Lipari & Vernetto (2020). To account for heavier elements, both in the CR composition and in the target, we consider a nuclear enhancement factor, ξ_N , computed following Mori (2009), by considering a solar-like composition for the target and the different measured energy spectra for the different CR elements. This implies that the nuclear enhancement factor depends on energy. We used the differential cross section given by Kafexhiu et al. (2014), choosing the SYBILL parametrization. Potential differences arising from the different choice of the cross-section parametrization are neglected in this work. The CR flux can be then computed as: $$F_{\gamma }\left(E_{\gamma }\right)=N_H\mathrm{\Delta }\mathrm{\Omega }{\xi }_N\left(E_{\gamma }\right){\displaystyle \int d}{E}_p\frac{d\sigma \left(E_{\gamma },E_p\right)}{d{E}_p}J\left(E_p\right).$$(C.1)

The density of the target, N_H, is derived from the Planck dust-opacity, τ_D , observed in the regions outside the M1_L and M2_L masks. The result is N_H = 7.85 × 10²¹ cm⁻² in ROI1 and N_H = 5.89 × 10²¹ cm⁻² in ROI2, when assuming a dust-to-H conversion factor $X_D^{-1}\equiv \frac{{\tau }_D}{N_H}=1.18\times {10}^{-26}{\text{cm}}^2$ (Planck Collaboration 2011). Finally, we consider the fact that the γ-ray emitting gas is located at different distances from us, and consequently, its emission can be affected by absorption due to the γ-γ interaction between the high-energy photons and the interstellar radiation fields. The effect of absorption is evaluated for the average Galactic radiation fields provided by Lipari & Vernetto (2018) as a function of the source distance. We here derive the gas distance by considering the measured radial velocity at each position. To do this we used the gas maps provided along with the Galprop code (Strong et al. 2009), already decomposed according to distance from the Galactic Center. This allows us to estimate what fraction of column density is located at which distance and weight for the interstellar absorption accordingly.

Fig. B.1 Contribution to the diffuse γ-ray emission from a synthetic population of YMSCs in ROI1 for the three different cases of diffusion coefficient considered: Kolmogorov (upper row), Kraichnan (central row) and Bohm-like (lower row). Each column shows the flux for the three masks implemented: 𝒜 (left column), ℬ (central column) and 𝒜 ∪ ℬ (right column). In each plot, the dashed line represents the median flux per energy bin after considering 100 different realisations of the Milky Way population of YMSCs, while the associated shaded region encloses the 25-75 percentile flux. The striped region instead encompasses the 25-75 percentile flux when WR stars are not included in the stellar mock populations. The violet solid lines in the right column and their associated shaded regions show the median and the 25-75 percentile diffuse γ-ray flux after the addition of emission from the CR sea.

Fig. B.2 Example of the three masks employed to estimate the diffuse γ-ray emission in ROI1. Upper panel: mask 𝒜 used to exclude the Local Spiral Arm and a large fraction of the inner Galactic Plane where most of the γ-ray sources are expected to be located. Central panel: mask ℬ used to remove the emission from the most luminous clusters that should have been detected by LHAASO (this example mask is the one built based for the specific source realization shown in Figure 2 and assuming the Kraichnan diffusion coefficient). Lower panel: mask 𝒜 ∪ ℬ. For comparison, in all panels, the M1L mask is shown as black contours.

Fig. B.3 Distribution of the number of synthetic YMSCs detectable by LHAASO in ROI1 (upper row) and ROI2 (bottom row) for the three diffusion coefficients considered in this analysis: Kolmogorov, Kraichnan and Bohm (from left to right). We label a YMSC as detected if the significance of its emission at 100 TeV is larger than 5σ. Colour filled (striped) histograms correspond to the scenarios where WRs are (are not) included in the mock stellar populations.

References

All Tables

All Figures

Fig. 1 Wind power as a function of time for a YMSC with 103 M⊙ (upper plot) and 104 M⊙ (lower plot) is shown as a solid line. The shaded region shows how the wind luminosity can statistically vary for different random sampling of the initial stellar mass function. The limits of the shaded region are calculated as the 10% and 90% percentile out of 100 possible samplings of the initial stellar mass function.

In the text

Fig. 2 Single realization of the Milky Way population of YMSCs. Top panel: face-on view of the Milky Way. The size and colour of each YMSC refer to its mass and age respectively, according to the scales indicated in the panel. Notice that the size distribution is continuous and the associations between mass and size reported in the panel are only for reference. The grey regions represents all the Galactic structures (spiral arms, Local Spur and Galactic bar) employed in the simulation of the Milky Way cluster population. The inset shows a zoomed view of the Milky Way center. The location of the Central Molecular Zone (CMZ) is also reported for reference. The star marks the position of the Solar System. Middle and bottom panels: edge on view of the Milky Way towards ROI1 (upper panel) and ROI2 (lower panel). Grey circles represent the projected sizes of the wind-blown bubbles surrounding YMSCs. The striped regions correspond to the 𝒜 ∪ ℬ masks employed in our analysis for this specific realization of the Milky Way (see Appendix B).

In the text

Fig. 3 Contribution to the diffuse γ-ray emission from a synthetic population of YMSCs in ROI1 (top row) and ROI2 (bottom row) for the three different cases of the diffusion coefficient considered: Kolmogorov (left column), Kraichnan (central column) and Bohm-like (right column). The spectrum is calculated after applying the mask 𝒜 ∪ ℬ. In each plot, the dashed line represents the median flux per energy bin after considering 100 different realisations of the Milky Way population of YMSCs, while the associated shaded region encloses the 25–75 percentile flux. The striped region instead encompasses the 25–75 percentile flux when WR stars are not included in the stellar mock populations. The violet solid lines in the right column and their associated shaded regions show the median and the 25–75 percentile diffuse γ-ray flux after the addition of emission from the CR sea.

In the text

Fig. B.1 Contribution to the diffuse γ-ray emission from a synthetic population of YMSCs in ROI1 for the three different cases of diffusion coefficient considered: Kolmogorov (upper row), Kraichnan (central row) and Bohm-like (lower row). Each column shows the flux for the three masks implemented: 𝒜 (left column), ℬ (central column) and 𝒜 ∪ ℬ (right column). In each plot, the dashed line represents the median flux per energy bin after considering 100 different realisations of the Milky Way population of YMSCs, while the associated shaded region encloses the 25-75 percentile flux. The striped region instead encompasses the 25-75 percentile flux when WR stars are not included in the stellar mock populations. The violet solid lines in the right column and their associated shaded regions show the median and the 25-75 percentile diffuse γ-ray flux after the addition of emission from the CR sea.

In the text

Fig. B.2 Example of the three masks employed to estimate the diffuse γ-ray emission in ROI1. Upper panel: mask 𝒜 used to exclude the Local Spiral Arm and a large fraction of the inner Galactic Plane where most of the γ-ray sources are expected to be located. Central panel: mask ℬ used to remove the emission from the most luminous clusters that should have been detected by LHAASO (this example mask is the one built based for the specific source realization shown in Figure 2 and assuming the Kraichnan diffusion coefficient). Lower panel: mask 𝒜 ∪ ℬ. For comparison, in all panels, the M1L mask is shown as black contours.

In the text

Fig. B.3 Distribution of the number of synthetic YMSCs detectable by LHAASO in ROI1 (upper row) and ROI2 (bottom row) for the three diffusion coefficients considered in this analysis: Kolmogorov, Kraichnan and Bohm (from left to right). We label a YMSC as detected if the significance of its emission at 100 TeV is larger than 5σ. Colour filled (striped) histograms correspond to the scenarios where WRs are (are not) included in the mock stellar populations.

In the text