Investigating the Nature of MGRO J1908+06 with Multiwavelength Observations

MGRO J1908+06 is an extended bright TeV source of unknown nature. It was first discovered by the Milagro water Cherenkov telescope in its sky survey results after seven years of operation (Abdo et al. 2007). MGRO J1908+06 was subsequently detected in the TeV range by the High Energy Stereoscopic System (H.E.S.S.; Aharonian et al. 2009), the Very Energetic Radiation Imaging Telescope Array System (VERITAS; Ward 2008; Aliu et al. 2014), the Astrophysical Radiation with Ground-based Observatory at YangBaJing (ARGO-YBJ; Bartoli et al. 2012), and recently by the High-Altitude Water Cherenkov Observatory (HAWC; Abeysekara et al. 2017). The TeV luminosity of MGRO J1908+06 is comparable to the Crab Nebula (Bartoli et al. 2012), making it one of the most luminous Galactic gamma-ray sources in the TeV range. MGRO J1908+06 is among the four sources detected above 100 TeV in a recent study by HAWC, indicating its possible ability to accelerate particles to PeV energy (Abeysekara et al. 2020). In addition to its possible PeVatron nature, MGRO J1908+06 is spatially associated with an IceCube neutrino emission hotspot (Aartsen et al. 2019, 2020), though the post-trial significance is low.

The nature of MGRO J1908+06 remains unknown. Searches for its multiwavelength counterparts in radio, X-ray, and GeV gamma-rays have been unsuccessful (Kong 2007; Abdo et al. 2010a; Pandel 2015; Duvidovich et al. 2020). MGRO J1908+06 is spatially associated with a middle-aged (20–40 kyr, Downes et al. 1980) supernova remnant (SNR) G40.5−0.5 and an energetic gamma-ray pulsar PSR J1907+0602, with a spin period of 107 ms, a spin-down luminosity of ∼2.8 × 10³⁶ erg s⁻¹, and a characteristic age of 19.5 kyr. The distance to PSR J1907+0602, estimated from its dispersion measure, is 3.2 ± 0.6 kpc (Abdo et al. 2010a), and distance estimates of SNR G40.5−0.5 place it either in a similar neighborhood (3.5 kpc, using CO observations; Yang et al. 2006) or a more distant region (using the Σ-D relation, 5.5–8.5 kpc; Downes et al. 1980; 6.1 kpc, Case & Bhattacharya 1998). PSR J1907+0602 was not likely born in SNR G40.5−0.5, considering the high estimated transverse velocity, and the lack of bow shock and trail morphology. The recent discovery of the radio pulsar PSR J1907+0631 (7.9 kpc; Lyne et al. 2017) in the projected center of the remnant promotes a larger distance estimation for SNR G40.5−0.5. ¹¹ In addition to PSR J1907+0602, GeV source 4FGL J1906.2+0631 is also spatially associated with MGRO J1908+06 (Abdollahi et al. 2020). Below we summarize these multiwavelength spatial associations.

• 1.  
  PSR J1907+0602, GeV and radio pulsar, distance ∼3.2 kpc.
• 2.  
  4FGL J1906.2+0631, unidentified GeV source, distance unknown.
• 3.  
  SNR G40.5−0.5, radio SNR, likely distance ∼8 kpc.
• 4.  
  PSR J1907+0631, radio pulsar, likely associated with SNR G40.5−0.5, distance ∼8 kpc.

The nature of MGRO J1908+06 is under debate. A leptonic pulsar wind nebula (PWN) scenario powered by PSR J1907+0602 and a hadronic scenario powered by SNR G40.5-0.5 have been proposed (Abdo et al. 2010a). Duvidovich et al. (2020) proposed that the TeV emission of MGRO J1908+06 originates from a combination of the two above scenarios, which play their own roles at different distances. To gain further insight into these scenarios, we have analyzed multiwavelength observations in the vicinity of MGRO J1908+06 and modeled the emission processes, and we report our results below.

Multiwavelength observations were investigated in this Letter, including CO observations from the Milky Way Imaging Scroll Painting (MWISP) project (Su et al. 2019), Fermi Large Area Telescope (LAT) observations of GeV gamma-rays, XMM-Newton observations in X-ray, and the Very Large Array Galactic Plane Survey (VGPS; Stil et al. 2006) in radio. The details of these observations are reported in the Appendices.

3.1. Molecular Clouds in the Region of MGRO J1908+06

