The contribution of winds of star clusters to the Galactic cosmic-ray population (2024)

Giada PeronUniversité de Paris, CNRS, Astroparticule et Cosmologie, F-75013 Paris, FranceMax Planck Institute für Kernphysik, Heidelberg, Germanye-mail: giada.peron@apc.in2p3.frSabrina CasanovaInstitute of Nuclear Physics, Krakow, PolandStefano GabiciUniversité de Paris, CNRS, Astroparticule et Cosmologie, F-75013 Paris, FranceVardan BaghmanyanLehrstuhl für Astronomie, Universität Würzburg, Würzburg, GermanyFelix AharonianMax Planck Institute für Kernphysik, Heidelberg, GermanyDublin Institute for Advanced Studies, Dublin, IrelandYerevan State University, 1 Alek Manukyan St, Yerevan 0025, Armenia

Cosmic rays are energetic nuclei that permeate the entire Galactic disk. Their existence requires the presence of powerful particle accelerators. While Galactic supernova explosions may supply the required energy, there is growing evidence that they cannot explain all of the observed properties of cosmic rays, such as their maximum particle energy and isotopic composition.Among Galactic objects, winds from stellar clusters meet the energetic requirement and provide a suitable environment for particle acceleration. The recent detection of some of these objects in $\gamma$ rays confirms that they indeed harbor high-energy particles.However, as most supernovae explode inside stellar clusters, it is difficult to distinguish the contribution of winds to particle acceleration. Here we report the detection of young star clusters in the nearby Vela molecular ridge star forming region. The young age of the systems guarantees an unbiased estimate of the stellar CR luminosity free from any supernova or pulsar contamination and allows us to draw conclusions on the acceleration efficiency and the total power supplied by these objects.We demonstrate that much more than 1% of the wind mechanical power is converted into CRs and consequently conclude that a small but non-negligible fraction $\sim 1-10\%$ of the CR population is contributed by stellar clusters.

Clusters of massive stars are thought to be powerful particle accelerators. Their energy supply resides in the winds blown by member stars. For our Galaxy the total mechanical power injected by winds of stellar clusters (SCs) has been estimated to be approximately $10^{41}$ erg s^-1, which is just a factor of a few smaller than the power supplied by supernova (SN) explosions[1].Such a large amount of energy may support a sizeable fraction of the $\sim 7\times 10^{40}$ erg s^-1 luminosity of Galactic cosmic rays (GCRs) [2, 3]. Moreover, allowing particle acceleration in stellar environment could explain the anomalous excess in the isotopic ratio $X_{CR}=^{22}$Ne/²⁰Ne$\sim 0.317$ observed in GCRs (the cosmic value being $X_{S}\sim 0.0735$)[4, 5, 6].

$\gamma$-ray observations of several SCs from GeV to PeV energies[7, 8, 9] confirmed that these objects are indeed effective particle accelerators[10, 11, 12, 13].However, the fraction of wind mechanical energy that is converted into GCRs is still unconstrained.This is because the $\gamma$-ray emitting clusters detected so far reside in complex regions, where several potential particle accelerators overlap: either because some SN already exploded in the cluster (this is the case for Westerlund 1[14] and Cygnus OB2[15]) or because of the superposition with other potential CR sources along the line of sight (as in the case of Westerlund 2 [16]).This prevents an observational estimate of the contribution of SC winds to the population of GCRs.

This problem can be solved by targeting very young SCs, i.e., younger than the typical age when SN explosions begin to occur[17]. Massive stars end their life with a SN explosion, and their lifetime decreases with the stellar mass ranging from $\sim$3 Myr for the most massive stars, $\sim 120$ M_⊙, to 10 Myr for the least-massive SN progenitors, $\sim 10$ M_⊙ [18].Earlier searches in this direction resulted in upper limits on the gamma-ray flux for a handful of young star clusters [17].

