© The Authors 2025

1. Introduction

The list of sources that cause the acceleration of Galactic cosmic rays (GCRs) is currently under revision because supernova remnants (SNRs) do not seem able to account for all the observed properties of CRs and because new classes of gamma-ray emitters were discovered that include clusters of young massive stars (e.g. Gabici et al. 2019; Amato & Casanova 2021). Star clusters (SCs) are promising CR contributors because their energetics might account for a large fraction of the power that is retained in GCRs, and because they might naturally explain the anomaly in the chemical composition of CRs with respect to the interstellar medium (ISM) (Tatischeff et al. 2021; Gabici 2024).

The gamma-ray emission originating from SCs is often ambiguous in nature: Some SCs seem to have dominant leptonic emission (e.g. Westerlund 1 Aharonian et al. 2022; Härer et al. 2023), others are convincingly hadronic accelerators (Peron et al. 2024), and the acceleration site and consequently the emission mechanism for a few objects are still debated, as in the case of the Cygnus OB2 association (Menchiari et al. 2024; Vieu et al. 2024). The nature of the emission is a fundamental issue when the contribution of this source class to the pool of GCRs is evaluated, however.

We present the analysis of gamma-ray data obtained from the M16 complex, and we argue that the emission is predominantly hadronic. M16, also known as the Eagle nebula, is one of the best-known H II regions within a few kiloparsec. It hosts and is heated by the young massive star cluster (YMSC) NGC 6611, which contains a total of 52 early-type stars, about a dozen of which are O stars. NGC 6611 has an estimated age of ∼1.3 ± 0.3 Myr (Bonatto et al. 2006). Because it is so young, no supernova is thought to have exploded in the region so far, meaning that the emission is free from contamination of these sources and their remnants. Notably, the cluster is also too young to host Wolf-Rayet (WR) stars, which appear in a later evolutionary phase. This is also confirmed by available catalogs (e.g. Rosslowe & Crowther 2015, and later updates), in which no WR is associated with NGC 6611. At the same time, because the cluster is still embedded in the gas cocoon from which it formed, a dense gas layer lies just outside it and serves as target for nuclear interactions and the production of hadronic gamma-ray emission. This is common in young stellar clusters, but the case of M16 is special in that an interacting giant molecular cloud, shaped by the wind of the stars, is unambiguously detected at the edge of the wind-blown bubble (Nishimura et al. 2021). We modeled the acceleration and propagation of particles from the termination shock to the cloud and derived clear constraints on the efficiency of particle acceleration in the region. This adds NGC 6611 to the handful of young star clusters for which an estimate of the acceleration efficiency has been possible.

The Letter is structured as follows: We describe in Sect. 2 the gamma-ray and gas measurements toward the region. We outline in Sect. 3 the acceleration and propagation model we applied to the system, and we derive constraints on the acceleration efficiency. We finally summarize our findings in Sect. 4 and discuss their implications.

2. Observations

We analyzed Fermi-LAT data collected in more than 16.5 years (from August 2008 to March 2025) from the region of M16. More precisely from a region centered at (l₀, b₀) = (17, 1)°. We applied standard data-quality cuts and standard models for background sources, as detailed in the Appendix A. In the region of M16, the unidentified source 4FGL J1818.0-1339c is listed in the 4FGL catalog. We removed it from the source model and remodeled the emission. The residual test-statistics (TS) map is shown in Fig. 1. We tested the position and extension of this source considering only events with energy > 1 GeV, and we obtained a best-fit position of (l_*, b_*) = (17.10 ± 0.05, 0.93 ± 0.05)° and a best-fit extension of 0.27_− 0.05°°^+ 0.07°, with a TS_ext of 17.8 and a total detection significance of the source of 7.5σ. We modeled the emission as a uniform disk with the reported extension and coordinates. This morphology was then used to extract the spectral points in the entire energy range (60 MeV–870 GeV). The resulting spectral points are shown in Fig. 3. They were derived with a standard procedure, namely by fitting the normalization of the flux in each considered energy bin¹.

Fig. 1. Left: Test statistics map in the region around M16 derived from Fermi-LAT observations. The position of the O stars in the cluster NGC 6611 is indicated as blue stars, and its termination and forward shock radii are indicated a solid and dotted blue circles. The dashed circle indicates the limit of the H II region. The light blue, cyan, and pink contours indicate CO gas column densities of 1.5 × 1022, 1.7 × 1022, and 2.3 × 1022 cm−2, which were obtained in the velocity range 15–30 km s−1. The magenta circles indicate the position of Fermi sources from the 4FGL catalog Ballet (2023) and their relative position uncertainties, evaluated as their 68% confidence radius. Center: Molecular gas distribution in the region as traced by CO line emission, integrated over the line of sight in the velocity range 15 km/s, 30 km/s. The black contours refer to Fermi-LAT test statistics and represent 3, 4, and 5σ levels. Right: Velocity-longitude distribution of CO. The map is integrated over the whole latitude range. The dot indicates the position of NGC 6611.

Fig. 3. Spectral energy distribution (SED) of the source detected by Fermi-LAT in correspondence of NGC 6611/M16. The error bars refer to the statistical 1σ uncertainties. The solid and dashed green lines represent the emission model in the case of Kraichnan turbulence, and the lines differ in the assumed wind luminosity, Lw, as indicated in the figure legend.