Molecular clouds (MCs) toward the SNR G40.5−0.5 region have been studied in the past (Yang et al. 2006; Duvidovich et al. 2020). Using data from the MWISP project (Su et al. 2019), we searched for MCs with ¹²CO (J = 1–0),¹³CO(J = 1–0), and C¹⁸O (J = 1−0) emission lines in a larger region covering MGRO J1908+06. MCs are discovered to be spatially associated with MGRO J1908+06 in the ¹²CO (J = 1–0) and ¹³CO (J = 1–0) maps between 46 and 66 km s⁻¹ (Figure 1; see also Appendix A for details), which is consistent with the possible distance, 3.2 kpc, to SNR G40.5−0.5 and PSR J1907+0602. A shell-like cavity encircling the radio morphology of SNR G40.5−0.5 is observed in both ¹²CO (J = 1–0) and ¹³CO (J = 1–0) maps, indicating a possible SNR swept-up shell (Figure 1).

Zoom In Zoom Out Reset image size

We show the ¹²CO (J = 1–0) and ¹³CO (J = 1–0) maps for five consecutive velocity ranges from 46 km s⁻¹ to 66 km s⁻¹ with a coverage of 4 km s⁻¹ in Appendix A, Figure 5. Different regions with MCs are studied in detail and their astrophysical properties are estimated in Appendix A. The mean density is estimated to be ∼45 cm⁻³. In the 46 km s⁻¹ to 66 km s⁻¹ velocity slices, there are apparent filaments of MCs positionally located next to PSR J1907+0602. We have inspected the ¹²CO (J = 1−0) and ¹³CO (J = 1−0) line profiles of the MCs toward the SNR G40.5−0.5 and PSR J1907+0602 region, to search for kinematic evidence for gas distribution due to external interaction (e.g., Zhou & Chen 2011; Liu et al. 2020). However, we do not find any such evidence of asymmetric broad profiles of the ¹²CO line (i.e., a wing part deviating from the main Gaussian in a ¹²CO profile where the signal-noise-ratio of the ¹³CO is less than 3) in the grid of the CO spectra. For a further search of interaction signals, we studied the data from the INT Photometric Hα Survey of the Northern Galactic Plane (IPHAS) covering the MGRO J1908+06 region. No Hα ascribable to interaction was detected.

3.2. GeV Counterpart of MGRO J1908+06

Though being very bright at TeV energies, no GeV counterpart for MGRO J1908+06 was identified in previous Fermi-LAT studies (Abdo et al. 2010a; Ackermann et al. 2011). Two gamma-ray sources PSR J1907+0602 and 4FGL J1906.2+0631 are spatially associated with MGRO J1908+06, as listed in the 4FGL catalog (Abdollahi et al. 2020). For a deeper search of its GeV counterpart, we analyzed more than 11 yr of Fermi-LAT data in the 0.1–300 GeV band. To minimize the contamination from the bright gamma-ray pulsar PSR J1907+0602, we carried out data analysis during the off-peak phases of this pulsar. A full description of our data analysis is presented in Appendix B. We discovered a previously undetected, extended GeV source spatially associated with MGRO J1908+06, hereafter referred as Fermi J1906+0626, which is shown in Figure 2 together with the TeV extension measured by H.E.S.S. and the TeV contours reported by VERITAS. The good morphological agreement between Fermi J1906+0626 and MGRO J1908+06 strongly suggests a common origin. To investigate this, we derived the morphological and spectral properties of Fermi J1906+0626.

Zoom In Zoom Out Reset image size

A morphology study was carried out in 1–300 GeV to estimate the extension of Fermi J1906+0626 assuming a power-law spectral model. We tested both the disk and Gaussian morphologies using Fermipy, and the disk morphology yields more significant extension detection. The disk morphology resulted in a center position of R.A. = 28661 ± 008, decl. = 644 ± 007, and a radius of 083 ± 005, yielding a TS_ext = 44 (see Appendix B). Comparing the likelihood values, the disk morphology is significantly preferred over a two-point-source model (PSR J1907+0602 plus 4FGL J1906.2+0631 from 4FGL) with a ΔTS = 23 (see Appendix B). We further considered a two-component spatial model, constituted by a template using VERITAS counts map (Aliu et al. 2014) plus the point source 4FGL J1906.2+0631. However, comparing the likelihood values, the disk morphology is again significantly preferred over this two-component spatial model with a ΔTS = 18. Additionally, in this two-component spatial model 4FGL J1906.2+0631 is not significantly detected either in 1–300 GeV or 0.1–300 GeV, leading to a TS value of 11 and 20, respectively, which are lower that the detection threshold (TS = 25, 4FGL). Very recently, Di Mauro et al. (2020) carried out an analysis of MGRO J1908+06 and detected extended emission, which is consistent with the results shown here.

