ABSTRACT

The Rosette Nebula is a young stellar cluster and molecular cloud complex, located at the edge of the southern shell of a middle-aged supernova remnant Monoceros Loop (G205.5+0.5). We revisited the GeV gamma-ray emission towards the Rosette Nebula using more than 13 yr of Fermi-LAT data. We tested several spatial models and found that compared to the result using the CO gas template only, the inclusion of the H ii gas template can significantly improve the likelihood fit. We performed spectral analysis using the new spatial template. With both the gamma-ray observation and CO+H ii gas data, we derived the cosmic ray spectrum of different components in the vicinity of the Rosette Nebula. We found the gamma-ray emissions from Rosette Nebula are substantially harder than previously reported, which may imply that Rosette Nebula is another example of a gamma-ray emitting young massive star cluster.

1 INTRODUCTION

As a complex of young massive clusters (YMCs) and molecular clouds (MCs), the Rosette Nebula is regarded as an important astrophysical laboratory for the studying early evolution of stars. A giant OB association NGC 2244 is located in the centre of the Rosette Nebula. It generates an expanding H ii region, which interacts with the surrounding the Rosette Molecular Cloud (RMC) and produces a photon-dominated region at their interface. To the west of NGC 2244, another star cluster NGC 2237 lies on the boundary between the H ii region and the cold neutral molecular material (Blitz & Thaddeus 1980; Li 2005; Román-Zúñiga et al. 2008). Bonatto & Bica (2009) suggested that NGC 2237 may be a young cluster located in the background of NGC 2244. The structure of NGC 2237 may contribute to the petal-like morphology of the Rosette Nebula (Wang et al. 2010). The ionized nebula interacts with the RMC, which is located to the east of the Rosette Nebula and has a collection of embedded young stellar clusters (Phelps & Lada 1997; Wang et al. 2009).

Along the line of sight, the Rosette Nebula partially overlaps with the edge of the southern shell of the Monoceros Loop (G205.5+0.5). Monoceros Loop is a well-studied middle-aged supernova remnant (SNR), and has a diameter of ∼3.8 (Katagiri et al. 2016), allowing detailed morphological studies in high-energy gamma-rays. By analysing about 5.5 yr of Fermi-LAT data, Katagiri et al. (2016) found extended gamma-ray emissions associated with the Monoceros Loop and the Rosette Nebula, and proposed that these gamma-rays have the same origin which is the interaction between the hadrons accelerated by the SNR and the interstellar gas.

However, we cannot formally rule out the different origins of gamma-rays in the Rosette Nebula and Monoceros SNR (Bourriche 2022). Indeed, in addition to SNRs, massive stars in YMCs were also proposed as major contributors to Galactic cosmic rays (CRs; Aharonian, Yang & de Oña Wilhelmi 2019). Recently, a series of observations found gamma-ray emissions that possibly originate from the CR nuclei accelerated in YMCs interacting with the surrounding gases, such as studies of Cygnus cocoon (Ackermann et al. 2011; Aharonian et al. 2019), Westerlund 1 (Abramowski et al. 2012), Westerlund 2 (Yang, de Oña Wilhelmi & Aharonian 2018), NGC 3603 (Yang & Aharonian 2017), 30 Dor C (H.E.S.S. Collaboration 2015), RSGC 1 (Sun, Yang & Wang 2020a), W40 (Sun et al. 2020b), Mc20 (Sun, Yang & Liang 2022), and NGC 6618 (Liu, Yang & Chen 2022). Furthermore, works have been done to determine the distances of NGC 2244 and the Rosette Nebula, and the results range from 1.4 to 1.7 kpc (Ogura & Ishida 1981; Hensberge, Pavlovski & Verschueren 2000; Park & Sung 2002). Recently, from the change of stellar extinction, Zhao et al. (2018) derived the distances to the Rosette Nebula and the Monoceros SNR are 1.55 and 1.98 kpc, respectively, implying no interaction between the nebula and the SNR. Thus, it is also possible that the gamma-rays in the Rosette Nebula originated from CRs accelerated by YMCs therein. With the accumulated Fermi-LAT data, we revisited the Rosette Nebula region to explore the origin of the possibly associated gamma-ray emissions to test such a hypothesis.

