Investigating the nature of MGRO J1908+06 with multiwavelength observations

The unidentified TeV source MGRO J1908+06, with emission extending from hundreds of GeV to beyond 100TeV, is one of the most intriguing sources in the Galactic plane. MGRO J1908+06 spatially associates with an IceCube hotspot of neutrino emission. Although the hotspot is not significant yet, this suggests a possible hadronic origin of the observed gamma-ray radiation. Here we describe a multiwavelength analysis on MGRO J1908+06 to determine its nature. We identify, for the first time, an extended GeV source as the counterpart of MGRO J1908+06, discovering possibly associated molecular clouds (MCs). The GeV spectrum shows two well-differentiated components: a soft spectral component below $\sim10$ GeV, and a hard one ($\Gamma\sim1.6$) above these energies. The lower-energy part is likely associated with the dense MCs surrounding the supernova remnant SNR G40.5$-$0.5, whereas the higher-energy component, which connects smoothly with the spectrum observed in TeV range, resembles the inverse Compton emission observed in relic pulsar wind nebulae. This simple scenario seems to describe the data satisfactorily, but raises questions about the interpretation of the emission at hundreds of TeV. In this scenario, no detectable neutrino flux would be expected.


Introduction
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). Besides its possible PeVatron nature, MGRO J1908+06 is spatially associated with an IceCube neutrino emission hotspot (Aartsen et al. 2019;Aartsen et al. 2020), though the post-trial significance is low.