We carried out the spectral analysis in 0.1–300 GeV. Adopting the disk morphology, we tested a power law (${dN}/{dE}={N}_{0}{(E/{E}_{0})}^{-{\rm{\Gamma }}}$ cm⁻² s⁻¹ MeV⁻¹) and a log-parabola spectral model (${dN}/{dE}={N}_{0}{(E/{E}_{b})}^{-(\alpha +\beta \mathrm{log}(E/{E}_{b}))}$ cm⁻² s⁻¹ MeV⁻¹). The log-parabola spectral model is preferred with a ΔTS = 89, indicating a significant spectral curvature. Low-energy turnovers are expected for proton–proton interactions, indicating a hadronic nature for MGRO J1908+06. With the disk morphology and log-parabola spectral model, Fermi J1906+0626 is significantly detected with a TS value of 190 and an energy flux of (7.44 ± 0.64) ×10⁻¹¹ erg cm⁻² s⁻¹ between 0.1 and 300 GeV. ¹² Adopting the disk morphology, the best-fitted parameters are shown in Table 1 and the gamma-ray spectral energy distribution (SED) of Fermi J1906+0626 is shown in Figure 2, top-right panel (in black). An excess beyond ∼5 GeV, deviating from the log-parabola spectrum, suggests the existence of another spectral component, which may originate from a second population of accelerated particles, or from a second emitting region.

Table 1. Best-fit Spectra Parameters of Fermi J1906+0626

		Off-peak Analysis from 0.1 to 300 GeV
Model	Γ	α	β	Energy Flux	TS
	(Power law)	(Log-parabola, Eb = 1 GeV)	(Log-parabola, Eb = 1 GeV)	(10−11 erg cm−2 s−1)
Log-parabola	⋯	2.90 ± 0.20	0.65 ± 0.11	7.44 ± 0.64	190
(083 disk)
Log-parabola + power law (fixed)	1.61 (fixed)	3.16 ± 0.24	0.78 ± 0.13	9.34 ± 0.64	197
(083 disk) (051 disk)
		Analysis from 30 GeV–1 TeV
Model	Γ	α	β	Energy Flux	TS
	(Power law)	(Log-parabola, Eb = 1 GeV)	(Log-parabola, Eb = 1 GeV)	(10−11 erg cm−2 s−1)
Power law	1.61 ± 0.16	⋯	⋯	3.05 ± 0.70	47
(051 disk)

Download table as:  ASCIITypeset image

For higher statistics, we studied Fermi-LAT data above 30 GeV further without pulsar gating. With a spectral cutoff at 2.9 GeV, the magnetospheric emission from PSR J1907+0602 is negligible above 30 GeV (Abdo et al. 2013), making pulsar gating unnecessary. In the 30 GeV to 1 TeV energy range Fermi J1906+0626 is detected as an extended source (Figure 2, bottom-right panel). The best-fitted disk model is centered at R.A. = 28688 ± 005, decl. = 629 ± 005 with a radius of 051 ± 002, yielding a TS of 47 assuming a power-law spectral model and a TS_ext = 46. Considering the possible spectral connection with the TeV range, we tested the template produced from VERITAS counts map (Aliu et al. 2014). Assuming a power-law model, Fermi J1906+0626 is detected from 30 GeV to 1 TeV with a TS of 43, which is a worse fit than a disk model. Adopting the disk model, Fermi J1906+0626 shows a spectral index of 1.61 ± 0.16 and an energy flux of (3.05 ± 0.70) ×10⁻¹¹ erg cm⁻² s⁻¹ (Table 1), confirming the existence of an additional spectral component at higher energies beyond the log-parabola component (Figure 2, top-right panel). We added this component to the analysis of Fermi J1906+0626 from 0.1 to 300 GeV in off-peak data, fixed it and constructed a two-component spectral model (log-parabola plus power law). The results are reported in Table 1 and shown in Figure 2, top-right panel. The multiwavelength SED of MGRO J1908+06 is shown in Figure 3. The GeV SED shows a natural continuity to the TeV range. The TS map in the 30 GeV to 1 TeV energy range is shown in Figure 2, bottom-right panel. No clear ¹²CO (J = 1–0) emission is associated with gamma-ray emission.

Zoom In Zoom Out Reset image size

Below 2 GeV, the SED of Fermi J1906+0626 is dominated by the log-parabola component. Assuming a disk morphology (Figure 2, bottom-left panel), the best-fitted values are R.A. = 28671 ± 0.07, decl. = 630 ± 0.07 with a radius of 068 ± 0.07, yielding a TS value of 109. The peak of gamma-ray emission is spatially consistent with the ¹²CO (J = 1–0) intensity in the 62–66 km s⁻¹ range, suggesting a hadronic origin. We produced a template from the 62–66 km s^{−1 12}CO (J = 1–0) intensity map within the best-fitted disk morphology. The template leads to a better fitting with TS value of 113, which is preferred over the disk morphology but not by a large margin.