The paper is organized as follows. In Section 2, we described the gas content in this region. The procedure and results of data analysis are described in Section 3. The calculation of CR spectrum and density is described in Section 4. Finally, we present a discussion of our results and former works and draw our conclusions in Section 5.

2 GAS CONTENT AROUND THE ROSETTE NEBULA

According to the previous study, the gamma-rays around the nebula are generated from the hadronic interactions between CRs and interstellar gas (Katagiri et al. 2016). As a result, the spatial distribution of gamma-rays would generally follow that of gas. Thus, we first studied the spatial distribution of the gas content around the Rosette Nebula. For the following analysis, we adopted a distance of 1.55 kpc for the Rosette Nebula and the gas content around it.

Both molecular gas (mainly H2) and ionized gas (mainly H ii) are detected in the vicinity of the Rosette Nebula. For the molecular gas, we used the composite CO survey data from Dame, Hartmann & Thaddeus (2001) and derived the distribution of H2, by assuming a linear relationship between the velocity-integrated brightness temperature of the CO 2.6-mm line, WCO, and the column density of molecular hydrogen, N(H2). As proposed by Lebrun et al. (1983), N(H2) can be obtained as XCO × WCO, where the conversion factor XCO = 2.0 × 1020 cm−2 K−1 km−1 s (Dame et al. 2001; Bolatto, Wolfire & Leroy 2013). Our CO template is the same as the work of Katagiri et al. (2016), and we integrated the data over velocities of 0 km s−1 < V < 20 km s−1 based on the distance to the Rosette Nebula from Earth. The derived column density map of molecular gas is shown in the left panel of Fig. 1. The data in a circular region near the centre of Rosette with obvious flux excess was retained to construct a spatial distribution template (hereafter referred to as the CO template), which is shown by the cyan circle in the left panel of Fig. 1. The region we selected was slightly different from the model generated by Katagiri et al. (2016). Given the distribution of column density, the total mass of this component was calculated to be |$1.2\times 10^5~\rm M_\odot$|⁠.

Maps of gas column densities in two gas phases. Left shows the H2 column density derived from the CO data (Dame et al. 2001). Right shows the H ii column density derived from the Planck free–free map assuming the effective density of electrons $n_{\rm e}=10~\rm cm^{-3}$. The regions selected to generate the spatial distribution templates are shown by the cyan circles. For details, see the context in Section 2.
Figure 1.

Maps of gas column densities in two gas phases. Left shows the H2 column density derived from the CO data (Dame et al. 2001). Right shows the H ii column density derived from the Planck free–free map assuming the effective density of electrons |$n_{\rm e}=10~\rm cm^{-3}$|⁠. The regions selected to generate the spatial distribution templates are shown by the cyan circles. For details, see the context in Section 2.

As for the H ii column density, we used the Planck free–free map (Planck Collaboration X 2016) and applied the conversion factor from Finkbeiner (2003) to transform the emission measure into free–free intensity. And the H ii column density is derived via equation (5) from Sodroski et al. (1997):

(1)

where Iν is the intensity of free–free emission, |$\nu = \rm 353\ GHz$| is the frequency, and an electron temperature of |$T_{\mathrm{ e}} =\rm 8000\ K$| is applied. For ne, an effective density |$10\ \rm cm^{-3}$| is adopted, which is the value suggested in Sodroski et al. (1997). The derived H ii column density map is shown in the right panel of Fig. 1. For the spatial distribution template of this component (hereafter referred to as the H ii template), we chose the data in a circular region same as the one we chose for molecular gas, which is shown by the cyan circle in the right panel of Fig. 1. The total mass of this component was calculated to be |$7.3\times 10^4 ~\rm M_\odot$|⁠.