To derive the gas distribution, we analyzed the CO(J = 2  →  1) line emission map provided by Dame et al. (2000), which traces the dense molecular gas. The results are shown in Figure 1 (central and right panels) in the longitude-velocity and longitude-latitude planes. At the selected coordinates, an isolated gas structure is located within the radial velocity range 15–30 km s⁻¹. This well-known structure corresponds to a molecular cloud associated with the cluster that appears to be compressed by the interaction with it (Nishimura et al. 2021). The gamma-ray emission (shown in the left panel of Fig. 1) is spatially coincident with this gas structure, which suggests that the emission originates from the interaction of high-energy nuclei with this cloud. From the CO map, we estimated the gas density associated with this system, which was used to compute the gamma-ray emissivity and to estimate the size of the bubble created by the stellar wind, as discussed in Sect. 3. We assumed a conversion factor from CO to H₂ of X_CO = 2 × 10²⁰ cm⁻² K km s⁻¹ (Bolatto et al. 2013) and used the all-sky map of H I line emission by Ben Bekhti et al. (2016) for the atomic component, with a conversion factor of X_HI = 1.8 × 10¹⁸ cm⁻² K km s⁻¹. The mass of the ionized component was neglected because it can be demonstrated that this is much more diluted than the neutral component. We used the radius of the H II region as traced by WISE (Anderson et al. 2014) as the extension of the ionized bubble. A dense layer of swept-up gas was assumed to occupy a shell in the immediate surroundings of the bubble. We assumed that the radius of the shell is ∼5% R_HII, which is consistent with the values obtained in simulations (Gupta et al. 2016). We therefore measured the mass contained within 1.05 R_HII. At the distance measured by Gaia of 1646 pc (Cantat-Gaudin et al. 2020), the resulting mass of the swept up material is 23.6 kM_⊙.

3. Derivation of the acceleration efficiency

No SNRs are expected to be present and accelerate particles in the system because the cluster is so young. As an alternative accelerator, we considered the collective wind termination shock from the NGC 6611 cluster and applied the steady-state model developed by Morlino et al. (2021) to describe the emission and infer the acceleration efficiency. We considered two different values for the wind luminosity, L_w, 1 = 6.3 × 10³⁶ erg s⁻¹ and L_w, 2 = 3.8 × 10³⁷ erg s⁻¹; the two values reflect the uncertainty on the stellar mass of the cluster, as argued in Appendix B. The wind luminosity powers the acceleration and injects turbulence into the system and thus contributes to an efficient confinement of the accelerated particles. To compute the diffusion coefficient, we assumed that a fraction η_B = 5% of the wind power is converted into magnetic turbulence, with a coherence length of l_c = 1 pc (roughly equal to the SC core radius), and we considered the three turbulence spectra that are most commonly adopted in astrophysics: the Bohm, Kraichnan, and Kolmogorov spectra. The magnetic field that results from this hypothesis is 4.7 μG at the wind-termination shock, which corresponds to a value of ∼13 μG downstream. The latter is largely below the upper limit found by Mackey & Lim (2011) based on magnetohydrodynamical (MHD) simulations of the famous pillar structures (the pillar of creation) observed in M16. We then modeled the geometry of the system as sketched in Fig. 2. We considered that most of the mass that we measure from the gas tracers is swept up to form a dense shell, whereas the mass inside the bubble is mostly contributed by the evaporation of the shell (Castor et al. 1975). The innermost parts of the bubble (r < R_TS) are populated by the stellar wind material, computed from Ṁ, which was in turn derived from the cluster mass as in Menchiari et al. (2025). In a 1D model, the shell is very thin because of radiative cooling. HD instabilities fragment the shell, however, which increases its effective thickness. According to El-Badry et al. (2019), the typical length scale to which the fragmentation instability grows over a time interval Δt is about $\lambda \sim 2\pi \mathrm{\Delta }t{v}_{\mathrm{rel}}{r}_{\rho }^{-1/2}$, where v_rel is the shear speed between the shell and the bubble, and $r_{\rho }=(\frac{{\rho }_{\mathrm{high}}}{{\rho }_{\mathrm{low}}})$ is their density ratio. When Δt is fixed to the age of the cluster and reasonable estimation for v_rel ∼ 10 km s⁻¹ is considered and r_ρ ∼ 10³, as inferred from the observed densities, the relevant instability scale is about 2 pc. In order to account for this, we assumed the mass of the shell to be distributed according to exp[−((R_fs − r)/0.5 pc)²] θ(R_fs − r) θ(R_fs − r − 2 pc), where θ is the Heaviside step function, and the last term takes the above estimate of the penetration length into account. We instead neglected the mass of larger but rarer structures such as the pillars; their estimated mass is ∼200 M_⊙ (McLeod et al. 2015), which amounts to only ∼1% of the swept-up mass of the shell.

Fig. 2. Distribution of the gas density as a function of the radial distance from the center of the cluster (black curve) and the radial distribution of accelerated particles at different energies (solid, 10 GeV; dotted, 100 GeV), considering propagation in Bohm (blue), Kraichnan (green), or Kolmogorov (red) turbulence. The highest value for the wind luminosity, Lw, 2, was assumed in the calculations.

Finally, we constrained the size of the wind-blown bubble, namely the location of the forward shock, R_fs, based on the observed size of the H II region. The latter is 0.312° (Anderson et al. 2014), yielding R_fs = 1.05 R_HII = 9.4 pc at the given distance. By comparing the observed radius with expectations from Weaver’s equation, which we modified to take losses into account,

$$\begin{array}{c}\hfill R_{\mathrm{fs}}=18.14(\frac{{\eta }_m{L}_w}{{10}^{37}{\mathrm{erg}\mathrm{s}}^{-1}}{)}^{\frac{1}{5}}(\frac{n_0}{197{\mathrm{cm}}^{-3}}{)}^{-\frac{1}{5}}(\frac{t}{1.3\mathrm{Myr}}{)}^{\frac{3}{5}}\mathrm{pc},\end{array}$$(1)

we derived the mechanical efficiency, η_m, which accounts for the fraction of the wind-luminosity that inflates the bubble (see, e.g. Vieu et al. 2022). The density n₀ was computed from the gas map described in Sect. 2 over the volume of the bubble, giving n₀ = 197 cm⁻³. The resulting value of η_m is then ∼8% for L_w, 1 and ∼1.3% for L_w, 2, which agrees well with the values found in the simulations of Yadav et al. (2017) for young systems in a dense medium. The value of η_m also defines the size of the termination shock, resulting in R_TS = 2.17 pc and R_TS = 5.6 pc for the two luminosity values (see Eq. (C.1)).