3.3. Search for MGRO J1908+06 Counterparts in X-Ray and Radio Wavelengths

The XMM-Newton X-ray satellite has covered MGRO J1908+06 with five observations (Obs. ID 0553640101, 0553640201, 0553640701, 0553640801, and 0605700201), providing a total exposure of 109 ks. Combining all available MOS data, we produced a particle background-subtracted and exposure-corrected count rate map for the MGRO J1908+06 region (Figure 4, left panel). No diffuse X-ray emission coincident with the gamma-ray emission region of MGRO J1908+06 could be detected in the 0.2–10 keV energy range. The morphology measured by H.E.S.S. (a Gaussian profile with σ of 034; Aharonian et al. 2009) is the only one in TeV range whose 1σ extension is fully covered by these XMM-Newton observations. Considering the 1σ range of H.E.S.S. morphology, excluding point sources within, assuming a spectral index of 2 and a H i column density of N_H = 1.3 × 10²² cm⁻² (Abdo et al. 2010a), we calculated a 95% unabsorbed upper limit for MGRO J1908+06 as 1.2 ×10⁻¹⁰ erg cm⁻² s⁻¹ (0.2−10 keV). From our CO observation, the local H i column density to the MGRO J1908+06 region is N_H = 2 × N ${}_{{{\rm{H}}}_{2}}\,\sim (2\mbox{--}7)\times {10}^{21}$ cm⁻² (Appendix A, Table 2), which is much lower than the total Galactic H i column density in this direction estimated using the HEASARC tool ¹³ and with Chandra (N_H = 1.3 × 10²² cm⁻²; Abdo et al. 2010a), suggesting that the local absorption does not likely lead to the lack of an X-ray counterpart.

Zoom In Zoom Out Reset image size

Table 2. Fitted and Derived Parameters for the MCs Around 50 km s⁻¹ in 4 Regions as Indicated in Appendix A, Figure 6

Region	N(H2)	n(H2)	M(H2)	FWHM	Line Center	Near/Far Distance
	(1021 cm−2)	(cm −3)	(103 M⊙)	( km s−1)	( km s−1)	(kpc)
1	1.2	$16{d}_{7.9}^{-1}$	$55.6{d}_{7.9}^{2}$	2.7	52.9	3.6/9.3
2	3.2	$16{d}_{7.9}^{-1}$	$166.9{d}_{7.9}^{2}$	3.8	57.9	4.0/8.9
3	1.0	$180{d}_{7.9}^{-1}$	$23.1{d}_{7.9}^{2}$	1.7	50.1	3.4/9.5
4	3.4	$24{d}_{7.9}^{-1}$	$125.0{d}_{7.9}^{2}$	3.5	56.9	3.9/9.0

Download table as:  ASCIITypeset image

With a 19 ks Chandra observation (Obs. ID 7049 on 2009 August 19), Abdo et al. (2010a) reported the detection of PSR J1907+0602 and hints of extension from 2−8 keV. However, because of the low statistics, no PWN associated with PSR J1907+0602 could be unambiguously identified. In the XMM-Newton data, Obs. ID 0605700201 (on 2010 April 26, 52.6 ks exposure), ID 0553640201 (on 2008 October 2, 24.7 ks exposure), ID 0553640701 (on 2009 March 19, 7.4 ks exposure) have covered PSR J1907+0602. However, PSR J1907+0602 is only detected in Obs. ID 0605700201, which provided the longest exposure (52.5 ks), as a point source, by combining MOS 1 and MOS 2 data. To check the previously claimed possible extension by Chandra (Abdo et al. 2010a), we produced counts map from 2−8 keV of PSR J1907+0602, combining MOS 1 and MOS 2 data. All counts are well located within the radius of ∼90% fractional encircled energy ¹⁴ (Figure 4, top-right panel). The detected emission is thus consistent with a point source and no associated PWN could be identified.

We checked for variability with spectral analysis of the eight-month separated XMM-Newton (Obs. ID 0605700201 on 2010 April 26) and Chandra (Obs. ID 7049 on 2009 August 19) observations. We adopted a simple power-law model plus absorption with H i column density fixed at N_H = 1.3 × 10²² cm⁻² (Abdo et al. 2010a) because of the low overall counts. The 1−10 keV XMM-Newton MOS 1&2 data yield a best-fit spectral index of ${1.79}_{-0.41}^{+0.42}$ with unabsorbed total energy flux of ${6.29}_{-1.18}^{+1.44}$ × 10⁻¹⁴ erg cm⁻² s⁻¹, whereas the Chandra ACIS spectrum is well fitted with compatible spectral parameters, a best-fit spectral index of ${1.28}_{-0.42}^{+0.45}$ and unabsorbed total energy flux of ${5.14}_{-0.94}^{+1.47}$ × 10⁻¹⁴ erg cm⁻² s⁻¹. No spectral or flux variability can be claimed. The XMM-Newton and Chandra spectra are shown in Appendix C.