3 Fermi-LAT DATA ANALYSIS

To study the gamma-ray emission towards the Rosette Nebula, we used the latest Fermi-LAT Pass 8 data from 2008 August 4 (MET 239557417) until 2022 March 3 (MET 667979611), and used the fermitools from Conda distribution1 for the data analysis. In this work, we chose a 20× 20 square region centred at the position of the Rosette Nebula (R.A. = 97.979, Dec. = 4.942) as the region of interest (ROI). We used gtselect to select ‘source’ class events (evtype = 3 and evclass = 128) with zenith angles less than 90 to filter out the background contamination from the Earth’s limb. We also selected the good time intervals by applying the recommended expression |$\rm (DATA\_QUAL \gt 0) \&\& (LAT\_CONFIG == 1)$|⁠. The instrument response functions P8R3_SOURCE_V3 were applied to analyse the source events. We performed a standard binned analysis following the official tutorial of binned likelihood analysis.2 We generated the source models using make4FGLxml.py,3 which consists of the sources in the LAT 12-yr Source Catalogue (4FGL-DR3, Abdollahi et al. 2020) those are within the ROI enlarged by 10, the Galactic diffuse background emission (gll_iem_v07.fits), and the isotropic emission background (iso_P8R3_SOURCE_V3_v1.txt). In addition, the spectral parameters of sources within 10 from the centre and the normalization factor of the Galactic diffuse background and isotropic background were set free.

3.1 Spatial analysis

First, we used the events from 0.1 to 500 GeV to study the spatial distribution of the gamma-ray emission near the nebula. The gamma-ray counts map in the 5 × 5 region around the Rosette Nebula is shown in Fig. 2. In the 4FGL-DR3 catalogue, the spatial models of Rosette Nebula and Monoceros Loop are adopted from the work of Katagiri et al. (2016), in which the nebula represented by a template derived from the CO gas distribution and the SNR is fitted as a Gaussian disc. Besides, there are three point-like sources in 4FGL-DR3 overlapping with the Rosette Nebula along the line of sight, each with a corresponding identified source (4FGL J0632.8+0550 and HESS J0632+057, 4FGL J0633.7+0632 and PSR J0633+0632, 4FGL J0631.8+0645 and PSR J0631+0646).

Gamma-ray counts map of around the Rosette Nebula. The magenta dashed circle shows the position of the Monoceros Loop. The point-like sources in 4FGL-DR3 are shown in green crosses. The green diamonds represent the positions of seven embedded clusters in the RMC identified by Phelps & Lada (1997). The cyan box shows the position of NGC 2237 and NGC 2244 (Román-Zúñiga et al. 2008). The black contours show the spatial distribution of the H2 gas while the brown contours show the distribution of the H ii gas.
Figure 2.

Gamma-ray counts map of around the Rosette Nebula. The magenta dashed circle shows the position of the Monoceros Loop. The point-like sources in 4FGL-DR3 are shown in green crosses. The green diamonds represent the positions of seven embedded clusters in the RMC identified by Phelps & Lada (1997). The cyan box shows the position of NGC 2237 and NGC 2244 (Román-Zúñiga et al. 2008). The black contours show the spatial distribution of the H2 gas while the brown contours show the distribution of the H ii gas.