Very young SCs are expected to be still embedded in the gas cocoon out of which they formed and where they spend 10%-20% of their main-sequence lifetime [19]. Even though the cocoon is optically thick to the light of the stars, the presence of embedded young SCs is revealed by the near-infrared emission that dust re-emits after being heated by the strong UV radiation of massive stars.Far infrared emission (from $\sim$ 10 $\mu$m to $\sim$ 100 $\mu$m ) acts as a tracer of early type (O/B) stars[19]. It has been demonstrated that using FIR luminosity is equivalent as using radio to identify the regions of ionized gas (Hii region) around massive stars [20]. That allowed the identification of more than 8000 Hii regions in the Galaxy using their 22-$\mu$m emission [21].Among these, several embedded young massive SCs are found in the near ($\sim$1-2kpc) Vela molecular cloud ridge (VMR) from which, making use of Fermi-LAT observations, we extracted the $\gamma$-ray spectrum and constrained the CR luminosity (see Table1). For 3 of such SCs in the VMR we report a firm detection of emission that we can relate to stellar winds due to the superposition of the high-energy and infrared emission. Their extracted energy spectra are shown in Fig. 2. All the targeted clusters are embedded in a Hii region, clearly identified both at IR, by the WISE-22 $\mu$m survey [21], and in optical by observing their characteristic H$\alpha$ emission [22].This further confirms that we are dealing with very massive stars.The typical age of embedded clusters is of the order of $\lesssim$ 1 Myr, allowing us to conclude that no SN event happened in these regions yet. The known supernova remnants in the region, Vela Y and Vela Z (also known as Vela Junior) are found in the foreground of the VMR and therefore can be easily masked to avoid any contamination to the $\gamma$-ray signal.Superposition of other possible accelerators is also excluded as we include in the sample only clusters which are not associated with any identified GeV source from the 4FGL Fermi-LAT (Large Area Telescope) catalog [23], nor with any pulsar from the ATNF (Australian Telescope National Facility) catalog [24]. We note, instead, that at the location of all of the considered SCs, unidentified Fermi-LAT sources are present and we suggest an association with the clusters.Even though no firm identification was proposed so far for these four Fermi-LAT sources, the association we put forward is consistent with tentative classifications reported in earlier works[25].

The GeV energy band is the optimal range to investigate the global energetic of CR particles, as in that part of the spectrum resides the bulk of their power.We analyzed data accumulated during 13 years by the Fermi-LAT in the energy interval between 500 MeV and 1 TeV as described in detail in the Methods section. In this energy range, the LAT point spread function ranges from 1^∘ to 0.15^∘ and its sensitivity for a 10-year long exposure is $\sim 10^{-12}$ erg cm^-2 s^-1.We detected $\gamma$-ray emission from the direction of the SCs: figures1 and 2 show the results of our analysis of the $\gamma$-ray observations carried out in the entire VMR region. Since the targeted objects are low-surface-brightness sources, we tested the stability of the results assuming three different models for the diffuse background emission: the details are reported in the Methods section where we show that the choice of the background model does not affect the detection, nor the spectral shape of the targeted sources, but only slightly (a factor 1.5 at most) the overall flux. We fitted the position and the extension of the sources. The results of the morphology tests are presented in Table 2.In a scenario where stars both heat the dust grains and accelerate CRs we expect a morphological similarity between the gamma-ray and the IR emitting regions. We therefore tested the morphology of the analyzed sources. We found a significant extension ($\sigma_{ext}>5$) only for the source RCW38, that is found to be comparable with the extension of the HII region. For the others, the extension significance is lower ($\sigma_{ext}\sim 3$), due to the limited angular resolution of Fermi-LAT, but the extension fit still converges to values similar to the extension of the Hii regions (see Table 2). In the latter cases, the extension should be regarded as an upper limit to the real extension. The centroid, instead in all cases is found to compare well with the centroid of the corresponding HII regions. The chances for random superposition are low, and are evaluated in the Methods.

The most prominent region in $\gamma$ rays corresponds to the location of RCW 38, which is detected with a statistical significance of $\sigma\sim 21$.RCW 38 is an ultra-compact Hii region that contains a young ($\sim$0.5 Myr) stellar cluster with more than 1000 stars, of which at least 30 early type stars are concentrated in a radius of $\sim$ 0.5 pc [26, 27].RCW 38 is embedded in the VelaB complex, the farthest region within the VMR. A recent precise determination of its distance has been obtained using Gaia extinction data that located RCW 38 at approximately 1600 pc from us [28]. The energy spectrum of RCW 38 is modeled as a power law, ${N_{0}(E/E_{0})^{\alpha}}$ ,with flux normalization $N_{0}=3.9\pm 0.2$ MeV^-1 cm^-2 s^-1 at $E_{0}$=1 GeV and index $\alpha_{\gamma}={-2.56\pm 0.05}$. Its luminosity above 1 GeV is $L_{\gamma}(>1\leavevmode\nobreak\ \mathrm{GeV})$=5$\times 10^{33}$ergs^-1 , assuming a distance of 1600 pc.The second brightest source in our sample coincides with RCW 36, embedded in the Vela C region, located 900 pc from us [28]. This region is slightly less luminous ( $N_{0}=3.8$ MeV^-1 cm^-2 s^-1 , $L_{\gamma}(>1\leavevmode\nobreak\ \mathrm{GeV})$= 1.5 $\times 10^{33}$ erg s^-1, detected at $\sigma\sim 16$) and slightly steeper ($\alpha_{\gamma}={-2.72\pm 0.06}$), compared to RCW 38. Similarly to RCW 38, also RCW 36 is ionized by a young ($\sim$ 1 Myr) and compact ($\sim$ 0.5 pc) stellar cluster [29].Finally, another statistically significant ($\sigma$>5) source is found to be coincident with the young star cluster RCW 32, in the Vela D region, localized at a similar distance than Vela C, $\sim 900$pc. The age of the RCW 32 is somewhat uncertain because the cluster seems to host two different populations of stars: one made of very massive stars at the core of the cluster with an age of $\lesssim$ 2 Myr, and a population of low-mass stars that formed earlier (5–13 Myr) [30]. However, this fact should not cause concern, as older low-mass stars do not end up in a SN explosion.