After we fixed the wind luminosity, the propagation regime, and the density profile, the gamma-ray emission only depended on the CR distribution, which we assumed to follow the solution of Morlino et al. (2021), namely a power law in momentum with a cutoff that depends on the turbulence spectrum ∝p^−α_CRe^−Γ(p) (see the extended functional form in Appendix C), and on the acceleration efficiency, η_CR, which is defined as the fraction of L_w that goes into accelerated particles. We only considered the hadronic contribution through pion decay because, as shown in Appendix E and discussed in the following section, a possible leptonic contribution is thought to be under-dominant. For each model of diffusive transport we considered (Kolmogorov, Kraichnan, and Bohm), we derived the best combination of α_CR and η_CR that fits the gamma-ray spectral energy distribution measured with Fermi-LAT. The goodness of the fit was evaluated with a χ²-test, which was computed over a grid of possible values of α_CR and η_CR. More details of this procedure are provided in Appendix D. An example of the best-fit SED is shown in Fig. 3, together with the measured spectral points. The calculation was made for a Kraichnan-type diffusion. The plot clearly shows that the assumption we made for the luminosity only affects the spectral shape at the highest energies, where a significant detection is unfortunately lacking. The results for the different cases are reported in Table 1, where we report the best-fit values for the acceleration efficiency and spectral slope, together with the corresponding reduced χ².

4. Discussion and conclusions

We have detected significant gamma-ray emission from the M16 region in the energy range from 60 MeV to 20 GeV. The emission is compatible with an extended source with a radius of ∼0.3°, and it is found in the vicinity of the star cluster NGC 6611, in coincidence with the location of a gas cloud that interacts with the wind-blown bubble of the system. This adds this cluster to the list of gamma-ray emitting young (< 3 Myr) star clusters whose emission must be powered by stellar winds because no SN has occurred so far. The spatial coincidence with dense target gas suggests a hadronic origin of the emission, which we interpreted as due to particles that were accelerated at the wind-termination shock of the cluster and then propagated to the cloud, where the interaction takes place. While formally, the coincidence with gas over-densities is expected also for bremsstrahlung emission, an explanation of the observations in terms of leptonic interactions requires in this case that electrons are accelerated more effectively than ions, and with a rather peculiar spectrum. Normally, in diffusive shock acceleration, the number of accelerated ions is much higher (100–1000 times) than the number of accelerated electrons, and hadronic emission therefore always dominates bremsstrahlung².

Even if the shock were to accelerate leptons alone, however, the emission can only be explained by lowering the highest electron energies to ≲2 TeV. In the case of Bohm or Kraichnan turbulence, this could in principle be self-consistently justified by invoking very large magnetic fields (on the order of hundreds of micro-Gauss) that would decrease the maximum energy via radiation losses. However, not only does this field strength exceed the upper limit found by Mackey & Lim (2011), it even requires the nonphysical condition η_B > 1. When Kolmogorov-type turbulence is assumed, the conditions on E_max are relaxed, but this requires a spectrum with α < 4, which we consider unrealistic. For these reasons, we argue that the emission must be predominantly hadronic, with a subdominant contribution of leptons that is discussed in Appendix E.

These considerations allowed us to derive the proton acceleration efficiency of the system, η_CR. We presented the results using upper and lower bounds for the wind kinetic luminosity derived from different estimates of the cluster mass from the literature, but we considered the upper bound (L_w, 2) closer to the real value because the lower value (L_w, 1) derives from an underestimation of the stellar mass. We also note, however, that the gamma-ray spectrum can also be explained using the lower bound and still reasonable values for the other parameters (see Table 1). The fit is mainly sensitive to the type of assumed turbulence: We note that for Kolmogorov-type diffusion, a very hard spectrum of accelerated particles is required to fit the data, with α_CR < 4, and an acceleration efficiency between 10% and 40% is needed, which is difficult to reconcile with theory of diffusive shock acceleration. The Bohm assumption instead requires a rather steep acceleration spectrum (α_CR ∼ 4.3), which agrees well with diffusive shock acceleration, because any harder spectral shape would violate the upper limits from observations (see also Appendix D). While we cannot exclude this injection spectrum, we remark that Bohm is generally considered an extreme case for turbulence, which provides an extremely efficient confinement. In a system with Bohm-like turbulence, an acceleration efficiency of 1.2% would be sufficient to explain the observations. For the considerations above, the latter can be regarded as a stringent lower limit to the acceleration efficiency in the system. Finally, the Kraichnan case falls somewhat in between, with a slope α_CR ∼ 4 and an efficiency of 3.6%. Our conclusion is that the acceleration efficiency can be constrained between 1.2% and 3.6%. These results depend only weakly on our choice of parameters, as discussed in Appendix C.1, where we show that the acceleration efficiency can change by a factor of 2 at most. Interestingly, the values of η_CR derived here are similar to the values derived by Peron et al. (2024) for the young embedded clusters RCW 32, RCW 36, and RCW 38 in the Vela molecular ridge, but the estimates in these cases were obtained under much rougher assumptions on the particle confinement inside the wind bubble. Following the same approach as we proposed there, we estimated the contribution of stellar winds to the production of CRs. We assumed that all stellar wind systems during the main sequence accelerate particles, on average, with the same efficiency as derived for NGC 6611. To do this, we considered that the total wind luminosity provided by OB stars throughout the Galaxy is estimated to be 1 × 10⁴¹ erg s⁻¹ (Seo et al. 2018) and that the total luminosity of cosmic rays is 7 × 10⁴⁰ erg s⁻¹ (Strong et al. 2010). With an efficiency of 3.6% (1.2%), the fraction of CR luminosity contributed by the winds of SC, ϵ_w = η_CRL_tot, SC/L_tot, CR, is 5.1% (1.7%), which agrees well with the estimate ϵ_w ∼ 6% derived by Tatischeff et al. (2021) from considerations of the chemical composition of CRs. In summary, this study of NGC 6611 adds evidence to the idea that SC winds provide a small but still important fraction of Galactic CRs.