With the sources in the 4FGL-DR3 catalogue, we conducted a binned likelihood analysis and applied the fitted 4FGL-DR3 model to calculate the residual test statistic (TS) map around the nebula. The TS value is defined as |$2\left(\log \mathcal {L}/\mathcal {L}_0\right)$|⁠, in which |$\mathcal {L}_0$| is the likelihood of null hypothesis and |$\mathcal {L}$| is the likelihood with a test source included at a given position. However, evident excess around the nebula was found from the residual TS map, indicating a more detailed analysis is required to study the gamma-ray emission from the Rosette Nebula. We then replaced the Rosette template in the 4FGL-DR3 model with the CO template which is described in Section 2, and applied the log-parabolic spectral function for the CO template, which is the same setting for the Rosette template in the 4FGL-DR3 catalogue (hereafter referred to as the CO model). The analysis results of the CO model are basically the same as that of the original 4FGL-DR3 model since both the CO template and the Rosette template in 4FGL-DR3 are derived from the same CO data (Katagiri et al. 2016). Both models obtained similar −|$\log ({\cal L})$| of about −2453860 and their residual TS maps show a similar distribution of the excess emission around the nebula. From the residual TS map of the CO model (as shown in the left panel of Fig. 3), we can see the distribution of residual gamma-rays mainly overlaps with the distribution of the H ii gas.

Residual gamma-ray TS maps obtained from different source models for energy range 0.1–500 GeV, smoothed with Gaussian filter of 0.3○. The black and brown contours show the spatial distribution of the CO template and the H ii template, respectively.
Figure 3.

Residual gamma-ray TS maps obtained from different source models for energy range 0.1–500 GeV, smoothed with Gaussian filter of 0.3. The black and brown contours show the spatial distribution of the CO template and the H ii template, respectively.

As we studied the H ii density in the same region in Section 2, the column density of H ii was found to be of the same order of magnitudes with molecular gas, and with quite a different distribution seemingly related to the excess in the residual TS map obtained using the CO model (see the left panel of Fig. 3). Thus, we added a new source of which the spatial model is the H ii template and the spectral function is log-parabola to the CO model (hereafter referred to as CO+H ii model), and performed a new likelihood analysis. As shown in the TS map derived from the CO+H ii model (the right panel of Fig. 3), the residual emission is negligible after adding the H ii component. Moreover, we also tried to use a radial Gaussian disc with the spectral function of log-parabola to represent the residual emission in the left panel of Fig. 3. The disc was centred at the TS peak position (R.A. = 98.531, Dec. = 4.992), with a sigma ranging from 0.1 to 0.9. The best-fitting result is the disc with a sigma of 0.2 (hereafter referred to as the Gaussian model). To find which model fits the data best, we applied the Akaike information criterion (AIC). AIC is calculated as −2log (likelihood) + 2k, in which k is the number of free parameters in the model. The −|$\log ({\cal L})$| and corresponding AIC value of the models are listed in Table 1. As shown in Table 1, the best fit for the data is the CO+H ii model.

Table 1.

Fermi-LAT data fitting results of different models.

ModelFree Parameters (k)|$\log ({\cal L})$|AIC
CO model77−2453860−4907566
CO+H ii model80−2453911−4907662
Gaussian model80−2453858−4907556
ModelFree Parameters (k)|$\log ({\cal L})$|AIC
CO model77−2453860−4907566
CO+H ii model80−2453911−4907662
Gaussian model80−2453858−4907556
Table 1.

Fermi-LAT data fitting results of different models.

ModelFree Parameters (k)|$\log ({\cal L})$|AIC
CO model77−2453860−4907566
CO+H ii model80−2453911−4907662
Gaussian model80−2453858−4907556
ModelFree Parameters (k)|$\log ({\cal L})$|AIC
CO model77−2453860−4907566
CO+H ii model80−2453911−4907662
Gaussian model80−2453858−4907556

3.2 Spectral analysis

Based on the above spatial analysis, we derived the spectral energy distribution (SED) of the Rosette Nebula of the CO model and the CO+H ii model respectively. Here, we divided the data into nine logarithmically spaced energy bins and obtained the SED by applying maximum-likelihood analysis for each bin. Systematic uncertainties of the flux for each energy bin were also estimated using the method described on the Fermi official website.4 Besides, we calculated the 99 per cent flux upper limit for energy bins with the source TS value less than 4. The derived SEDs of different components in different models are shown in Fig. 4.