Other SCs are located in the VMR (see Table 2). We report a hint of $\gamma$-ray emission ($\sigma\gtrsim 5$) from RCW 35, RCW 37, RCW 40, and RCW 41. Their morphology and spectra have larger uncertainties and therefore we prefer to exclude them from the discussion, but their spectral energy distributions are reported in the Methods section. Gamma ray emission in proximity of IRS 31 is also detected at a 11 $\sigma$ level. The latter source is not included in the RCW catalog of H$\alpha$ regions [22] but it is a well known IR source[31] coinciding with the WISE region G264.124+01.926. Most likely, it belongs to the population of radio-quiet HII regions[21] which, due to their small size are believed to be optically thick to recombination lines (at both optical and radio frequencies), suggesting that also in this case the $\gamma$-ray emission could be connected to massive stars. Even if the mass in the region of IRS 31 is quite uncertain (see Table 3), this source emerges above the tested backgrounds, but with a reduced significance when changing the background and therefore we ignore it in the discussion.A significant $\gamma$-ray source with extension and position matching RCW 27 is also detected. However, the superposition with a pulsar makes it hard to associate the $\gamma$-ray emission to the embedded stars, and therefore we leave this source for a more detailed investigation. No significant emission emerged from RCW 34, this could be due to the vicinity of the source to RCW 36 and IRS 31, but also to its slightly larger distance, $\sim 2.5\pm 0.2$ kpc [32]

The detection of these three young SCs in $\gamma$-rays is unique, having no contamination from other sources and a sufficiently precise estimate of their distance and age. This allows us to estimate the fraction of the mechanical power of stellar winds that is converted into CRs.To estimate the mechanical power of winds in a cluster, $L_{*}$, we assume that the fitted size of the $\gamma$-ray source corresponds to the radius $R$ of the bubble inflated in the interstellar medium by the SC winds. The evolution of the radius of the bubble with the age of the SC, $t$, follows the well known relation[33] $R\sim(\xi L_{*}/n)^{1/5}t^{3/5}$, where $n$ is the density of the ambient gas, that can be estimated from gas column density measurements performed by making use of tracers such as CO lines and dust emission (see Methods section).More precisely, the mass determined in this way is assumed to be concentrated within spherical volume of radius equal to $R$, i.e., the measured extension of each $\gamma$-ray source. The depth of the molecular gas along the line-of-sight is certainly larger than $R$, but the dust profiles obtained by Gaia in this region puts an upper limit on the gas extent in Vela to $\sim$ 100 pc [34]. What we obtain in this way is an upper limit on the gas density, that can be converted into an upper limit for the wind mechanical power: $L_{*}<R^{5}t^{-3}n/\xi$.The parameter $\xi<1$ accounts for the fact that, due to radiative losses, only a fraction of the mechanical energy $L_{*}$ is actually used to create the bubble [35]. We constrain the value of $\xi$ to match the mechanical luminosity ($\dot{M}v^{2}$/2) estimated for RCW38, for which a mass-loss rate of $\dot{M}=2\times 10^{-4}$ M_⊙, and an average wind velocity 1000 km s^-1 are estimated [36] and find $\xi=0.07$, which is consistent with the values found in hydro-dynamical simulations of very young SCs embedded in a dense ($\sim$ 1000 cm^-3) medium [35].

The mechanical luminosity obtained in this way is compared to the observed CR luminosity, $L_{CR}$, derived from $\gamma$-ray observations, to compute the acceleration efficiency, $\eta=L_{CR}/L_{*}$.We consider a scenario where the observed $\gamma$-rays are the result of the decay of neutral pions which are in turn generated in the interactions of accelerated CR nuclei with the ambient gas. We consider leptonic processes due to CR electrons either scattering off the soft ambient photons (inverse Compton scattering) or due to bremsstrahlung less plausible (see extended discussion in the Methods section). Synchrotron losses largely dominate in these systems, given the enhanced magnetic field compared to the average in the ISM: for RCW 38 estimations suggest a value of $\sim$ 40 $\mu$G [37], that could be even larger in similar systems[38]. In such a large magnetic field it is hard to explain the observed spectral points up to 100 GeV with inverse Compton (IC) scattering or bremmstrahlung, as high energy electrons would quickly cool by emitting synchrotron photons in the radio domain[39].On the other hand, the dense gas that constitutes the VMR $n\sim$ 1000 cm^-3 provides thick target for CR hadronic interactions.