During this XMM-Newton observation, the PN camera was operating in small window mode, providing sufficient time resolution (5.7 ms) to search for X-ray pulsations. We extracted photons from PN data using a radius of 20'' in 0.2−10 keV and 3−10 keV with barycenter correction applied. Using Tempo2 (Hobbs et al. 2006) and the photons plugin ¹⁵ and contemporaneous Fermi-LAT gamma-ray ephemeris, we assigned pulsar rotational phase to each extracted photon. No significant X-ray pulsation was detected.

The VGPS at 1420 MHz (Stil et al. 2006) has covered MGRO J1908+06 region and the image is shown in Figure 4, bottom panel. PSR J1907+0602 and its possible associated PWN are not detected, but SNR G40.5−0.5 is clearly visible (Abdo et al. 2010a). No diffuse radio large-scale emission associated with MGRO J1908+06 is seen.

Detailed analysis of the Fermi-LAT data revealed that the gamma-rays from the direction of MGRO J1908+06 follow a two-component spectrum. Based on the TS maps and multiwavelength observations, we propose two accelerators on the field: one region related to SNR G40.5−0.5, and a second one on the direction of the gamma-ray pulsar PSR J1907+0602. The high-energy component discovered in Fermi-LAT data shows a hard spectrum of α_he ∼ 1.6 and connects smoothly with the very-high-energy spectrum measured by the Cherenkov instruments (Figure 3; Aharonian et al. 2009; Aliu et al. 2014; Abeysekara et al. 2020). Previous studies (Abdo et al. 2010a; Abeysekara et al. 2017) have tentatively associated MGRO J1908+06 with a PWN powered by PSR J1907+0602 at a distance of 3.2 kpc. The high-energy (>30 GeV) component measured by Fermi-LAT, combined with the very-high energy spectrum measured by the Cherenkov instruments, do indeed resemble the spectral signature associated with inverse Compton emission from GeV/TeV PWNe (e.g., Crab, Abdo et al. 2010b; MSH 15-52, Abdo et al. 2010c; HESS J1825-137, Grondin et al. 2011). The low-energy component (below a few GeV) is described by a soft spectrum, similar to the ones observed on evolved SNRs (Acero et al. 2016), ¹⁶ and it shows a significant peak coincident with an enhancement of molecular material (see Figure 2, bottom-left panel, and results shown in the Appendices), implying a tentative hadronic origin. Within 50 pc region of SNR G40.5−0.5 and PSR J1907+0631, there are located several MCs that have compatible distances (Table 2, far distance). Considering that PSR J1907+0631 is not a gamma-ray pulsar, we attempt to associate the emission to the CR-MC interaction for the low-energy contribution at ∼8 kpc, whereas the high-energy part of the spectrum is associated to leptonic emission associated to the putative PWN of PSR J1907+0602 at ∼3.2 kpc.

To model the total emission in the region we assumed a combination of hadronic and leptonic emission, as is shown in Figure 3. For the hadronic pp collision, the gamma-ray flux is calculated using the parameterized formulae provided by Kamae et al. (2006), for a proton spectrum of the form of a power-law function in energy space, with a slope s_p . Assuming a distance of 8 kpc, the steep gamma-ray spectrum at low energy can be well reproduced by employing s_p = 2.8, which is consistent with the measured CR slope in the local interstellar medium. Such a steep slope may reflect a diffusion-modified spectrum, i.e., the injection spectrum is softened by the energy-dependent diffusion of CR protons with a diffusion coefficient D(E) ∝ E^δ . We do not further specify the value of the injection spectral slope of CR protons and the index of the diffusion coefficient δ, but just note that an injection spectral slope of 2.3 − 2.4 and δ = 0.4–0.5 could be a reasonable combination of these two parameters, where the former one is motivated by the discovery of the pionic gamma-ray components from two other SNRs, i.e., W44 and IC 443 (Ackermann et al. 2013), and the latter one is based on the observation and modeling of secondary-to-primary CR ratios (Aguilar et al. 2016; Génolini et al. 2019; Huang et al. 2020). We assume an average hydrogen density of n = 45 cm⁻³ in the surrounding MC, resulting in a pp collision cooling timescale of ${t}_{{pp}}\simeq {\left(0.5{\sigma }_{{pp}}{nc}\right)}^{-1}=7\times {10}^{5}{(n/100{\mathrm{cm}}^{-3})}^{-1}\,$ yr. We multiply a factor of 1.4 to the obtained spectra of secondaries to account for the contribution of nuclei in molecular materials. The total proton energy needed to account for the gamma-ray emission is W_p ≃ 2 × 10⁵⁰ erg, which is well consistent with the reach of the usual 10% of the kinetic energy released in SNRs (Aharonian et al. 2004). Note that if the distance of the SNR and MCs is 3.2 kpc, the inferred gas density of the MC would be 2.5 times higher (see Appendix A, Figure 6 and Table 2), and this change reduces the requirement for the proton energy to 10⁴⁹ ergs.