Gamma-ray spectral energy distributions of the Rosette Nebula and the Monoceros SNR obtained from different source models. For each data point, the error bar indicates the combined statistical and systematic uncertainty. Data points with TS values lower than 4 are replaced by their upper limits. The dashed lines show the corresponding pion-decay emission derived from the best-fitting CR spectra (see Section 4 for details.).
Figure 4.

Gamma-ray spectral energy distributions of the Rosette Nebula and the Monoceros SNR obtained from different source models. For each data point, the error bar indicates the combined statistical and systematic uncertainty. Data points with TS values lower than 4 are replaced by their upper limits. The dashed lines show the corresponding pion-decay emission derived from the best-fitting CR spectra (see Section 4 for details.).

4 CR CONTENT IN THE VICINITY OF THE ROSETTE NEBULA

We assumed that the gamma-rays are mainly generated from the proton–proton interaction between CR protons and the gases, which enabled us to derive the spectra of CR in the vicinity of the Rosette Nebula using the information of the gamma-ray spectra as well as the derived gas masses. We used the gamma-ray production cross-section from Kafexhiu et al. (2014) and assumed a power-law proton spectrum F = fiEΓ to perform a likelihood fitting to the gamma-ray data. Here, F is the flux density of protons in unit of |$\mathrm{cm}^{-2}\, \mathrm{s}^{-1}\, \mathrm{GeV}^{-1}$|⁠, E is the proton energy in unit of GeV, and Γ stands for the spectral index. Given the gas masses obtained in Section 2, we derived the absolute CR fluxes associated with the CO and H ii components, respectively. Especially, for Monoceros, we used an average gas density of 3.6 cm−3 based on the H i observations (Xiao & Zhu 2012). The best-fitting results are listed in Table 2 and illustrated in Fig. 5. The spectral index of CR protons associated with the CO gas is a bit harder after adding the H ii template and is much harder than that of the Monoceros SNR. Besides, the dashed lines in Fig. 4 represent the corresponding gamma-ray emissions derived from the CR fluxes and gas masses. As shown in Fig. 5, the fluxes of CR protons that are associated with the nebula and the SNR are higher than the local CR proton flux measured by AMS-02 (Aguilar et al. 2015).

Derived CR spectra of the Rosette Nebula and the Monoceros SNR for different models. The dashed lines show the best-fitting CR spectra and the shaded regions represent the 1σ errors including both the statistics and the systematic uncertainty. The black dots are the AMS-02 observation data acquired from Aguilar et al. (2015).
Figure 5.

Derived CR spectra of the Rosette Nebula and the Monoceros SNR for different models. The dashed lines show the best-fitting CR spectra and the shaded regions represent the 1σ errors including both the statistics and the systematic uncertainty. The black dots are the AMS-02 observation data acquired from Aguilar et al. (2015).

Table 2.

Fitted spectral parameters of CR proton for different gas contents of the Rosette Nebula and the Monoceros SNR.

ModelComponent|$f_i(10^{-10}\, \mathrm{cm}^{-2}\, \mathrm{s}^{-1}\, \mathrm{GeV}^{-1})$|Γ
CO modelCO8.59 ± 0.35−2.46 ± 0.02
Monoceros93.17 ± 6.23−2.62 ± 0.04
CO+H ii modelCO3.57 ± 0.30−2.30 ± 0.03
H ii7.63 ± 1.24−2.46 ± 0.10
Monoceros94.72 ± 10.04−2.68 ± 0.04
ModelComponent|$f_i(10^{-10}\, \mathrm{cm}^{-2}\, \mathrm{s}^{-1}\, \mathrm{GeV}^{-1})$|Γ
CO modelCO8.59 ± 0.35−2.46 ± 0.02
Monoceros93.17 ± 6.23−2.62 ± 0.04
CO+H ii modelCO3.57 ± 0.30−2.30 ± 0.03
H ii7.63 ± 1.24−2.46 ± 0.10
Monoceros94.72 ± 10.04−2.68 ± 0.04
Table 2.

Fitted spectral parameters of CR proton for different gas contents of the Rosette Nebula and the Monoceros SNR.