Starting from the observed $\gamma$-ray spectra (Fig.2), we derived the spectra of the parent CR protons using the naima software package [40].Then, from the derived CR spectral distribution we evaluated the CR luminosity as: $L_{CR}\approx{W_{p}}/{t_{pp}}$, where $W_{p}$ is the total CR proton energy stored within the SC and $t_{pp}$ is the proton-proton interaction timescale [41].Such an estimate of the CR luminosity is based on a calorimetric assumption: all CRs accelerated at the SC lose all their energy due to interactions with the ambient gas over a time which is shorter than both the escape time of CRs from the region and the age of the SC. While we know the age ($\sim 10^{13}$ s), the escape time can be computed from the observed $\gamma$-ray and CR luminosity as a function of $\eta$, the fraction of SC mechanical energy that goes into CRs (the derivation is reported in Methods section). For RCW 38 the escape time for 1-GeV particles results $t_{esc}=7.6\times 10^{11}\eta^{-1}({n}/{1000\leavevmode\nobreak\ \mathrm{cm^{-3}}})^{-1}\leavevmode\nobreak\ \mathrm{s}$ which is always smaller than the proton-proton interaction timescale $t_{pp}=1.6\times 10^{12}({n}/{1000\leavevmode\nobreak\ \mathrm{cm^{-3}}})^{-1}\leavevmode\nobreak\ \mathrm{s}$ for $\eta>0.005$. Similarly, the diffusion coefficient results (see Methods section) $D\sim 2\times 10^{28}\leavevmode\nobreak\ \mathrm{cm^{2}s^{-1}}\eta\leavevmode\nobreak\ (E/1\leavevmode\nobreak\ \mathrm{GeV})^{1/3}({n}/{1000\leavevmode\nobreak\ \mathrm{cm^{-3}}})$. Thus, even if a strong suppression of $D$ must take place in these systems, compared to the average Galactic value, $D_{ISM}(1\leavevmode\nobreak\ \mathrm{GeV})\sim 10^{28}\leavevmode\nobreak\ \mathrm{cm^{2}s^{-1}}$, we expect that CRs are efficiently escaping from the systems, meaning that the actual CR luminosity is much higher than the estimations obtained in the calorimetric hypothesis. Despite the large target gas density, we don’t see $\gamma$-ray emission in the surroundings of the targeted SCs. However, since the SCs themselves have a low surface brightness, on a comparable level to the diffuse emission, it is also possible that escaped particles are too diluted to cause a detectable enhancement. The only hint comes from a gas region near RCW36, approximately at $l,b\sim 266^{\circ},0.9^{\circ}$, but the uncertainties on its mass prevent us from a firm conclusion on its nature (Figure 1).

The CR luminosity calculated in the calorimetric assumption allows us to derive a strict lower limit for $L_{CR}$, as it is based on the most effective conversion of CR energy into $\gamma$-ray photons. It follows that our estimate of $\eta=L_{CR}/L_{*}$, based on a lower limit for $L_{CR}$ and on a upper limit for $L_{*}$, is a strict lower limit for the CR acceleration efficiency (see Table1).The values that we obtain are consistent with the upper limits for $\eta$ obtained from the non-detection of a number of embedded star clusters [17].

Taking the average value of ${\eta}_{min}\approx 0.5$% as representative of all SCs in the Galaxy, stellar winds would accelerate CRs at a rate at least equal to $5\times 10^{38}$ erg s^-1, corresponding to the 0.7% of the total energy supply of GCRs. Remembering that $\eta_{min}$ is a strict lower limit of the CR acceleration efficiency, we conclude that SCs provide a sizable fraction, $\epsilon_{w}\gg$ 1%, of the total power of GCRs.

Remarkably, the lower-bound value for $\epsilon_{w}=1$% that we obtained from $\gamma$-ray observations implies that an excess in the ²²Ne/²⁰Ne ratio in CRs must exist. Applying this value, in fact, we find a lower limit for the isotopic ratio $X_{CR}\sim\epsilon_{w}X_{w}+(1-\epsilon_{w})X_{S}=0.09>X_{s}$, where $X_{w}\sim 1.56$ is the estimated value for the isotopic ratio in the accelerated stellar wind material[42]. This argument can be reversed to say that, in order to reproduce the observed CR isotopic ratio $X_{CR}\sim 0.317$, the fraction of GCR particles originating from accelerated wind material has to be $\epsilon_{w}\sim(X_{CR}-X_{S})/(X_{w}-X_{S})\sim 16$%.To guarantee that this value is not overshoot, a maximum fraction $\eta_{max}\sim\epsilon_{w}\leavevmode\nobreak\ L_{tot,CR}/L_{tot,*}\sim 10\%$ of the mechanical energy of the stars may be released in the form of CR nuclei. Our results indicate that young SCs are non-negligible contributors of galactic CRs and their ability to produce high-energy particles must not be disregarded, especially for their close interconnection to the interstellar medium.