For the leptonic component, we used an electron/positron broken-power-law distribution, i.e., ${dN}/{{dE}}_{e}\propto {E}_{e}^{-{s}_{e,1}}$ for E_e < E_b and ${dN}/{{dE}}_{e}\propto {E}_{e}^{-{s}_{e,2}}$ for E_e ≥ E_b , which is usually chosen to describe the SED of PWNe (Tanaka & Takahara 2010; Bucciantini et al. 2011; Martín et al. 2012; Torres et al. 2013, 2014). The time-independent spectra of synchrotron radiation and IC radiation are calculated following Blumenthal & Gould (1970). The IC emissivity is calculated in the optically thin case, which is appropriate for these objects, using the general Klein–Nishina differential cross section. We adopt the interstellar radiation field modeled in Popescu et al. (2017) as well as the cosmic microwave background as the target photon field for the IC scattering. The Fermi data above 10 GeV together with the HAWC data can be reproduced with s_e,1 = 1.5, s_e,2 = 3.0 and a break energy at E_b = 8.2 TeV, as well as a total energy of W_e ≃ 4 × 10⁴⁷ erg in the emitting electrons/positrons. ¹⁷ Assuming the age of the system equal to the characteristic age of PSR J1907+0602, i.e., 20 kyr, we can roughly evaluate the average magnetic field of a PWN from such a large breaking energy, via equating the age of the system to the cooling timescale of the electrons at the breaking energy ${t}_{\mathrm{cool}}\simeq 40{(E/8.2\mathrm{TeV})}^{-1}{\left[{({{\rm{U}}}_{\mathrm{ph}}+{{\rm{U}}}_{B})}^{-1}/1\mathrm{eV}{\mathrm{cm}}^{-3}\right]}^{-1}$ kyr. Given U_ph = 1 eVcm⁻³ and U_B = B²/8π, the inferred magnetic field strength is B ≃ 6 μ G. Such a magnetic field is consistent with the X-ray upper limit posed by XMM-Newton. The low magnetic field strength is similar to some other relic nebulae in the TeV regime (Aharonian et al. 2006; H.E.S.S. Collaboration et al. 2012; Liu et al. 2019) associated to intermediate-aged pulsars. Some of those PWNe display energy-dependent morphology in the TeV regime (Aharonian et al. 2006; H.E.S.S. Collaboration et al. 2012, 2019). Above 30 GeV, the best fit to the LAT data is the template adopted from VERITAS data, consistent with the morphology obtained in TeV range. If the emission is indeed related to the pulsar, we expect the nebula to be more compact and closer to the pulsar at TeV energies. The energy-dependent morphology could be the key to understanding the transport mechanism of particles within the PWN and the evolution of the PWN (H.E.S.S. Collaboration et al. 2019; Liu & Yan 2020). Deep observations with TeV observatories such as H.E.S.S., HAWC, or LHAASO will provide crucial input to disentangle the origin of the gamma-ray emission observed.

We also note that there are actually many relevant physical processes that can influence the modeling, such as the particle injection history, particle spectral evolution, and particle transport. We leave such a more realistic modeling to the future study and here we simply test the feasibility of the hybrid interpretation. In the considered hybrid scenario, the neutrino emission from MGRO J1908+06 would not be detectable by current instruments that are operating above 100 GeV (e.g., IceCube). This is because the neutrino spectrum arising from pp collisions generally resemble that of the pionic gamma-ray spectrum, which is important only below ∼10 GeV and drops quickly with energy. However, we should also note that it is not currently clear whether the gamma-ray spectrum above 100 TeV would decline as in our model expectation. IceCube Collaboration (Aartsen et al. 2020) found that the best-fit ${\nu }_{\mu }+{\bar{\nu }}_{\mu }$ number in this region is 4.2 and the best-fit spectral slope is −2. This is translated to a 90% C.L. upper limit for the neutrino flux as ${dN}/{{dE}}_{\nu }\,=5.7\,\times {10}^{-13}{({E}_{\nu }/1\mathrm{TeV})}^{-2}\,{\mathrm{TeV}}^{-1}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$. In the pp collision, the pionic gamma-ray flux is about twice that of the ${\nu }_{\mu }+{\bar{\nu }}_{\mu }$ flux. Therefore, an upper limit for the accompanying pionic gamma-ray component can be obtained as shown in the dashed gray line in Figure 3. If a hardening of the gamma-ray spectrum beyond 100 TeV presents in the future observation by HAWC or LHAASO, a second hadronic gamma-ray component, as well as neutrino detection, may then be confirmed.