Multiwavelength observations
Multiwavelength observations were investigated in this paper, 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 Appendices.
We show the 12 CO (J= 1-0) and 13 CO (J= 1-0) maps for five consecutive velocity ranges from 46 km s −1 to 66 km s −1 with a coverage of 4 km s −1 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 −3 . In the 46 km s −1 to 66 km s −1 velocity slices, there are apparent filaments of MCs positionally located next to PSR J1907+0602. We have inspected the 12 CO (J=1−0) and 13 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;. However, we do not find any such evidence of asymmetric broad profiles of the 12 CO line (i.e., a wing part deviating from the main Gaussian in a 12 CO profile where the signal-noise-ratio of the 13 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.

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 years 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.
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.=286.61 • ±0.08 • , Decl.= 6.44 • ±0.07 • , and a radius of 0.83 • ±0.05 • , 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 ) −Γ cm −2 s −1 MeV −1 ) and a log-parabola spectral model . 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 −11 erg cm −2 s −1 between 0.1−300 GeV 2 . Adopting the disk morphology, the best-fitted parameters are shown in Table 1 and the gammaray 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.
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.=286.88 • ±0.05 • , Decl.= 6.29 • ±0.05 • with a radius of 0.51 • ±0.02 • , 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-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 −11 erg cm −2 s −1 (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-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  The peak of gamma-ray emission is spatially consistent with the 12 CO (J= 1-0) intensity in the 62 -66 km s −1 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.

Search for MGRO J1908+06 counterparts in X-ray and radio wavelengths
The XMM-Newton X-ray satellite has covered MGRO J1908+06 with 5 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 backgroundsubtracted 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 to 10 keV energy range. The morphology measured by H.E.S.S. (a Gaussian profile with σ of 0.34 • ; 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 22 cm −2 (Abdo et al. 2010a), we calculated a 95% unabsorbed upper limit for MGRO J1908+06 as 1.2 × 10 −10 erg cm −2 s −1 (0.2−10 keV).   Off-peak Analysis from 0.1-300 GeV  Figure 1. The x and y axes are R.A. and decl. (J2000) in degrees. Top right: GeV SED of Fermi J1906+0626. Data points in black and blue are from the PSR J1907+0602 off-peak and phase averaged analyses, respectively (see text for detail). The dashed red and blue curves on two instances indicate the two-component spectral modelling (Table 1)   From our CO observation, the local H i column density to the MGRO J1908+06 region is N H = 2×N H 2 ∼ (2 − 7) × 10 21 cm −2 (Appendix A, Table 2), which is much lower than the total Galactic H i column density in this direction estimated using the HEASARC tool 3 and with Chandra ( We checked for variability with spectral analysis of the eight-month separated XMM-Newton (Obs. ID 0605700201 on 2010-04-26) and Chandra (ObsID 7049 on 2009-08-19) observations. We adopted a simple power-law model plus absorption with H i column density fixed at N H = 1.3×10 22 cm −2 (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.42 −0.41 with unabsorbed total energy flux of 6.29 +1.44 −1.18 ×10 −14 erg cm −2 s −1 , whereas the Chandra ACIS spectrum is well fitted with compatible spectral parameters, a best-fit spectral index of 1.28 +0.45 −0.42 and unabsorbed total energy flux of 5.14 +1.47 −0.94 ×10 −14 erg cm −2 s −1 . No spectral or flux variability can be claimed. The XMM-Newton and Chandra spectra are shown in Appendix C, Figure 8.
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 extract photons from PN data using a radius of 20 arcsec in 0.2−10 keV and 3−10 keV with barycenter correction applied. Using Tempo2 (Hobbs et al. 2006) and the photons plugin 5 and contemporaneous Fermi-LAT gamma-ray ephemeris, we have assigned pulsar rotational phase to each extracted photon. No significant X-ray pulsation is detected.

Discussion
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 Cherenkov instruments (Figure 3 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-modifiled 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 modelling of secondaryto-primary CR ratios (Aguilar et al. 2016;Genolini et al. 2019, Huang et al. 2020. We assume an average hydrogen density of n = 45 cm −3 in the surrounding MC, resulting in a pp collision cooling timescale of t pp 0.5σ pp nc −1 = 7 × 10 5 (n/100cm −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 50 ergs, 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 Table 2), and this change reduces the requirement for the proton energy to 10 49 ergs. for E e ≥ E b , which is usually chosen to describe the SED of PWNe (Tanaka & Takahara 2010;Bucciantini et al. 2011;Martin et al. 2012;Torres et al. 2013;Torres et al. 2014). The time-independent spectra of synchrotron radiation and IC radiation are calculated following Blumenthal (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 modelled in Popescu et al. (2017) as well as the cosmic microwave background as the target photon field for the IC scattering. The Fermi-LAT 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 47 erg in the emitting electrons/positrons 7 . 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 cool 40(E/8.2TeV) −1 (U ph + U B ) −1 /1 eV cm −3 −1 kyr. Given U ph = 1 eV cm −3 and U B = B 2 /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 . 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 understand the transport mechanism of particles within the PWN and the evolution of the PWN (H.E.S.S. Collaboration 2019; Liu & Yan 2020). Deep observations with TeV observatories such H.E.S.S., HAWC or LHAASO will provide crucial input to disentangle the origin of the gammaray emission observed.
We also note that there are actually many relevant physical processes which can influence the modelling, such as the particle injection history, particle spectral evolution and particle transport. We leave such a more realistic modelling 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 which 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 clear for the time being that whether the gamma-ray spectrum above 100 TeV would decline as our model expectation. IceCube Collaboration (Aartsen et al. 2020) found that the best-fit ν µ +ν µ 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 ν = 5.7×10 −13 (E ν /1TeV) −2 TeV −1 cm −2 s −1 . In the pp collision, the pionic gamma-ray flux is about twice that of the ν µ +ν µ flux. Therefore, an upper limit for the accompanying pionic gammaray component can be obtained as shown in the dashed grey 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 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.
the molecular column density of N(H 2 ) ≈ 7 × 10 5 N( 13 CO) (Frerking et al. 1982) has been used. The mean density of each region is calculated by dividing the column density toward the 13 CO emission peaks by the cloud size along the line-of-sight, which is estimated from the FWHM of the 13 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 −3 .

Appendix B: Fermi-LAT data analysis
The analysis shown in this paper uses more than 11 years of Fermi-LAT P8R3 data, from 2008 August 4 (MJD 54682) to 2019 November 09 (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 0.1 • 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 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 ln(L point /L 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  Four regions delineated in green and labelled with roman numerals "1" to "4" are used to estimate the astrophysical parameters for the molecular gases (see Table 2). Lower panels: 12 CO (J= 1-0 black) and 13 CO (J= 1-0, blue) spectra of region "1" to "4". 3.5 56.9 3.9/9.0 Chandra ACIS spectrum. The best-fitted models and post-fit residuals are also shown in each case.