Detection of diffuse gamma-ray emission towards a massive star forming region hosting Wolf-Rayet stars

Isotopic and elemental abundances seen in Galactic cosmic rays imply that $\sim20\%$ of the cosmic-ray (CR) nuclei are probably synthesized by massive Wolf-Rayet (WR) stars. Massive star clusters hosting WR and OB-type stars have been proposed as potential Galactic cosmic-ray accelerators for decades, in particular via diffusive shock acceleration at wind termination shocks. Here we report the analysis of {\em Fermi} Large Area Telescope's data towards the direction of Masgomas-6a, a young massive star cluster candidate hosting two WR stars. We detect an extended $\gamma$-ray source with $\rm{TS}=183$ in the vicinity of Masgomas-6a, spatially coincident with two unassociated {\em Fermi} 4FGL sources. We also present the CO observational results of molecular clouds in this region, using the data from the Milky Way Imaging Scroll Painting project. The $\gamma$-ray emission intensity correlates well with the distribution of molecular gas at the distance of Masgomas-6a, indicating that these gamma rays may be produced by CRs accelerated by massive stars in Masgomas-6a. At the distance of $3.9{\rm \ kpc}$ of Masgomas-6a, the luminosity of the extended source is $(1.81\pm0.02)\times 10^{35}{\rm \ erg \ s^{-1}}$. With a kinetic luminosity of $\sim 10^{37}{\rm erg \ s^{-1}}$ in the stellar winds, the WR stars are capable of powering the $\gamma$-ray emission via neutral pion decay resulted from cosmic ray $pp$ interactions. The size of the GeV source and the energetic requirement suggests a CR diffusion coefficient smaller than that in the Galactic interstellar medium, indicating strong suppression of CR diffusion in the molecular cloud.