Acknowledgments

The authors are thankful to the participant of the TOSCA-24 workshop for enlightening discussion. GP is supported by the Inaf Astrophysical Fellowship initiative. This work was partially funded by the European Union – NextGenerationEU RRF M4C2 1.1 under grant PRIN-MUR 2022TJW4EJ and by INAF under Theory Grant 2024 Star Clusters As Cosmic Ray Factories II. SM acknowledges financial support from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033.

References

Appendix A: Fermi-LAT analysis details

For the Fermi-LAT analysis, we considered data in a region of interest (RoI) of 8° around (l₀, b₀) = (17, 1)°. We performed a binned analysis, considering spatial steps of 0.1° × 0.1°, and 8 energy bins per decade, logarithmically spaced in the range between 60 MeV and 870 GeV. We applied the recommended data-quality cuts for source analysis (DATA_QUAL> 0 && LAT_CONFIG= = 1, with maximum zenith angle set to 90°) and selected only events of type evclass = 128 converted both at the front and at the back of the detector (evtype = 3). As a starting model for the emission, we used the latest released catalog of sources (4FGL-DR4,Ballet (2023), Abdollahi et al. (2022)), together with the Galactic and extra-galactic emission templates provided by the Fermi collaboration (namely gll_iem_v07.fits and iso_P8R3_SOURCE_V3_v1.txt). After a first optimization of the parameters, we included two additional point-like sources to account for two hot spots in the residual test-statistic (TS) maps that exceeded a TS value of 9. The additional sources, located at (l₁, b₁) = (17.933, −2.75)° and (l₂, b₂) = (17.949, −2.02)° are well separated (> 1°) from the source of interest, and hence are of no concern in terms of systematic uncertainties. We remove then the unidentified source 4FGL J1818.0-1339c from the model; the resulting TS map can be seen in Fig. 1 and A.1. We performed an extension test, which positively detected an extension of 0.27${}_{-0.05}^{+0.07}$ at a significance level of TS_ext = 17.5. We found no appreciable difference in modeling the source as a Gaussian or as a disk, since the log-likelihood of the overall fit in the two cases differs by less than 0.1σ. We further tested a possible curvature in the spectral energy distribution, by comparing the likelihood of a power-law model with the likelihood of a log-parabola, and calculated TS_curv = 2(log L_LP − log L_PL). The resulting TS_curv = 8.37 is well above the threshold (TS_curv > 4) set in the 4FGL catalog Abdollahi et al. (2022) for preferring a curved spectrum over a power-law spectrum. This means that in the 4FGL catalog our source would be flagged as a log-parabola, even if the curvature significance is below 5, as ${\sigma }_{\mathrm{curv}}=\sqrt{T{S}_{\mathrm{curv}}}\lesssim 3$. This relatively low threshold is adopted in the catalog because the power-law spectrum tends to overestimate the flux at the edge of the sensitivity range, and this is particularly important when catalog sources are used as background. Given that our aim is to characterize the flux using a more sophisticated modeling, based on acceleration and interaction of protons, we find the power-law description satisfactory for the overall fit. The resulting best fit parameters for the power-law are: α = 2.14 ± 0.06 and $F_0(E_0=1\mathrm{GeV})=(1.7\pm 0.2)\times {10}^{-12}{\mathrm{MeV}}^{-1}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$. We also report, for completeness, the best-fit parameters obtained with a log-parabola function, namely:

$$\begin{array}{c}\hfill F(E)=F_0(\frac{E}{E_b}{)}^{(-\alpha +\beta ln(\frac{E}{E_b}))}\end{array}$$

with $F_0(E_0=1\mathrm{GeV})=(1.8\pm 0.3)\times {10}^{-12}{\mathrm{MeV}}^{-1}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$,α = 1.77 ± 0.13 and β = 0.151049 ± 0.000005, while E_b is fixed to 1 GeV.

As a last step, we extracted the spectral energy distribution, by fitting in every energy bin the normalization of a power-law spectrum with spectral index 2. We checked that this choice does not affect the spectral points, by using the locally measured index, namely, the spectral index derived from the overall fit, for each pixel. We further repeated the procedure with both the power-law and the log-parabola models and by varying the size of the energy bins, finding no appreciable difference in the results. We evaluated systematic uncertainties using the bracketing method suggested by the Fermi collaboration³ but they resulted smaller than the statistical errors.

Fig. A.1. Test statistic map of the region after modeling all sources except for 4FGL J1818.0-1339c. All 4FGL in the map are marked as magenta circles. The extension of the circle represents their localization uncertainty.

Appendix B: Mass and luminosity of the cluster

The wind luminosity of NGC 6611 has been estimated by Celli et al. (2024) to be L_W,1 = 6.3 × 10³⁶ erg s⁻¹, which corresponds to a total stellar mass M_∗,1 = 1697 M_⊙. However, the latter is considered by the same authors to be a strict lower limit, because of the limitation of their method in resolving the stellar population in young clusters (we refer to Sec. 4 and 5 in Celli et al. (2024) for a more detailed explanation). A similar value for the cluster mass is derived independently by Bonatto et al. (2006), who also argue that the value should be considered as a lower limit, because it results from the distribution of stars less massive than 5 M_⊙. A later study (Selim et al. 2016) derived a mass of 2260 M_⊙ based on observations of stars up to 35 M_⊙, while estimates by Stoop et al. (2023) suggest a mass range between 5 and 10 kM_⊙, when extrapolating the mass distribution up to stars 100 M_⊙. We therefore considered a second value for the wind luminosity, L_W,2 = 3.8 × 10³⁷ erg s⁻¹, which we derived following the method presented by Menchiari et al. (2025): we simulated the stellar population of the cluster, assuming the initial mass function reported by Kroupa (2001) and fixing the total mass of the cluster to M_*, 2 = 3 M_*, 1. The choice of M_*, 2 = 3 M_*, 1 is motivated by the estimation of the systematic uncertainty that Celli et al. (2024) argue to affect their sample. Moreover, the estimated value of M_*, 2 is consistent with the mass range derived by Stoop et al. (2023).