The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration (NASA) and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.

J.L. acknowledges the support from the Alexander von Humboldt Foundation and the National Natural Science Foundation of China via NSFC-11733009. R.-Y.L. acknowledges the support from the National Natural Science Foundation of China via NSFC-U2031105. E.O.W. acknowledges the support from the Alexander von Humboldt Foundation. The work of D.F.T. has been supported by the grants PGC2018-095512-B-I00, SGR2017-1383, and AYA2017-92402-EXP. Q.C.L. acknowledges support from the program A for Outstanding PhD candidate of Nanjing University. Work at NRL is supported by NASA. This research made use of the data from the Milky Way Imaging Scroll Painting (MWISP) project, which is a multi-line survey in 12CO/13CO/C18O along the northern galactic plane with PMO-13.7m telescope. We are grateful to all the members of the MWISP working group, particularly the staff members at PMO-13.7m telescope, for their long-term support. MWISP was sponsored by National Key R&D Program of China with grant 2017YFA0402701 and CAS Key Research Program of Frontier Sciences with grant QYZDJ-SSW-SLH047.

The CO data used in this work are part of the MWISP project (Su et al. 2019), including ¹²CO (J = 1–0), ¹³CO (J = 1–0), and C¹⁸O (J = 1–0), which were observed simultaneously using the 13.7 m millimeter-wavelength telescope of the Purple Mountain Observatory at Delingha. A new 3 × 3 pixel Superconducting Spectroscopic Array Receiver was used as the front end (Shan et al. 2012). The bandwidth of the spectrometer is 1 GHz and the half-power beamwidth of the telescope is about $50^{\prime\prime}$ at 115 GHz. The spectral resolution is 61 KHz, corresponding to velocity resolutions of 0.16 km s⁻¹ for ¹²CO and 0.17 km s⁻¹ for ¹³CO and C¹⁸O. The typical rms noise level is about 0.5 K for ¹²CO (J = 1–0) and 0.3 K for ¹³CO (J = 1–0) and C¹⁸O (J = 1–0).

MCs are observed to be spatially associated with the TeV emission (Figure 1). Figure 5 displays the ¹²CO (J = 1–0) and ¹³CO (J = 1–0) maps for five consecutive velocity ranges from 46 to 66 km s⁻¹, with a coverage of 4 km s⁻¹.

Zoom In Zoom Out Reset image size

We made an estimation of the astrophysical properties of the MCs in four regions and have parameterized the distance as d = 7.9d_7.9 kpc. We estimated the kinematic distance to the MCs using the Milky Way's rotation curve suggested by Brand & Blitz (1993), assuming the Sun's Galactocentric distance to be 8.5 kpc and orbital speed to be 220 km s⁻¹. Therefore, the velocity of each MC could indicate two candidate kinematic distances, the near side one and the far side one. ¹²CO(J = 1–0) and ¹³CO(J = 1–0) spectra of the four regions are extracted and shown in Figure 6. The estimated results are summarized in Table 2. Among the parameters, the column densities are estimated by using ¹³CO lines, under the assumption of local-thermal-equilibrium condition, optically thin conditions for ¹³CO (J = 1–0) line, and optically thick conditions for ¹²CO (J = 1–0). The excitation temperature is assumed to be T_ex = 12 K, 12 K, 21 K, and 15 K, respectively, the value estimated from the maximum ¹²CO (J = 1–0) line emission. Also, the conversion relation for the molecular column density of N(H₂) ≈ 7 × 10⁵ N(¹³CO) (Frerking et al. 1982) has been used. The mean density of each region is calculated by dividing the column density toward the ¹³CO emission peaks by the cloud size along the line of sight, which is estimated from the FWHM of the ¹³CO line with Larson's law (Larson 1981). Furthermore, a mean density of the four regions, when taking the weight of the volume into consideration, is estimated to be ∼45 cm⁻³.

Zoom In Zoom Out Reset image size