INTRODUCTION
The elemental composition of Galactic cosmic rays (GCRs) has been investigated with several experiments (IMP7 (Garcia-Munoz et al. 1979), ISEE-3 (Wiedenbeck & Greiner 1981), Voyager (Lukasiak et al. 1994), ACE-CRIS (Binns et al. 2005), DAMPE , AMS (Aguilar et al. 2020) and etc.). For most elements, the isotopic ratios seen in GCRs are similar to those in the solar wind (Wiedenbeck et al. 1979). The most notable exception to this is 22 Ne/ 20 Ne, where the CR value is about 5 times that of the solar wind (Binns et al. 2005(Binns et al. , 2008. According to the most recent models of nucleosynthesis, a large amount of 22 Ne is generated in Wolf-Rayet (WR) stars (Higdon & Lingenfelter 2003;Kalyashova et al. 2019). The interacting winds of massive stars have been recognized as potential CR accelerators as early as in the 1980s (Bykov & Toptygin 1982;Cesarsky & Montmerle 1983). The acceleration could take place in the vicinity of the stars (Bykov 2014;Aharonian et al. 2019), in the colliding wind region of WR and OB binary sys-tems (Benaglia & Romero 2003), or in superbubbles, multi-parsec structures caused by the collective activity of massive stars (Parizot et al. 2004;Bykov 2014). Higdon & Lingenfelter (2003) found that the observed value of 22 Ne/ 20 Ne can be achieved with a mixture of ∼ 20% WR star material and ∼ 80% material with standard composition. Using state-of-the-art stellar evolution models, Kalyashova et al. (2019) calculated the amount of neon isotopes produced by interacting winds of OB-type and WR stars and also concluded that massive star clusters can produce a significant fraction of GCRs and 22 Ne.
Spectroscopic data show that WR stars, along with OB-type stars, have powerful stellar winds with velocities of 1000-3000 km/s (Hamann et al. 2019). Due to their small size and high star density (Portegies Zwart et al. 2010), young massive star clusters are likely to have strong colliding magnetohydrodynamic shock flows, where efficient particle acceleration can take place. Recently, Seo et al. (2018) estimated the total power of OB-type and WR stars' winds in the Galaxy as ∼ 1.1 × 10 41 erg s −1 . If (1 − 10)% of the wind luminosity is converted to GCRs, these stellar winds can provide a significant contribution to the GCR production below the knee. The conversion of stellar wind power to CRs might be traced by secondary γ-rays, the products of interactions of CRs with the surrounding gas. The diffuse GeV γ-rays detected by Fermi-LAT telescope around some compact clusters, such as Cygnus OB2 (Ackermann et al. 2011), NGC 3603 (Yang & Aharonian 2017), Westerlund 2 (Yang et al. 2018) and RSGC 1 (Sun et al. 2020), can be naturally interpreted within this scenario. Detection of TeV or higher energy photons from some massive star clusters (e.g., the Cygus region) may suggest that cosmic rays are accelerated up to energy of 100 TeV to PeV (Aharonian et al. 2002;Abeysekara et al. 2021;Amenomori et al. 2021) Recently, Ramírez Alegría et al. (2018) identified two groups of massive stars aligned in the l ∼ 38 • Galactic direction. They find more than 20 massive stars, including two Wolf-Rayets (WR122-11 and the new WR122-16), OB-type stars and one transitional object (the Luminous Blue Variable (LBV) candidate IRAS 18576+0341 (Pasquali & Comerón 2002)). The individual distances and radial velocities of these massive stars indicate two populations of massive stars in the same line of sight: Masgomas-6a and Masgomas-6b. Masgomas-6a is reported as a massive star cluster candidate with a total mass (lower limit) of (9.0 − 1.3) × 10 3 M ⊙ , located at 3.9 +0.4 −0.3 kpc, contains both Wolf-Rayets and several OB-type stars. Masgomas-6b, located at 9.6 ± 0.4 kpc, has a total mass (lower limit) of (10.5 − 1.5) × 10 3 M ⊙ and hosts a LBV candidate and an evolved population of supergiants. For both objects, the presence of evolved massive stars (the WN8-9h Wolf-Rayet in Masgomas-6a and the LBV candidate in Masgomas-6b) sets an age limit of 5 Myr (Ramírez Alegría et al. 2018).
In this paper, we report the detection of an extended γ-ray source towards the direction of Masgomas-6 with the Fermi Large Area Telescope. We also present the CO observational results of molecular clouds (MCs) in this region, using the data from the Milky Way Imaging Scroll Painting (MWISP) project, which is a multi-line survey in 12 CO/ 13 CO/C 18 O along the northern Galactic plane with the Purple Mountain Observatory (PMO) 13.7 m telescope. The correlation between the γray intensity and the molecular gas at the distance of Masgomas-6a suggests that the GeV γ-ray source is likely related with the massive stars in Masgomas-6a. The paper is organized as follows. We first present the Fermi/LAT data analysis result of the γ-ray source in §2. We then present the MWISP observational results of molecular clouds in this region in §3. In §4, the interpretation of the GeV source and the implication is presented. Finally, we give a summary in §5.

FERMI /LAT DATA ANALYSIS
The Fermi-LAT is a wide field-of-view, high-energy γ-ray telescope, covering the energy range from below 20 MeV to more than 300 GeV, and it has continuously monitored the sky since 2008 (Atwood et al. 2009). The Pass 8 data taken from 2008 August 4 to 2022 April 9 are used to study the GeV emission towards Masgomas-6, where two Fermi sources (4FGL J1900.4+0339 and 4FGL J1858.0+0354) are found in the vicinity. All γ-ray photons within a 14 • × 14 • region of interest (ROI) centered on the position of Masgomas-6 are considered for the binned maximum likelihood analysis. The publicly available software Fermipy (ver.1.0.1) and Fermitools (ver.2.0.8) are used to perform the data analysis. The event type FRONT + BACK and the instrument response functions (IRFs) (P 8R3 SOU RCE V 3) are used. We only consider the γ-ray events in the 3 − 500 GeV energy range, with the standard data quality selection criteria "(DAT A QU AL > 0)&&(LAT CON F IG == 1)". To minimize the contamination from γ-rays from the Earth limb, the maximum zenith angle is adopted to be 90 • . We include all sources listed in the fourth Fermi-LAT catalog (Abdollahi et al. 2020) and the diffuse Galactic interstellar emission (IEM, gll iem v07.f its), isotropic emission ("iso P 8R3 SOU RCE V 3 v1.txt" ) in the background model. The spectral parameters of the sources within 4 • of Masgomas-6, Galactic and isotropic diffuse emission components are all set free. The significance of γ-ray sources are estimated by maximum likelihood test statistic (TS), which is defined by TS= −2(ln L 0 − ln L 1 ), where L 1 and L 0 are maximum likelihood values for the background with target source and without the target source (null hypothesis). The square root of the TS is approximately equal to the detection significance of a given source.

Morphological analysis
To obtain a sufficiently good angular resolution for the morphology analysis, we use the LAT data in 3-500 GeV. Figure 1 shows the 3-500 GeV TS map in the the vicinity of Masgomas-6 with the binned likelihood analysis. The two 4FGL point-like sources near Masgomas-6, 4FGL J1900.4+0339 and 4FGL J1858.0+0354, are located 0.16 • and 0.52 • away from the center of Masgomas-6, respectively. We plan to investigate what appear to be diffuse emissions beyond the two 4FGL point-like sources. Since these two point-like sources are unassociated with any known counterparts and there is diffuse emission in the vicinity, we investigate whether the γ-ray emission in this region represents a spatially extended source. We use the Fermipy tool to quantitatively evaluate the extension and location of this γ-ray emission. A uniform disk model and a two-dimensional (2D) Gaussian model are used to evaluate the extension. The results are shown in Table 1. We find that the uniform disk model describes well the GeV morphology, with the best-fit extension of 0.43 •+0.02 • −0.03 • and the best-fit position at (R.A., Decl.)=(284.81 • ± 0.08 • ,3.84 • ± 0.04 • ) in the energy band above 3 GeV. This position is 0.22 • away from the position of Masgomas-6. The extension represents the radius containing 68% of the intensity. The TS ext is estimated to be 129.80 for the disk model, which is defined as TS ext = 2(ln L ext − ln L ps ), where L ext is the maximum likelihood value for the extend model and L ps is the maximum likelihood value for the point-like model (Ackermann et al. 2017). According to Ackermann et al. (2017), the criterion to define a source as being extended is TS ext ≥ 16 and TS ext ≥ TS 2pts (here TS 2pts = 2(ln L 2pts − ln L ps ) represents the improvement when adding a second point source). We find that the γ-ray source, assuming a uniform disk model, meets the criterion with TS ext = 129.80 and TS 2pts = 106.62 (see Table 1). We also test the 2D Gaussian model and find that it is also better than the the point-like model. Therefore, the GeV emission in this region is described better by an extended source.
Furthermore, we compare the residual TS maps after subtracting the disk component and the two point-like source component. We find that there is bright (maximum TS ∼ 25) emission (referred to as Src A, see Figure 1) at the top left corner of the residual TS map in the two point-like source model, but the significance of this emission is significantly reduced in the disk model. This implies that the extended source model fits the GeV emission better than the two point-like sources model, consistent with our above analysis.

Energy Spectrum
After the morphology was fixed, the extended disk model is used to study the spectrum of the source. The spectral energy distribution (SED) of the source in the energy band > 0.1 GeV, shown in Figure 2, is derived by gtlike assuming the best-fit uniform disk extension. When the TS value of spectra data point is less than 4, an upper limit is calculated at 95% confidence level using a Bayesian method (Helene 1983). Compared to a single power-law, the spectrum is better described by a log parabola empirical function . The spectral index and break energy are found to be α = 2.29 ± 0.01, β = 0.19 ± 0.01 and E b = 1268.47 ± 2.81MeV. The energy flux is (9.92 ± 0.12) × 10 −11 erg cm −2 s −1 in 0.1-500 GeV. Assuming a distance of d = 3.9 kpc, the total γ-ray luminosity of the source in 0.1-500 GeV is (1.81 ± 0.02) × 10 35 erg s −1 . The preference of a curved spectrum over a single power-law is supported by TS curve = 102.75 (10.14σ), where TS curve is defined as TS curve = 2(ln L curved spectrum −ln L power−law ). In the energy band greater than 3 GeV, however, the spectra of the disk model is well described by a single power-law function with TS curve < 9 (TS curve = 9 corresponding to 3σ (Abdollahi et al. 2020)).

CO OBSERVATIONS
We make use of the data from the Milky Way Imaging Scroll Painting (MWISP 1 ) project, which is a multiline survey in 12 CO/ 13 CO/C 18 O observed simultaneously using the 13.7 m millimeter-wavelength telescope of the Purple Mountain Observatory at Delingha. The detailed observing strategy, the instrument, the spectral resolution, and the quality of the CO observations can be found in Su et al. (2019). In this work, we present the results of the MWISP CO survey for 2 • × 2 • regions centering at (l, b) = (37.27 • , −0.23 • ). We estimated the kinematic distance to the MCs using the Milky Way's rotation curve suggested by Sofue (2015), assuming the Sun's Galactocentric distance to be 8.0 kpc and Solar rotation velocity to be 238 km s −1 .
The distance of Masgomas-6a is estimated to be 3.9 +0.4 −0.3 kpc (Ramírez Alegría et al. 2018), corresponding to a velocity V LSR in the range of ∼ 60 − 70 km s −1 for the MC CO emission. Note that, the velocity of each MC could indicates two candidate kinematic distances, the near side one and the far side one. The far distance of the MC in ∼ 60 − 70 km s −1 is 8.4-9.0 kpc. The distance of Masgomas-6b is estimated to be 9.6 ± 0.4 kpc, corresponding to a velocity V LSR in the range of ∼ 45 − 59 km s −1 . The 12 CO(J = 1-0) maps for the velocity range 60 to 70 km s −1 and 45 to 59 km s −1 are, respectively, shown in the left and right panel of Figure 3.
The intensity maps of the gamma-ray emission correlate more strongly with the gas distribution in the velocity range of 60 to 70 km s −1 . Particularly, the two point-like 4FGL sources (4FGL J1900.4+0339 and 4FGL J1858.0+0354) and Src A each coincide with one of the densest regions of the gas distribution. On the other hand, γ-ray emission inten-sity correlates poorly with gas distribution in the 45-59 km s −1 velocity range compared to that in the 60-70 km s −1 velocity range. This indicates that the extended γ-ray source is more likely to be related to Masgomas-6a. Miville-Deschênes et al. (2017) produced a catalog of 8107 molecular clouds that covers the entire Galactic plane and includes 98% of the 12 CO emission observed within b ± 5 • , using a hierarchical cluster identification method applied to the result of a Gaussian decomposition of CO data. We note that a molecular cloud in this catalog (named "[MML2017] 1165") 2 is located in the velocity range of V LSR = (63.52 ± 6.41) km s −1 , consistent with our measurement. The angular radius of this cloud is 0.518 • (see the magenta circle in the left panel of Figure 3). Interestingly, the position and size of this MC agree well with those of the extended GeV source, supporting the association between the GeV source and the MC.
The spectra of 12 CO and 13 CO emission toward the two 4FGL sources, Masgomas-6a and Src A are shown in Figure 4. We use a Gaussian function to fit the spectra of 12 CO emission around 65 km s −1 toward the two 4FGL sources, Src A and Masgomas-6a. The line center and the full width at half maximum (FWHM) in the four regions are, respectively, v c = 64.06 km s −1 with FWHM = 5.17 km s −1 for 4FGL J1858.0+0354, v c = 66.03 km s −1 with FWHM = 8.63 km s −1 for 4FGL J1900.4+0339, v c = 65.25 km s −1 with FWHM = 6.65 km s −1 for Masgomas-6a and v c = 63.70 km s −1 with FWHM = 2.42 km s −1 for Src A. Thus, one peak at ∼ 64 km s −1 is seen in the direction of 4FGL J1858.0+0354 and Src A, and another peak at ∼ 65 − 66 km s −1 is seen in the direction of 4FGL J1900.4+0339 and Masgomas-6a. This strengthens our confidence on the aforementioned correlation between the intensity of γ-ray emission and the gas distribution in the velocity range of 60 to 70 km s −1 . In addition, We do not find any significant evidence of broadenings or asymmetries in 12 CO line with respect to the narrow 13 CO line, indicating no shock interaction. Figure 5 displays the 12 CO(J = 1-0) and 13 CO(J = 1-0) maps for five consecutive velocity ranges from 60 to 70 km s −1 , with a step of 2 km s −1 for each map. We find that the molecular gas are observed to be spatially associated with the GeV emission in all the velocity range. In particular, the correlation is the best in the velocity ranges of 62 − 64 km s −1 and 64 − 66 km s −1 .
Adopting the mean CO-to-H 2 mass conversion factor X CO = 2 × 10 20 cm −2 K −1 km −1 s (Bolatto et al. 2013 we estimate that the total mass of gas within the 0.52 • disk of the GeV source is about 3.13 × 10 5 d 2 3.9 M ⊙ . If we assume the γ-ray emission region is spherical in geometry, then the radius is estimated to be R ∼ 36d 3.9 pc. Then the average atom density of the H 2 cloud in this region is about n ≃ 70d −1 3.9 cm −3 .

INTERPRETATION OF THE GEV EMISSION
The correlation between the γ-ray intensity and the molecular gas density favors the hadronic origin of the γ-ray emission, where GeV γ-rays are produced by the pp interaction between the CRs and the gas. We fit the γ-ray spectra of the source with the hadronic model using the Markov Chain Monte Carlo fitting routines of Naima, a package for the calculation of nonthermal emission from relativistic particles (Zabalza 2015). First, we assume the parent proton spectrum is a single power-law given by dN/dE = A(E/E 0 ) −α . The derived parameters are α = 2.57 +0.03 −0.04 , and the total energy is W p = 1.81 × 10 49 erg for the protons above 1 GeV, assuming the gas density is 70 cm −3 . The best-fit result is shown by the red dashed line in Figure 2. Below 1 GeV, the hard spectrum is due to the π 0 bump. At higher energies (> 10 GeV), the theoretical flux exceeds the LAT data slightly (but still within the 2σ uncertainty). This could indicate a possible cutoff in the proton spectrum around 100 GeV . Then, we consider a cutoff power-law function for the parent proton spectrum, i.e. dN/dE = A(E/E 0 ) −α exp(−(E/E cut ) β ). The derived parameters are α = 2.18 ± 0.05 and W p = 1.34 × 10 49 erg for the protons above 1 GeV, if we take β = 1, d = 3.9 kpc, n = 70 cm −3 and E cut = 100 GeV. The fitting result is shown by the blue line in Figure 2.
The young massive star cluster candidate Masgomas-6a contains two Wolf-Rayet stars (WR122-11 and WR122-16) and several OB-type stars. WR122-11 is classified as WN6 and WR122-16 is classified as WN8-9h (Ramírez Alegría et al. 2018). These WR stars have a mass loss rate 2 − 3 × 10 −5 M ⊙ yr −1 and the wind velocity is 1000 − 2000 km s −1 , so the kinetic luminosity of the stellar wind is about 10 37 erg s −1 , which is sufficient to power the γ-ray emission. As the mass-loss of WR stars are quasi-continuous over the lifetime T , one could expect a quasi-continuous injection of CRs into the interstellar medium over the age of WR stars (T 10 6 yr). The stellar wind from Masgomas-6a may drive a termination shock by interacting with the surrounding gas. The termination shock then accelerates cosmic rays, which produce gamma-ray emission by colliding with the molecular cloud.
The efficiency of production of γ-rays in the cloud is determined by the ratio between the diffusion timescale and the energy loss timescale of cosmic rays, i.e., η pp = t diff /t pp . The energy loss time of cosmic ray protons colliding with the molecular gas is t pp = 1/(σ pp )K pp nc = 10 6 yr(n/70cm −3 ) −1 , where σ pp = 3 × 10 −26 cm 2 is the cross section of pp collisions, K pp = 0.5 is the inelasticity of pion production and n is the gas number density. The diffusion time depends on the diffusion coefficient in the cloud, i.e. t diff = R 2 /4D, where R is the size of the source and D is the diffusion coefficient. The diffusion coefficient for CRs propagating in the molecular cloud is poorly known. In the interstellar medium, the diffusion coefficient is of order of 10 29 cm 2 s −1 for 10-100 GeV protons (see e.g. (Strong et al. 2007)) responsible for 1-10 GeV γ-rays. We assume a homogeneous diffusion coefficient and parameterize it as D = χD ISM = χ10 29 cm 2 s −1 for CRs with energy of 10−100 GeV, with χ being the ratio between the CR diffusion coefficient in the cloud and the average one of the Galactic ISM. The diffusion time is then given by t diff = 10 3 χ −1 (D ISM /10 29 cm 2 s −1 ) −1 yr for R = 36 pc. Thus, the required CR injection rate to power the γ-ray emission is L p = L γ /η pp = 1.8 × 10 38 χ erg s −1 . To account for the CR injection rate with the stellar wind power of WRs, χ 0.1 is required for a reasonable CR acceleration efficiency, implying that the diffusion coefficient is significantly suppressed in the molecular cloud. We also note that a significant part of CRs could have already escaped out of the γ-ray emission region. The relation between the total injected CR energy and the CR energy remained in the γ-ray emission region can be expressed as L p,tot /L p = (r diff /R) 2 = 100(χ/0.1)(T /10 6 yr), where r diff is the the total diffusion length of CRs corresponding to the age T of the WR stars.
Supernova remnants (SNRs) could also be a potential accelerator of cosmic rays in this region, since the energy release in CRs from a single SNR is 10 49 − 10 50 erg.
The CRs escaping from SNRs can produce γ-rays by interacting with the surrounding gas (Aharonian & Atoyan 1996;Rodriguez Marrero et al. 2008;Aharonian et al. 2008;Gabici et al. 2009;Uchiyama et al. 2012;Zhang et al. 2021). No recorded SNRs is found in this region and future multi-wavelength observations may be helpful to determine this scenario. The radial distribution of γrays can provide key information for the injection history of CRs, which is useful to distinguish between the continuous stellar wind scenario and the single SNR scenario (Aharonian et al. 2019). Pulsar wind nebulae (PWN) are another potential sources for the γ-ray emission. There is a middle-aged pulsar PSR J1858+0346 located positionally close to the center of the γ-ray sources, as marked by the grey dot in Figure 1, with a characteristic age of 2 Myr and a distance of 5.5 kpc (Manchester et al. 2005). If the γ-ray source is related to the PWN of the pulsar, the γ-ray luminosity would be about 3.5 × 10 35 erg/s, and cannot be explained based on the current spin-down power of PSR J1858+0346, which is 4.8×10 33 erg/s. However, we note that the cooling timescale of GeV-emitting electrons, t c ≈ 5(E e /100GeV) −1 (B/5µG) −2 Myr, is comparable to or even longer than the age of the pulsar assuming a typical interstellar magnetic field strength. As a result, considering a much higher spin-down power of the pulsar at its early age, we cannot exclude the pulsar/PWN origin of the γ-ray source.

SUMMARY
Massive star clusters hosting WR and OB-type stars have been suggested to be accelerators of Galactic CRs. We analyzed the Fermi Large Area Telescope's data towards the direction of Masgomas-6a, a young massive star cluster candidate hosting two WR stars, which is located at a distance of 3.9 kpc. An extended γ-ray source with a radius of 0.52 degree is found in the vicinity of Masgomas-6a.
Our analysis shows a correlation between the intensity of γ-ray emission and the density of molecular gas in the region. Combined with the fact that the distance of the molecular cloud inferred from the velocity distribution is consistent with that of Masgomas-6a, the correlation indicates that the γ-ray emission is related with Masgomas-6a. A natural scenario is that the two WR stars in Masgomas-6a accelerate cosmic rays, which produce the γ-ray emission through inelastic proton-proton collisions with the molecular gas. With a kinetic luminosity of ∼ 10 37 erg s −1 in the stellar winds, the two WR stars are capable of powering the γ-ray emission. The radius of the GeV source is about 36 pc for a source distance of 3.9 kpc. Assuming a normal CR diffusion coefficient, the vast majority of CRs accelerated by WR stars in their lifetime must have escaped out of the source as the injection of CRs is expected to last million years. The energy budget constraint suggests a CR diffusion coefficient smaller than that in the Galactic ISM, indicating a strong suppression of CR diffusion when they are propagating in the molecular cloud. The spectral data shown in Figure 2 seems to suggest the presence of a cutoff in the proton spectrum around 100 GeV, but we note that the errors of the data at high energies are rather large, so one can not exclude a single power-law spectrum. Thus, the maximum energy of cosmic rays accelerated by Masgomas-6a should be least ∼ 100 GeV.

ACKNOWLEDGMENTS
We would like to thank the referee for the constructive report. This work was supported by the National KeyR & D program of China under the grant 2018YFA0404203, the NSFC grants 12121003 and U2031105, and China Manned Spaced Project (CMS-CSST-2021-B11). This research made use of the data from the Milky Way Imaging Scroll Painting (MWISP) project, which is a multi-line survey in 12 CO/ 13 CO/C 18 O along the northern galactic plane with PMO-13.7m telescope. Thanks to NASA Fermi Collaboration for providing the excellent data and tools that made this work possible. We are grateful to all the members of the MWISP working group, particularly the staff members at     Note: PS, Disk and Gaussian represent, respectively, a point-like source model, an uniform disk model and a two-dimensional Gaussian model. PS+PS represents the two point source model with source positions at the two 4FGL sources. Extension in the disk and Gaussian models represents the radius containing 68% of the intensity of the tested models. N dof represents the number of degrees of freedom for each model. 1-500 GeV. The red dash line represents the modeling assuming a single power-law proton spectrum. Model parameters are α = 2.57 +0.03 −0.04 (the proton spectral index) and Wp = 1.81 × 10 49 erg (the total energy in protons above 1 GeV). The blue solid line represents the modelling assuming a power-law with a high-energy cutoff for the proton spectrum. Model parameters are α = 2.18 ± 0.05, β = 1, Ecut = 100 GeV, and Wp = 1.34 × 10 49 erg for the protons above 1 GeV. The distance of the source is taken to be d = 3.9 kpc and the gas number density is taken to be n = 70 cm −3 .      Figure 5. 12 CO(J=1-0; left column) and 13 CO(J=1-0; right column) intensity map for five consecutive velocity ranges from 60 to 70 km s −1 , with a step of 2 km s −1 for each map. White contours correspond to Fermi-LAT significance map starting from 3σ to 8σ by 1σ steps.