ModelComponent|$f_i(10^{-10}\, \mathrm{cm}^{-2}\, \mathrm{s}^{-1}\, \mathrm{GeV}^{-1})$|Γ
CO modelCO8.59 ± 0.35−2.46 ± 0.02
Monoceros93.17 ± 6.23−2.62 ± 0.04
CO+H ii modelCO3.57 ± 0.30−2.30 ± 0.03
H ii7.63 ± 1.24−2.46 ± 0.10
Monoceros94.72 ± 10.04−2.68 ± 0.04
ModelComponent|$f_i(10^{-10}\, \mathrm{cm}^{-2}\, \mathrm{s}^{-1}\, \mathrm{GeV}^{-1})$|Γ
CO modelCO8.59 ± 0.35−2.46 ± 0.02
Monoceros93.17 ± 6.23−2.62 ± 0.04
CO+H ii modelCO3.57 ± 0.30−2.30 ± 0.03
H ii7.63 ± 1.24−2.46 ± 0.10
Monoceros94.72 ± 10.04−2.68 ± 0.04

5 DISCUSSION AND CONCLUSION

In this paper, we analysed Fermi-LAT data towards the Rosette Nebula with a time span of more than 13 yr. We found that the extended gamma-ray emission in the vicinity of the nebula can be well modelled by the CO gas distributions, which is consistent with the former results (Katagiri et al. 2016). However, we found that adding the H ii gas template can significantly improve the fitting results.

After including this H ii gas template, the derived gamma-ray emissions associated with the Rosette Nebula (CO template and H ii template) and the SNR Monoceros show clear low-energy break below |$1~\rm GeV$|⁠, which are strong hints that the gamma-ray emissions from both regions are from the pion-decay process in the interaction of CR protons with ambient gas. Such results are also similar to that in the former study (Katagiri et al. 2016). However, the gamma-ray emission from the Rosette Nebula is now significantly harder than that from Monoceros. And the CR associated with the ionized gas has a similar index (about −2.46) to that associated with the molecular gas (about −2.30), harder than that of the Monoceros SNR (from −2.6 to −2.7). This is substantially different from the results in the former study in Katagiri et al. (2016), in which the similar spectral shape in both regions is used implying that the gamma-ray emissions from both regions are from the same CR population accelerated by the SNR Monoceros.

The different gamma-ray spectra in Rosette and Monoceros cannot formally rule out the possibility that the gamma-ray emissions in the nebula are produced by the CRs that are accelerated and injected by the SNR. The harder spectrum in the nebula may be formed by propagation effects, that is, due to the energy-dependent diffusion only the high-energy CRs reached the Rosette Nebula. Assuming the Monoceros SNR and the Rosette Nebula are located at the same distance (1.55 kpc), the distance between the centre of Monoceros SNR and Rosette Nebula is about |$60~\rm pc$|⁠. In case the harder spectrum of the Rosette Nebula is attributed to that the CR below about |$10~\rm GeV$| has not reached Rosette Nebula yet and taken into account the age of the SNR of 3 × 104 yr (Katagiri et al. 2016), the required diffusion coefficient is about |$10^{28}~\rm cm^2/s$| at |$10~\rm GeV$|⁠, which is about one order of magnitude smaller than the value in the interstellar medium (Strong & Moskalenko 1998). Such a low diffusion coefficient is also observed in other SNR–MC systems (e.g. Fujita, Ohira & Takahara 2010; Li & Chen 2012; Liu et al. 2015), which may be caused by a higher turbulence level near the CR accelerators. However, if the distance between the SNR and the nebula is ∼0.4 kpc as the recent measurement made by Zhao et al. (2018), the diffusion coefficient is about |$10^{30}~\rm cm^2\, s^{-1}$| at 10 GeV for the protons accelerated by the SNR to reach the nebula, which is about ten times larger than the value in the interstellar medium.