Methods

Analysis summary

We analyzed more than 13 years of Fermi-LAT PASS8 data, collected between August 8th 2009 (MET 239557417) and December 14th 2021 (MET 632287927). We selected source class events (evclass=128) reconstructed both at the front and back of the detector (evtype=3), imposed data quality DATA_QUAL==1 && LAT_CONFIG==1 and a maximum zenith angle of 90^∘, and enabled the energy dispersion evaluation (edisp=True).Due to the large extension of the VMR, we chose a region of interest 17^∘-wide centered at $(l,b)=(265.5,0)^{\circ}$. In the same region of interest we can find the Vela Junior supernova remnant (4FGL J0851.9-4620e), partially overlapping the cloud complex and the Vela X pulsar (4FGL J0835.3-4510) with its nebula (4FGL J0834.3-4542e). We included in the analysis all the sources from the 4FGL catalog[23] and the standard galactic (gll_iem_v07.fits) and extragalactic (iso_P8R3_SOURCE_V3_v1.txt) background provided by the Fermi-LAT collaboration. We then tested our detection with different galactic backgrounds based on the map of CO, Hi and dust, as explained in the next section.

We ran the analysis over the energy range 500 MeV–1 TeV. To investigate the morphology of the sources in the region, we restricted our analysis to data with energy higher than 3 GeV. In this energy range, the Fermi-LAT point spread function ($\sim 0.2^{\circ}$) is much better than at lower energies, where it could be as large as 1^∘. We deleted from the model all 4FGL sources which spatially coincided with known Hii regions, as reported in Table 2 and we re-modeled the emission in this regions. In particular, we fitted the position and the extension of these sources using the extension routine of the fermipy software package. Results are collected in Table 2: we include the best-fit position and extension and their relative uncertainties, along with the $TS_{ext}$ value which gives the significance of the extended source hypothesis over the point-like hypothesis, namely $TS_{ext}=2\log(L_{ext}/L_{ps})$, where $L_{ext,ps}$ [47] are the maximum likelihood of the extended and point-like model respectively. Although we cannot claim that all the investigated sources are extended, as they are conventionally defined extended if $TS_{ext}>16$ [48], we see that the fitted values well agrees with the extension of the corresponding Hii regions. We then used the optimized morphology, obtained at high energies, to extract the spectrum in the whole considered energy range. spectral index where the extension hypothesis is indistinguishable from the point-like hypothesis, the fitted extension can be regarded as an upper limit to the real extension. In the table, we also report the significance of the sources from our analysis, determined in the range 500 MeV-1 TeV as $\sigma=\sqrt{-2ln(L_{0}/L_{1})}$ , where $L_{1}$ and $L_{0}$ are the value of maximum likelihood obtained for a model with or without the considered source [49]. We checked that the total flux and significance did not depend on the uncertainties of the response functions. These can be accounted for by calculating the weighted likelihood [50].

Name	WISE name	4FGL source (type)	Fitted position $(l,b)$	WISE size	Size (TSext)	$\sigma$
			[∘]	[∘]	[∘]
RCW 27	G259.771+00.541	J0838.4-3952 (psr)	259.68±0.10, 0.76±0.09	0.933	0.82 111This value was computed in the energy range 500 MeV-1 TeV, because the source is not significant at higher energeies.${}^{+0.09}_{-0.08}$ ( 46.4)	9.01
RCW 32	G261.515+00.984	J0844.9-4117 (unid)	261.51±0.06, 0.95±0.08	0.440	0.26${}^{+0.06}_{-0.05}$ (10.7)	7.48
RCW 34	G264.343+01.45	–	–	0.108	–	–
RCW35	G264.681+00.272	J0853.1-4407 (unid)	264.86±0.08, -0.01±0.12	0.155	n.c.	6.73
	G264.220+00.216			0.5
RCW 36	G265.151+01.454	J0859.3-4342 (unid)	265.09±0.05, 1.36±0.04	0.224	0.27${}^{+0.07}_{-0.06}$ (10.8)	16.14
RCW 37	–	J0900.2-4608 (unid)	266.97 ±0.04, 0.01 ±0.05	–	point-like	4.03
RCW 38	G267.935-01.075	J0859.2-4729 (unid)	267.91 ±0.03,-1.03±0.03	0.155	0.21${}^{+0.04}_{-0.04}$ (26.781)	21.04
RCW 40	G269.174-01.436	–	268.57±0.06, -0.73±0.04	0.115	–	–
RCW 41	G270.310+00.851	J0917.9-4755 (unid)	270.12 ±0.13, 0.67±0.11	0.248	0.34${}^{+0.15}_{-0.09}$ (6.5)	6.96
IRS 31	G264.124+01.926	J0857.7-4256c (unid)	264.28±0.05, 1.82 ±0.04	0.075	0.14${}^{+0.06}_{-0.05}$ (4.9)	11.73
Gas core	–	J0900.5-4434c (unid)	${266.13^{\pm 0.05},0.86^{\pm 0.05}}$	–	${0.39^{+0.05}_{-0.04}(57)}$	21.22
		J0900.5-4434c (unid)