Appendix C: Modeling the acceleration and propagation within the wind-blown bubble

The model for particle acceleration and propagation adopted here is the one proposed by Morlino et al. (2021) and applied by Menchiari et al. (2024) to the Cygnus region. The model assumes that particles are accelerated at the wind-termination shock and then propagate by advection and diffusion in the low-density bubble, made turbulent by the development of MHD instabilities. We refer to the literature cited above for a full description of the solution. Here, we just recall the basic ingredients of the model. We further recall that the model is a steady-state model; this approximation is good enough in SCs, as the wind luminosity stays basically constant throughout the entire main-sequence phase (see e.g. Menchiari et al. 2025; Vieu et al. 2022), while the size of the bubble and of the termination shock grows slowly. The termination shock radial position is calculated as:

$$\begin{array}{cc}\hfill R_{\mathrm{ts}}=26{\eta }_m^{-\frac{1}{5}}& (\frac{\dot{M}}{{10}^{-4}{\mathrm{M}}_{\odot }y{r}^{-1}}{)}^{\frac{3}{10}}(\frac{V_w}{2000km{s}^{-1}}{)}^{\frac{1}{10}}\times \hfill \\ \hfill & \times (\frac{n_0}{10{\mathrm{cm}}^{-3}}{)}^{-\frac{3}{10}}(\frac{t}{10\mathrm{Myr}}{)}^{\frac{2}{5}}\mathrm{pc}\hfill \end{array}$$(C.1)

,while the expression for the forward shock is reported in Eq. 1.

The distribution function is obtained by solving the steady state transport equation in radial symmetry using two boundary conditions: net zero flux at the center of the bubble, matching of the Galactic CR distribution f_gal (tuned on near-Earth observations from Aguilar et al. 2015) at large distance. We can distinguish three regions, labeled from 1 to 3: the cold wind, the shocked wind (i.e. the hot bubble) and the region external to the bubble (see Fig 2). In the cold wind the CR distribution can be approximated as

$$\begin{array}{c}\hfill f_1(r,p)\simeq f_{\mathrm{ts}}(p){\mathrm{e}}^{-v_w(R_{\mathrm{ts}}-r)/D_1},\end{array}$$(C.2)

inside the bubble as:

$$\begin{array}{cc}\hfill f_2(r,p)& =f_{\mathrm{ts}}(p){\mathrm{e}}^{\alpha (r)}\frac{1+\beta ({\mathrm{e}}^{\alpha (R_b)}{\mathrm{e}}^{-\alpha (r)}-1)}{1+\beta ({\mathrm{e}}^{\alpha (R_b)}-1)}\hfill \\ \hfill & +f_{\mathrm{gal}}(p)\frac{\beta ({\mathrm{e}}^{\alpha (r)}-1)}{1+\beta ({\mathrm{e}}^{\alpha (R_b)}-1)},\hfill \end{array}$$(C.3)

while outside as:

$$\begin{array}{c}\hfill f_3(r,p)=f_b(p)\frac{R_b}{r}+f_{\mathrm{Gal}}(p)(1-\frac{R_b}{r})\end{array}$$(C.4)

where f_b(p) = f₂(R_b, p) is the particle distribution at the bubble boundary. The functions α and β are:

$$\begin{array}{cc}\hfill \alpha (r,p)& =\frac{u_2{R}_{\mathit{ts}}}{D_2(p)}(1-\frac{R_{\mathit{ts}}}{r})\hfill \end{array}$$(C.5)

$$\begin{array}{cc}\hfill \beta (p)& =\frac{D_3(p)R_b}{u_2{R}_{\mathit{ts}}^2},\hfill \end{array}$$(C.6)

where the diffusion coefficients are assumed to be spatially constant within each of the regions: D₁ in the cold wind, D₂ in the shocked wind region and D₃ in the ISM. D₃ is assumed equal to the average Galactic diffusion coefficient, while D₁ and D₂ are determined by the level of magnetic turbulence inside the bubble, which is largely unknown. For this reason, we explore different types of turbulence spectra (Kolmogorov, Kraichnan and Bohm) normalized in such a way that a fraction η_B = 5% of the wind kinetic energy is converted into magnetic energy, namely

$$\begin{array}{c}\hfill 4\pi r^2{v}_w\frac{\delta B_w^2}{4\pi }={\eta }_B\frac{1}{2}\dot{M}{v}_w^2.\end{array}$$(C.7)

For Kolmogorov and Kraichnan we assume an injection scale fixed to 1 pc (roughly corresponding to the size of the stellar cluster core) while the Bohm case implies that the injection also occurs on smaller scales with similar power.

The solution at the termination shock is:

$$\begin{array}{c}\hfill f_{\mathrm{ts}}(p)=\frac{{\eta }_{\mathrm{CR}}}{2}\frac{\sigma n_1{v}_w^2}{4\pi \mathrm{\Lambda }{(m_{p}c)}^3{c}^2}{\left(\frac{p}{m_{p}c}\right)}^{-{\alpha }_{\mathrm{CR}}}{e}^{-\mathrm{\Gamma }(p/p_{max})}.\end{array}$$(C.8)