On the other hand, it is also possible that the gamma-rays from the Rosette Nebula are produced by CRs other than those accelerated by Monoceros SNR. The recent estimation of the distance of the Rosette Nebula also suggests that it doesn’t interact with the SNR (Zhao et al. 2018). A natural acceleration site is the star cluster NGC 2244 inside the Rosette Nebula. NGC 2244 is a YMC with an age of about |$2~\rm Myr$| (Mužić et al. 2022), and harbours more than 20 stars with a spectral type early than B3 (Wang et al. 2009). The total wind power can be estimated to be at the order of |$10^{37} ~\rm erg\, s^{-1}$|⁠. Assuming an acceleration efficiency of 10 per cent, the total CR accelerated in NGC 2244 can be estimated as |$10^{50}~\rm erg$|⁠. The gamma-ray luminosity in the Rosette Nebula is about |$10^{34}~\rm erg\, s^{-1}$|⁠, taken into account of the gas density of about |$10~\rm cm^{-3}$|⁠, and the total energy in the CR content is |${\sim} 10^{49} ~\rm erg$|⁠. Therefore, NGC 2244 is powerful enough to account for the gamma-ray emissions in the Rosette Nebula. Furthermore, the hard gamma-ray spectrum is also similar to the gamma-ray emission in other YMC systems, such as Cygnus cocoon (Ackermann et al. 2011; Aharonian et al. 2019), NGC 3603 (Yang & Aharonian 2017), W40 (Sun et al. 2020b), W43 (Yang & Wang 2020), NGC 6618 (Liu et al. 2022).

In conclusion, we revisited the GeV gamma-ray emission towards the Rosette Nebula and found the inclusion of an H ii gas template can significantly improve the likelihood of the fitting results. With the new spatial template, we found the gamma-ray emissions from the Rosette Nebula are substantially harder than previously reported. Although we cannot rule out the possibility that the gamma-ray emission from this region is produced by CRs accelerated by the nearby SNR Monoceros. A more natural explanation is that the CR accelerated by YMC NGC 2244 illuminated both molecular and ionized gases in the Rosette Nebula. In this case, the Rosette Nebula is another example of the gamma-ray-emitting YMC system.

ACKNOWLEDGEMENTS

RzY is supported by the NSFC under grant 12041305 and the national youth thousand talents programme in China. BL acknowledges the support from the NSFC under grant 12103049.

DATA AVAILABILITY

The Fermi-LAT data used in this work is publicly available, which is provided online by the NASA-GSFC Fermi Science Support Center5. We make use of the CO data6 to derive the H2. The data from Planck legacy archive7 are used to derive the H ii gas content. The H i data are from the HI4PI8.

Footnotes

References

Abdollahi
S.
et al. ,
2020
,
ApJS
,
247
,
33

Abramowski
A.
et al. ,
2012
,
A&A
,
537
,
A114

Ackermann
M.
et al. ,
2011
,
Science
,
334
,
1103

Aguilar
M.
et al. ,
2015
,
Phys. Rev. Lett.
,
114
,
171103

Aharonian
F.
,
Yang
R.
,
de Oña Wilhelmi
E.
,
2019
,
Nat. Astron.
,
3
,
561

Blitz
L.
,
Thaddeus
P.
,
1980
,
ApJ
,
241
,
676

Bolatto
A. D.
,
Wolfire
M.
,
Leroy
A. K.
,
2013
,
ARA&A
,
51
,
207

Bonatto
C.
,
Bica
E.
,
2009
,
MNRAS
,
394
,
2127

Bourriche
N.
,
2022
,
Star forming regions as hadronic cosmic ray accelerators
,
Master’s thesis
,
Leopold-Franzens Universität
,
Innsbruck

Dame
T. M.
,
Hartmann
D.
,
Thaddeus
P.
,
2001
,
ApJ
,
547
,
792

Finkbeiner
D. P.
,
2003
,
ApJS
,
146
,
407