Name	Mass222d=1 kpc (dust/CO)	$N$ (dust/CO)
	[ 104 M⊙ ]	[103 cm-3]
RCW 32	0.6 /0.5	2.4/1.9
RCW 36	0.5 /0.5	2.9/2.6
RCW 38	0.4/0.5	1.9/2.1
IRS 31	0.5/0.1	7.3/ 2.5
Gas core	3.3/1.5	3.4/1.6

Systematic uncertainties

The largest source of uncertainty in the analysis of Fermi-LAT data derives from modeling of the Galactic diffuse emission. We used three different background templates to assess the validity of our results. The Galactic diffuse emission includes the pion decay emission arising from interactions of CR nuclei with the interstellar medium and the inverse Compton (IC) radiation of CR electrons with the interstellar radiation fields. Models for IC emission are usually derived with particles propagation codes, assuming a certain distribution of sources and certain boundary conditions. For what concerns pion emission instead, we tested different models including the latest released galactic background gll_iem_v07.fits. In this model the pion emission is based on maps of HI and CO and is adjusted to the data in order to match the dust emission measured by Planck and the measured $\gamma$-ray diffuse emission. The fitting procedure is done independently for different rings at different galacto-centric distances, in order to account for possible variation in the CR density towards the Galactic center. For the IC component, the standard background provided by the Fermi collaboration uses a model computed with the galprop propagation code [51]. We do the same for our customized background, choosing the output configuration ${}^{\textrm{S}}$Y^Z10^R30^T150^C2, which assumes the source distribution from Yusifov and Kucuk[52], and the boundary conditions on the scale of the galactic disk: R=30 kpc, h=10 kpc. We constructed, instead, the model for the pion emission by using spatial maps of the gas derived by CO[44] and Hi [45] line emission, and by dust opacity at 353 Hz [46]. The spectrum is modeled as a broken power law and is fitted to the data. Differently from the standard background, a single spectrum is assumed here for the entire line of sight. This is justified, since the gas in the region of Vela is all concentrated around the distance of $\sim$ 1 kpc, with the only exception of the region of Vela B which is located at $\sim$ 2 kpc. This means that the entire molecular cloud complex is concentrated in $\lesssim$ 1 kpc and therefore one should not expect a significant variation in the CR spectrum here. Spectra obtained with the different backgrounds are reported in Fig. 3. The detection of the sources and the general trend is confirmed with all the backgrounds with the exception of the region corresponding to the peak of gas, which is detected only with the standard and CO based background, but disappears when using a dust template. The values of flux calculated in the different cases vary by a factor 0.2-0.5 when using a gas template and of 0.3-0.8 when using dust. As a further check, we compared the results of our analysis to the values reported in the 4FGL-DR3 catalog [23] and saw that the resulting fluxes differ at most by 20% compared to the tabulated values. Meanwhile, the new version of the catalog (DR4) has been announced and is available in a preliminary version[53]. We checked and confirm that all sources discussed are still detected. While the significance of the sources slightly changes from one to another catalog, probably due to different modeling of the background sources, their flux is confirmed within uncertainties. Other instrumental systematic uncertainties are much smaller, and therefore can be disregarded.

Emissivity

We computed (Table 3) the mass from the observed column density of gas and dust using the conversion factors $X_{CO}=2\times 10^{20}$ cm^-2 K^-1 km^-1 s [54] and $X_{dust}=(\tau_{D}/N_{H})^{-1}=8.3\times 10^{25}$ cm^-2 [46]. Notice that in general the amount of gas and dust in each region is consistent, apart from the gas-core region and IRS 31, which are 3-times more massive in the dust template. In principle, this could suggest a component of dark gas, i.e. a region where gas tracers are saturated and therefore fail to account for all the gas. To understand if this gamma-ray excess in correspondence of the gas core could be due to a mis-calculation of the mass of the gas, we extracted the spectrum of this region and compared it to the spectrum of the large-scale region associated to the VMR. We used the customized template (based on HI and CO) and isolated from the background the region of the VMR, selected with the condition that the column density exceeded 10²² cm^-2. This region has a total mass of 2.8$\times 10^{6}$ M_⊙. The resulting emissivity, defined as the flux per hydrogen atom, is shown in Fig. 4 together with the emissivity derived from diffuse emission in the region within 8 and 10 kpc from the Galactic center [55]. The gas region of the VMR perfectly matches the emissivity of the diffuse emission, with a spectral index of ${-2.557\pm-0.003}$. For the gas core instead, the spectrum is much harder ($\alpha_{\gamma}=-2.4{\pm 0.05}$) which suggests a different origin for the emission than untraced gas. We compared also the emission of the Hii regions with the spectrum of the CR sea. The spectra, normalized to the mass of the hydrogen in that regions are shown in Fig 4. These need to be interpreted as the level of excess of these sources over the diffuse emissivity. This is because the underlying VMR is part of the background and is therefore subtracted from the flux of these regions. We see a significant variation of the spectra compared to the one of the diffuse gas. To explain the radiation in terms of a dark component, a localized enhancement of the gas density would be needed as the enhancement varies case by case from 0.5 to 3 times (see Fig.4).