In the above equation n₁ and v_w are the plasma density and speed upstream of the termination shock, while σ is the compression ratio at the shock, m_p is the proton mass and c the speed of light and p_max = E_max/c is the maximum momentum of accelerated particles. Notice that f_ts is normalized by requiring that the cosmic-ray luminosity is a fraction of the wind luminosity, namely η_CRL_w = L_CR ≡ 4πR_ts²v_w/σ∫Ef_ts(E)dE. This determines the normalization constant $\mathrm{\Lambda }={\displaystyle {\int }_{x_{\mathrm{inj}}}^{\infty }}{x}^{2-{\alpha }_{\mathrm{CR}}}[{(x^2+1)}^2-1]{\mathrm{e}}^{-\mathrm{\Gamma }(x)}dx$. Finally, the exponential function e^−Γ depends on the ratio p/p_max and on the diffusion coefficient chosen. E_max is derived from the assumptions we made on the magnetic field intensity and on the diffusion coefficient, equating the diffusion length to R_TS

$$\begin{array}{c}\hfill D_1(E_{max})/V_w=R_{\mathrm{TS}},\end{array}$$(C.9)

where we recall that D₁ is the diffusion coefficient in the wind region. The result is: E_max= 18.8 TeV, 115.3 TeV, and 708.5 TeV in the Kolmogorov, Kraichnan and Bohm regime, respectively. These energies correspond to gamma-rays out of the energy range probed by Fermi. Therefore, E_max cannot be constrained using Fermi-LAT data alone. As a consequence, our model depends only on two parameters, the acceleration efficiency (η_CR) and the spectral slope at injection (α_CR), which we determine by fitting the gamma-ray emission.

C.1. Sensitivity to model parameters

When fitting our model to the data, we fix some parameters to benchmark values. Here we estimate how much their variation from reference values can affect the calculation of the acceleration efficiency. We considered in particular the fraction of wind luminosity that is converted into magnetic energy, η_B, and the coherence length of the turbulence, l_c. Assessing their value is only possible through numerical simulations of the wind-bubble system, which is much beyond the aim of this paper. Nevertheless, we can use some educated guess on their range: for η_B we obtain a minimum value by imposing that the magnetic field in the bubble is equal to the average magnetic field in the ISM (3 μG), that corresponds to η_B = 0.2%, while the maximum value can be estimated as η_B = 10% based on the results of numerical simulations (see e.g. Badmaev et al. 2022); for l_c, we consider as minimum value the average distance between the stars in the cluster core (considering the values reported by Hunt & Reffert 2023), namely 0.5 pc, and the maximum value comparable with the extension of the termination shock radius, namely ∼6 pc. We then compute, for different combinations of values in these ranges, the quantity Φ(E, r)≡f_CR(E, r)/f_ts(E), namely the CR distribution function normalized to its value at the termination shock, and compare it with Φ⁰ ≡ Φ(η_B = 0.05, l_c = 1 pc), namely the same quantity computed for reference values of the parameters. Φ depends on energy and on the distance from the cluster core, but, since most of the gamma-ray emission derives from the shell region, we limit our investigation to this location, namely at r = 9.6 pc. The resulting ratios of Φ/Φ⁰ at the shell location are shown in Fig. C.1 as a function of η_B and l_c for two different particle energies, and for different turbulence spectra. For the Kraichnan case the plots show that there is no combination of parameters that cause a variation larger than a factor of two simultaneously for both energy bins. This means that the acceleration efficiency can vary at most by a factor of two, by changing η_B and l_c. For Bohm, instead, the allowed variation is much smaller, of the order of 2%.

Fig. C.1. Variation of Φ/Φ0 as a function of ηB and lC (see the text for a complete description). The white contour visible in the upper right panel corresponds to Φ/Φ0 = 0.5. In each panel the circle marks the reference parameters (lc = 1pc, ηB = 0.05)

Appendix D: Model fitting

To derive the unconstrained parameters, we fit the model described in Appendix C to the derived spectral data, with the caveat of imposing f_Gal(p) = 0, to account that this component is already subtracted in the analysis. The two free parameters of our model, namely the acceleration efficiency η_CR and the injection spectral slope α_CR, are obtained by a χ² fit of the data. Having 5 data points and 2 parameters, our system has 3 degrees of freedom. We perform a grid-search to find the parameters that minimize the reduced χ². For the spectral index we considered the intervals ${\alpha }_{\mathit{CR}}^{\mathrm{Kol}}\in [2(2),4.1(4.3)]$, ${\alpha }_{\mathit{CR}}^{\mathrm{Kra}}\in [3(3.2),4.9(4.9)]$, and ${\alpha }_{\mathit{CR}}^{\mathrm{Bohm}}\in [3.5(3.8),4.9(4.8)]$, while for the efficiency we used ${log}_{10}{\eta }_{\mathit{CR}}^{\mathrm{Kol}}\in [-1.5(-2.5),0(-0.3)]$, ${log}_{10}{\eta }_{\mathit{CR}}^{\mathrm{Kra}}\in [-2(-2.8),0(-0.4)]$, and ${log}_{10}{\eta }_{\mathit{CR}}^{\mathrm{Bohm}}\in [-2.5(-3),-0.4(-0.9)]$, where the pairs correspond to L_W, 1 (L_W, 2). We used equally spaced grids with bin sizes of δα_CR = 0.2 and δlog₁₀η_CR = 0.2, for the injection spectral slope and the efficiency respectively. In addition to the χ² evaluation, we further excluded the region of the parameter space that violates the upper limits derived from the Fermi-LAT analysis. The latter appears as a dashed area in Fig. D.1, where the distribution of χ² is shown for all six cases considered (two values of the wind luminosity times three models for the diffusion coefficient). If the best fit set of values fell in the region of parameter space prohibited by the upper limits, we would then choose a second minimum in the allowed area. The corresponding SEDs are shown in D.2.

Fig. D.1. Parameter space explored with the fitting procedure. The darker zones are values with a smaller χ2, while the shaded regions represent those parameters excluded by the gamma-ray upper limits. The best fit value is indicated as circle, unless it is found in the forbidden parameter region, in which case the new best fit value is indicated as a cross.