The analysis shown in this Letter uses more than 11 years of Fermi-LAT P8R3 data, from 2008 August 4 (MJD 54682) to 2019 November 9 (MJD 58796). All gamma-ray photons within a circular region of interest (ROI) of 15° radius centered on PSR J1907+0602 are considered in our analysis. The gamma-ray events from "Pass 8" event reconstruction are analyzed using the Fermi Science Tools 11-07-00 release. In the data reduction, a zenith angle threshold of 90° is adopted to reject contamination from gamma-rays from the Earth's limb. The selected Fermi-LAT instrument response functions (IRFs) is "P8R3 V2 Source." Known gamma-ray sources from the Fermi Large Area Telescope Fourth Source Catalog (4FGL; Abdollahi et al. 2020) within 20° of PSR J1907+0602 were included in the spectral-spatial model, as well as Galactic ("gll_iem_v07.fits") and isotropic diffuse emission components ("iso_P8R3_SOURCE_V2_v1.txt"). The spectral parameters of the sources within 4° of PSR J1907+0602, Galactic and isotropic diffuse emission components were all left free. PSR J1907+0602 and 4FGL J1906.2+0631 are spatially associated with MGRO J1908+06, thus not included in the spectral-spatial model. The spectral parameters of sources with larger angular separations were fixed at the 4FGL values. The spectral analysis was performed using a binned maximum likelihood fit (spatial bin size 01 and 30 logarithmically spaced bins in the 0.1–300 GeV range) For the analysis in 30 GeV–1 TeV, the spectral parameters of the sources within 4° of PSR J1907+0602 and the Galactic diffuse emission component are fixed to 4FGL values except for the prefactor (spectral normalization) because of low statistics.

The significance of the sources was evaluated by the Test Statistic (TS). This statistic is defined as TS =$-2\mathrm{ln}({L}_{\max ,0}/{L}_{\max ,1})$, where ${L}_{\max ,0}$ is the maximum likelihood value for a model in which the source studied is removed (the "null hypothesis"), and ${L}_{\max ,1}$ is the corresponding maximum likelihood value with this source being incorporated. The square root of the TS is approximately equal to the detection significance of a given source. The significance of source extension was defined as TS_ext = $-2\mathrm{ln}({L}_{\mathrm{point}}/{L}_{\mathrm{ext}}$), where L_ext and L_point are the gtlike global likelihood of the extended source hypotheses and the point source, respectively. The threshold for claiming the source to be spatially extended is set as TS_ext > 16, which corresponds to a significance of ∼4σ. The source localization, extension fitting, and TS maps production were carried out using the Fermipy analysis package (version 0.17.4; Wood et al. 2017). Energy dispersion correction has been applied in the analysis. The SEDs are computed assuming a power-law shape with spectral index fixed at 2.

PSR J1907+0602 is a bright gamma-ray pulsar spatially associated with MGRO J1908+06. To minimize contamination from this pulsar, we carried out data analysis during the off-peak phases of PSR J1907+0602. Using Tempo2 (Hobbs et al. 2006) with the Fermi plugin (Ray et al. 2011), we have assigned pulsar rotational phases for each gamma-ray photon that passed the selection criteria, adopting the most updated ephemeris for PSR J1907+0602. The pulse profile of PSR J1907+0602 is shown in Figure 7. We followed the off-peak definition in Li et al. (2020), which is ϕ = 0.0−0.136 and 0.697−1.0. To account for the off-peak phase selection, the prefactor parameter of all sources were scaled by 0.439, the width of the off-peak interval.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the off-peak analysis of the 0.1–300 GeV band, we searched for significant TS excess beyond Fermi J1906+0626 in the TS map within 4° of PSR J1907+0602. A new point source (hereafter PS 1) is located at R.A. = 28719 ± 006, decl. = 707 ± 003 (Figure 2). Assuming a power-law spectral shape, the likelihood analysis of PS 1 resulted in a TS value of 31, spectral index of 2.27 ± 0.10, and an energy flux of (1.77 ± 0.40) × 10⁻¹¹ erg cm⁻² s⁻¹.

Since VERITAS observations have the deepest exposure on MGRO J1908+06 in TeV range and provided most detailed TeV morphology, we included the VERITAS counts map as a template in our GeV morphology analysis. HAWC observations provided the only SED data points on MGRO J1908+06 above 50 TeV. Thus HAWC data are adopted in the multiwavelength SED modeling. We noticed that in ∼1–10 TeV HAWC data points have higher flux than H.E.S.S. (Figure 3). This may due to the fact that Imaging Atmospheric Cherenkov Telescopes (e.g., VERITAS, H.E.S.S.) use blank sky region near gamma-ray sources as background. In case of a large source extension (e.g., MGRO J1908+06), there might be dim emissions in the rim taken as background, leading to a lower source flux level.

XMM-Newton data sets were reduced with the Science Analysis System (SAS, version 16.1.0). Standard pipeline tasks emproc for MOS data were used to process the raw observation data files. XMM-Newton data were also filtered to avoid the periods of hard X-ray background flares.

The Chandra data were reduced using CIAO version 4.7 and CALDB version 4.7.7. We reprocessed the Chandra data to level = 2 and removed periods of high background or flaring appearing in the observations. The XMM-Newton and Chandra spectra of PSR J1907+0602 are shown in Figure 8.