Chance-correlation evaluation

Previous attempt of understanding the nature of the unidentified regions that we targeted were inconclusive. However, in a recent analysis Malyshev and collaborators[25] investigated the 4FGL unidentified sources by studying their possible association with typical source classes. They created 6 classes each including different types of sources and, using random forest (RF) and neural networks (NN) algorithms, they calculate the probability of each unidentified source of belonging to a certain class. Remarkably the sources we investigate have a large probability of belonging to class 6, that includes galactic sources, namely: binaries, pulsars, star-forming regions, supernova remnants and pulsar wind nebulae. In all the cases considered for the discussion the probability of belonging to this class exceeded 16%, in the two brightest cases the probability was larger than 40%. This argument supports our identification with stellar sources, although it should be noticed that these type of classification based on the spectral shape can not predict a new source class, and that embedded star clusters, as the ones we target could have very different spectral characteristic compared to other identified star-forming regions.

In the special case of the HII regions in Vela however, the morphological similarity between the gamma and the IR emission is an additional argument in favor of this identification. We evaluated if the superposition between Fermi and WISE sources in the region of Vela could be a chance coincidence. The 4FGL-DR3 catalog contains 2082 sources with no association. In the RoI, they match with 10 sources from the WISE catalog. We consider it a match if the Fermi source centroid falls within the radius of the HII region. We then simulated 1000 synthetic populations of 2082 sources, by performing Monte Carlo extraction over the longitude and latitude distribution of the unidentified Fermi sources. We repeated the matching procedure with the simulated sources, the distribution of the number of matches is plotted in Figure 5. We notice that with the simulations we never obtain a value as high as 10, meaning that the chances for random correlation are smaller than 0.1%. The average of the matches with the simulated populations is instead 0.4.

Leptonic scenarios

The detected sources are located in an environment with a large target gas density ($\sim 10^{3}$ cm^-3) and high magnetic field ($B\sim 40\leavevmode\nobreak\ \mu$G ), therefore nuclear interactions are naturally expected to dominate the emission, having a shorter timescale and being the electrons suppressed by cooling. One can see that in these conditions the observed spectra cannot be associated with a leptonic scenario. In the case of inverse Compton scattering on a photon field of energy $\epsilon$, for a given injection rate of electrons, $q(t)=q_{0}E^{-\alpha_{e}}$, the observed gamma-ray spectrum would be:

$$Q_{\gamma}(E_{\gamma})\sim q(E)\tau_{cool}\bigg{(}\frac{dE}{dt}\bigg{)}_{IC}\frac{dE}{dE_{\gamma}}\frac{1}{E_{\gamma}}$$

We assume as a characteristic timescale $\tau_{cool}=\frac{E}{(\frac{dE}{dt})_{IC+SYN}}$ the cooling timescale due to IC and synchrotoron interactions. The latter is proportional to the radiation and magnetic energy density respectively: $\big{(}\frac{dE}{dt}\big{)}_{IC+SYN}=A\omega_{R/B}.$ The fraction of energy transferred from electrons to $\gamma$-rays is: $E_{\gamma}=\frac{4}{3}\bigg{(}\frac{E}{mc^{2}}\bigg{)}^{2}\epsilon\propto E^{2}$, therefore:

$$Q_{\gamma}=Q_{0,\gamma}E^{-\alpha_{\gamma}}=q_{0}\leavevmode\nobreak\ E\leavevmode\nobreak\ \frac{A\omega_{R}}{A\omega_{R}+A\omega_{B}}\frac{dE}{dE_{\gamma}}\frac{1}{E_{\gamma}}\propto E^{-(\frac{\alpha_{e}}{2}+1)}_{\gamma}.$$