Fujita
Y.
,
Ohira
Y.
,
Takahara
F.
,
2010
,
ApJ
,
712
,
L153

H.E.S.S. Collaboration
,
2015
,
Science
,
347
,
406

Hensberge
H.
,
Pavlovski
K.
,
Verschueren
W.
,
2000
,
A&A
,
358
,
553

Kafexhiu
E.
,
Aharonian
F.
,
Taylor
A. M.
,
Vila
G. S.
,
2014
,
Phys. Rev. D
,
90
,
123014

Katagiri
H.
et al. ,
2016
,
ApJ
,
831
,
106

Lebrun
F.
et al. ,
1983
,
ApJ
,
274
,
231

Li
H.
,
Chen
Y.
,
2012
,
MNRAS
,
421
,
935

Li
J. Z.
,
2005
,
ApJ
,
625
,
242

Liu
B.
,
Chen
Y.
,
Zhang
X.
,
Zhang
G.-Y.
,
Xing
Y.
,
Pannuti
T. G.
,
2015
,
ApJ
,
809
,
102

Liu
B.
,
Yang
R.-z.
,
Chen
Z.
,
2022
,
MNRAS
,
513
,
4747

Mužić
K.
,
Almendros-Abad
V.
,
Bouy
H.
,
Kubiak
K.
,
Peña Ramírez
K.
,
Krone-Martins
A.
,
Moitinho
A.
,
Conceição
M.
,
2022
,
A&A
,
668
,
A19

Ogura
K.
,
Ishida
K.
,
1981
,
PASJ
,
33
,
149

Park
B.-G.
,
Sung
H.
,
2002
,
AJ
,
123
,
892

Phelps
R. L.
,
Lada
E. A.
,
1997
,
ApJ
,
477
,
176

Planck Collaboration X
,
2016
,
A&A
,
594
,
A10

Román-Zúñiga
C. G.
,
Elston
R.
,
Ferreira
B.
,
Lada
E. A.
,
2008
,
ApJ
,
672
,
861

Sodroski
T. J.
,
Odegard
N.
,
Arendt
R. G.
,
Dwek
E.
,
Weiland
J. L.
,
Hauser
M. G.
,
Kelsall
T.
,
1997
,
ApJ
,
480
,
173

Strong
A. W.
,
Moskalenko
I. V.
,
1998
,
ApJ
,
509
,
212

Sun
X.-N.
,
Yang
R.-Z.
,
Liang
E.-W.
,
2022
,
A&A
,
659
,
A83

Sun
X.-N.
,
Yang
R.-Z.
,
Wang
X.-Y.
,
2020a
,
MNRAS
,
494
,
3405

Sun
X.-N.
,
Yang
R.-Z.
,
Liang
Y.-F.
,
Peng
F.-K.
,
Zhang
H.-M.
,
Wang
X.-Y.
,
Aharonian
F.
,
2020b
,
A&A
,
639
,
A80

Wang
J.
,
Feigelson
E. D.
,
Townsley
L. K.
,
Broos
P. S.
,
Román-Zúñiga
C. G.
,
Lada
E.
,
Garmire
G.
,
2010
,
ApJ
,
716
,
474

Wang
J.
,
Feigelson
E. D.
,
Townsley
L. K.
,
Román-Zúñiga
C. G.
,
Lada
E.
,
Garmire
G.
,
2009
,
ApJ
,
696
,
47

Xiao
L.
,
Zhu
M.
,
2012
,
A&A
,
545
,
A86

Yang
R.-z.
,
Aharonian
F.
,
2017
,
A&A
,
600
,
A107

Yang
R.-z.
,
de Oña Wilhelmi
E.
,
Aharonian
F.
,
2018
,
A&A
,
611
,
A77

Yang
R.-Z.
,
Wang
Y.
,
2020
,
A&A
,
640
,
A60

Zhao
H.
,
Jiang
B.
,
Gao
S.
,
Li
J.
,
Sun
M.
,
2018
,
ApJ
,
855
,
12

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)