Fig. D.2. Three different fit results to the spectral energy distribution (SED) of the source detected by Fermi-LAT in correspondence of NGC 6611/M16. The solid and dashed curve differ in the assumed wind luminosity, Lw, as indicated in the figure legend. The three different panels refer to the three turbulence spectra assumed: Kolmogorov (red), Kraichnan (green), and Bohm (blue).

Appendix E: The leptonic contribution

In this appendix we provide an estimate of a possible leptonic contribution to the gamma-ray emission as due to non-thermal bremmstrahlung and inverse Compton (IC). The former, like pion decay, depends only on the gas density, while IC depends on the photon radiation field density. We considered both CMB radiation (with a radiation energy density of 0.261 eV cm⁻³) and an infrared radiation field due to the local environment. We inferred the latter from observations by integrating the best-fit SED derived from the shell region of M16 (Flagey et al. 2011). We found an IR radiation density of ∼9 eV cm⁻³ for the warm dust component at 70 K. We also estimate the radiation field contributed by the massive O stars of the cluster, by summing up the bolometric luminosity of each star, L_i, reported by Guarcello et al. (2012). Considering an average temperature of 3 × 10⁴ K, we then compute the radiation density as ∑_iL_i(4πr²c)⁻¹ ∼ 50 eV cm⁻³(r/5 pc)⁻², where r is the distance from the center of the cluster. We evaluate the radiation density at R_TS and the resulting IC contribution. This choice corresponds to an upper limit on the radiation density, that rapidly decreases outwards and, at the same time, to a lower limit on the stellar temperature, both tending to maximize the IC intensity. However, even with these values of temperature and radiation density, this IC component turns out to be negligible compared to IC on infrared light and bremsstrahlung.

As for the electron population, we considered the same injection spectrum as for protons (see eq. C.8) but with a different normalization and a different maximum energy, i.e.

$$\begin{array}{c}\hfill f_{\mathrm{ts},\mathrm{e}}(p)=K_{\mathrm{ep}}{f}_{\mathrm{ts},\mathrm{p}}(p){\mathrm{e}}^{-p/p_{max,e}}.\end{array}$$(E.1)

The normalization is scaled to that of protons by a factor K_ep, left as a free parameter. Notice that typical values of K_ep for the case of SNR shocks range between ∼10⁻⁴ and few times 10⁻³.

The electron maximum energy is determined by comparing the acceleration time with the energy-loss time due to bremsstrahlung, IC and synchrotron emission. For the latter, the assumed magnetic field strength inside the bubble is ∼13 μG, which corresponds to the value estimated downstream of the termination shock assuming η_B = 5%, as explained in Sec. 3. Fig. E.1 shows a comparison of all relevant time scales: advection, diffusion, and energy losses for both protons and electrons. The maximum energy can be estimated by comparing the acceleration time scale with the combination of escape and loss times, namely $T_{\mathrm{acc}}^{-1}=T_{\mathrm{adv}}^{-1}+T_{\mathrm{diff}}^{-1}+T_{\mathrm{loss}}^{-1}$. This comparison is shown for the cases of Kraichnan and Bohm diffusion in Fig.E.1. One can see that the resulting electron maximum energy is a factor ≲2 (4) smaller than for protons in the Kraichnan (Bohm) case. To account for the energy losses during propagation inside the bubble, we use the same spatial profile as for protons (Eq. (C.2)-(C.4)) but rescaled by the ratio of the characteristic loss/escape timescale of the two species:

$$\begin{array}{c}\hfill f_e(r,p)=K_{\mathrm{ep}}{f}_p(r,p){\mathrm{e}}^{-p/p_{max,e}}\times \frac{{(T_{\mathrm{esc}}^{-1}+T_{\mathrm{IC}}^{-1}+T_{\mathrm{syn}}^{-1}+T_{\mathrm{brem}}^{-1})}_e^{-1}}{{(T_{\mathrm{esc}}^{-1}+T_{\mathrm{pp}}^{-1})}_p^{-1}},\end{array}$$(E.2)

where $T_{\mathrm{esc}}^{-1}=T_{\mathrm{age}}^{-1}+T_{\mathrm{adv}}^{-1}+T_{\mathrm{diff}}^{-1}$. This approach is an approximation to the real solution, however it is accurate enough for the purpose of this work.

Fig. E.1. The upper panel displays a comparison between the characteristic timescales of all relevant processes inside the bubble: advection, diffusion and energy losses for both protons and electrons. For diffusion, different types of turbulence are considered as detailed in the figure legend. The shaded area represent timescales larger of the age of the system (black line). The lower panel shows the maximum energy of electrons and protons for the Bohm (dashed) and the Kraichnan (solid) cases, obtained as the intersection between the acceleration time and the time-scale that results from the combination of escape ($T_{\mathrm{esc}}={(T_{\mathrm{adv}}^{-1}+T_{\mathrm{diff}}^{-1})}^{-1}$), energy losses (Tloss) and the age of the system Tage.

Once the electron spectrum is determined, we calculate the leptonic emission and we fit the observed gamma-ray data with a combination of Bremsstrahlung, IC and pion decay by changing the value of K_ep and η_CR. First, we consider the case where protons are completely negligible (this case is labeled K_ep ≫ 1) and try to fit the SED using only leptonic processes (top panel of Fig.E.2). We found that, when assuming Kolmogorov-type diffusion, we can in fact reproduce the data, but only assuming a very hard (α = 3.8) electron spectrum. The other two considered diffusion regimes do not allow a good fit of the data in any case. In addition, they require a very low maximum energy of the electrons in order not to exceed the highest energy upper limits. The only way to lower the maximum energy in a self-consistent way by invoking either a very small magnetic field, slowing down the acceleration rate, or a very large one, increasing the radiation loss rate. In the first case, the required magnetic field to obtain a low enough maximum energy must be of the order of ∼10 nG for both Kraichnan and Bohm, as one can easily see from Eq. (C.9). In the second case, by comparing the radiative timescales to the acceleration timescale, one finds that to lower the maximum energy enough a magnetic field of 100-200 μG is required: this exceeds the upper limit set in other works Mackey & Lim (2011), and would also imply η_B ≳ 1. The conclusion is that both the Kraichnan and Bohm cases require unreasonably low or high magnetic field strength.