This implies that, in order to obtain the observed gamma-ray spectrum $\alpha_{\gamma}\gtrsim$ 2.6, an injection spectral index of $\alpha_{e}\gtrsim 3.2$ is required, in the assumption where electrons cool via synchrotron and inverse Compton radiation. The timescale for IC cooling depends on energy and on the radiation field energy density. We estimated the radiation density in the NIR, by integrating the maps provided by WISE at 22$\mu$m regions and we found values spanning from 0.5 to 50eVcm^-3. Assuming the latter as nominal value, we obtain a cooling time of $6\times 10^{7}\mathrm{\leavevmode\nobreak\ yr}(E/1\leavevmode\nobreak\ \mathrm{GeV})^{-1}$, which would exceed the age of the SCs.In the absence of cooling, the injection spectrum should be even softer to justify the observed emission. In this case, in fact, $\tau=\tau_{age}$ and one obtains $\alpha_{\gamma}=(\alpha_{e}+1)/2$ which translates to $\alpha_{e}=4.2$ for $\alpha_{\gamma}=2.6$ .

On the other hand bremsstrahlung in a dense environment is as fast as proton-proton interaction, $\sim 5\times 10^{5}\mathrm{\leavevmode\nobreak\ yr}$, for densities of 1000 cm^-3 , but in order to explain the emission via bremsstrahlung, a large electron/proton rate would be required. We can roughly estimate the ratio of electron over proton needed to explain the observed gamma-ray luminosity in terms of bremsstrahlung rather than protons. In the bremsstrahlung/hadronic dominated phase we would have $L_{\gamma}=W_{e}/\tau_{bremss}$, and in the hadronic-dominated scenario we would have $L_{\gamma}=W_{p}/\tau_{pp}$ , where

$$W_{e}=A_{e}\int_{E_{e,min}}^{E_{e,max}}E\leavevmode\nobreak\ E^{-\alpha_{e}}$$

and

$$W_{p}=A_{p}\int_{E_{p,min}}^{E_{p,max}}E\leavevmode\nobreak\ E^{-\alpha_{p}}$$

are the total energy in electron and protons respectively. By equating the two quantities and assuming that in the observed energy range $E_{p,min/max}=10\leavevmode\nobreak\ E_{e,min/max}$ and that the gamma-ray distribution follows the distribution of the parent particles: $\alpha_{e}=\alpha_{p}=\alpha_{\gamma}$ one can obtain $K_{ep}=A_{e}/A_{p}$. With our spectral index we obtain values from 0.05 to 0.25, which is much larger that what is estimated in the contest of DSA theory, namely $10^{-2}-10^{-3}$ [56]

Cosmic ray derivation

Naima [40] is a software package that allows the derivation of the parent CR distribution given a certain $\gamma$-ray SED. Here we extracted the proton spectrum from the SEDs of the SCs seen by Fermi-LAT and shown in Fig. 2, assuming that pion decay is the main radiation mechanism. By default, the software uses the parametrization of the proton-proton cross section derived by Kafexhiu and collaborators [57]. We assumed a Power-law spectrum for the protons’ total energy, $E$:

$$\frac{dN}{dE}=N_{0}\bigg{(}\frac{E}{E_{0}}\bigg{)}^{-\alpha_{CR}}$$

The fitted parameters are reported in Table 4. The total kinetic energy of protons, $W_{p}$, is also reported. The latter is calculated as $\int dT\leavevmode\nobreak\ T\leavevmode\nobreak\ \frac{dN}{dT}$, remembering that $E=T+m_{p}c^{2}$.

	$N_{0}$	$\alpha_{CR}$	$W_{p}$
	[1035 eV-1]		[1046 erg]
RCW 32	2.33${}^{+0.2}_{-0.3}$	-2.5${}^{+0.07}_{-0.05}$	1.3
RCW 36	2.8${}^{+0.4}_{-0.3}$	-2.86${}^{+0.08}_{-0.06}$	0.63
RCW 38	4.3${}^{+0.3}_{-0.3}$	-2.71${}^{+0.07}_{-0.05}$	2.71

Diffusion timescale

We compute the diffusion coefficient from the derived proton energy distribution $\frac{dN}{dE}$ derived with naima. In a situation where particle escape the system, the escape time shapes the observed distribution as:

$$\frac{dN}{dE}=q(E)\tau_{esc}=q(E)\frac{R^{2}}{D_{0}(\frac{E}{E_{0}})^{\delta}}.$$

This expression can be used to compute the CR luminosity, $L_{CR}=\eta L_{mecc}$, which we assume to be a fraction $\eta$ of the mechanical luminosity:

$$L_{CR}=\eta L_{mecc}=\int dEE\frac{dN}{dE}$$

and therefore derive $\tau_{esc}$ and $D_{0}$ as a function of $\eta$. We compute this for RCW 38, using its mechanical luminosity which is constrained both by wind-bubble theory[33] and by simulations. We find for this system $D_{0}=\eta\leavevmode\nobreak\ 1.7\times 10^{28}$ cm² s^-1 and $\tau_{esc}=\eta^{-1}\leavevmode\nobreak\ 7.6\times 10^{11}$ s.