In the main paper we disregarded the leptonic component when computing the acceleration efficiency, so we now look for the maximum leptonic contribution allowed by the data when it is combined with the hadronic component. In other words, we look for the maximum value of K_ep in Eq.(E.1) that allows us to fit the gamma-ray spectrum by reducing the proton normalization accordingly. To do that, we fix both the proton and electron spectral slope to the best-fit values reported in Table 1. The results are shown in the lower panel of Fig. E.2, for Kraichnan and Bohm diffusion. Kolmogorov, as argued in the main paper, does not reproduce the data, unless assuming a very-unlikely hard spectrum. We find that K_ep must be very small to avoid overshooting the data ($K_{\mathrm{ep}}^{max}={10}^{-3}(0.015)$ for Bohm (Kraichnan)). Even when maximizing the leptonic contribution, the normalization of protons needs to be reduced by a factor of 2 at most. We can therefore safely assume that the derived acceleration efficiency will not change more. Interestingly, in the Kraichnan case, the leptonic bremmstrahlung improve the agreement with the lowest energy data point.

Fig. E.2. Leptonic contribution to the SED of NGC 6611, considering two different ratio of electron-to-protons (Kep) at injection. In the upper panel the case where no protons are injected is considered (labeled as Kep ≫ 1), in the lower panel we consider the maximum allowed Kep values that do not overshoot the data. The emission is described in terms of bremsstrahlung and inverse Compton. More details on the modeling can be found in the text.

All Tables

All Figures

Fig. 1. Left: Test statistics map in the region around M16 derived from Fermi-LAT observations. The position of the O stars in the cluster NGC 6611 is indicated as blue stars, and its termination and forward shock radii are indicated a solid and dotted blue circles. The dashed circle indicates the limit of the H II region. The light blue, cyan, and pink contours indicate CO gas column densities of 1.5 × 1022, 1.7 × 1022, and 2.3 × 1022 cm−2, which were obtained in the velocity range 15–30 km s−1. The magenta circles indicate the position of Fermi sources from the 4FGL catalog Ballet (2023) and their relative position uncertainties, evaluated as their 68% confidence radius. Center: Molecular gas distribution in the region as traced by CO line emission, integrated over the line of sight in the velocity range 15 km/s, 30 km/s. The black contours refer to Fermi-LAT test statistics and represent 3, 4, and 5σ levels. Right: Velocity-longitude distribution of CO. The map is integrated over the whole latitude range. The dot indicates the position of NGC 6611.

In the text

	Fig. 3. Spectral energy distribution (SED) of the source detected by Fermi-LAT in correspondence of NGC 6611/M16. The error bars refer to the statistical 1σ uncertainties. The solid and dashed green lines represent the emission model in the case of Kraichnan turbulence, and the lines differ in the assumed wind luminosity, Lw, as indicated in the figure legend.
In the text

Fig. 2. Distribution of the gas density as a function of the radial distance from the center of the cluster (black curve) and the radial distribution of accelerated particles at different energies (solid, 10 GeV; dotted, 100 GeV), considering propagation in Bohm (blue), Kraichnan (green), or Kolmogorov (red) turbulence. The highest value for the wind luminosity, Lw, 2, was assumed in the calculations.

In the text

	Fig. A.1. Test statistic map of the region after modeling all sources except for 4FGL J1818.0-1339c. All 4FGL in the map are marked as magenta circles. The extension of the circle represents their localization uncertainty.
In the text

	Fig. C.1. Variation of Φ/Φ0 as a function of ηB and lC (see the text for a complete description). The white contour visible in the upper right panel corresponds to Φ/Φ0 = 0.5. In each panel the circle marks the reference parameters (lc = 1pc, ηB = 0.05)
In the text

	Fig. D.1. Parameter space explored with the fitting procedure. The darker zones are values with a smaller χ2, while the shaded regions represent those parameters excluded by the gamma-ray upper limits. The best fit value is indicated as circle, unless it is found in the forbidden parameter region, in which case the new best fit value is indicated as a cross.
In the text

	Fig. D.2. Three different fit results to the spectral energy distribution (SED) of the source detected by Fermi-LAT in correspondence of NGC 6611/M16. The solid and dashed curve differ in the assumed wind luminosity, Lw, as indicated in the figure legend. The three different panels refer to the three turbulence spectra assumed: Kolmogorov (red), Kraichnan (green), and Bohm (blue).
In the text

Fig. E.1. The upper panel displays a comparison between the characteristic timescales of all relevant processes inside the bubble: advection, diffusion and energy losses for both protons and electrons. For diffusion, different types of turbulence are considered as detailed in the figure legend. The shaded area represent timescales larger of the age of the system (black line). The lower panel shows the maximum energy of electrons and protons for the Bohm (dashed) and the Kraichnan (solid) cases, obtained as the intersection between the acceleration time and the time-scale that results from the combination of escape ($T_{\mathrm{esc}}={(T_{\mathrm{adv}}^{-1}+T_{\mathrm{diff}}^{-1})}^{-1}$), energy losses (Tloss) and the age of the system Tage.

In the text

Fig. E.2. Leptonic contribution to the SED of NGC 6611, considering two different ratio of electron-to-protons (Kep) at injection. In the upper panel the case where no protons are injected is considered (labeled as Kep ≫ 1), in the lower panel we consider the maximum allowed Kep values that do not overshoot the data. The emission is described in terms of bremsstrahlung and inverse Compton. More details on the modeling can be found in the text.

In the text