VER J2227+608: A Hadronic PeVatron Pulsar Wind Nebula?

TL;DR

This paper reports the detection of GeV gamma-ray emission from VER J2227+608, suggesting it may be a PeVatron powered by a pulsar wind nebula, with implications for cosmic ray acceleration up to hundreds of TeV.

Contribution

It provides evidence supporting a hadronic origin of gamma rays and identifies VER J2227+608 as a potential PeVatron candidate, advancing understanding of cosmic ray sources.

Findings

Gamma-ray emission is extended and coincident with molecular clouds.

The gamma-ray spectrum fits a single power-law with index 1.90.

The hadronic model suggests proton cutoff energy > 400 TeV.

Abstract

We report the detection of GeV $\gamma$-ray emission from the very-high-energy (VHE) $\gamma$-ray source VER J2227+608 associated with the "tail" region of SNR G106.3+2.7. The GeV $\gamma$-ray emission is extended and spatially coincident with molecular clouds traced by CO emission. The broadband GeV to TeV emission of VER J2227+608 can be well fitted by a single power-law function with an index of 1.90$\pm$0.04, without obvious indication of spectral cutoff toward high energies. The pure leptonic model for the $\gamma$-ray emission can be marginally ruled out by the X-ray and TeV data. In the hadronic model, the low energy content of CRs and the hard $\gamma$-ray spectrum, in combination with the center-bright source structure, suggest that VER J2227+608 may be powered by the PWN instead of shocks of the SNR. And the cutoff energy of the proton distribution needs to be higher than $\sim$ 400 TeV, which makes it an attractive PeVatron candidate. Future observations by the upcoming Large High Altitude Air Shower Observatory (LHAASO) and the Cherenkov Telescope Array in the north (CTA-North) could distinguish these models and constrain the maximum energy of cosmic rays in supernova remnants.

Figures

Full text

gamma rays: ISM - ISM: cosmic rays - individual: (G106.3+2.7 = VER J2227+608) - radiation mechanisms: non-thermal

\AuthorCallLimit

=1 \fullcollaborationNameThe Friends of AASTeX Collaboration

1 Introduction

It is widely believed that supernova remnants (SNRs) are the dominant accelerators of Galactic cosmic rays (CRs). The acceleration of electrons to hundreds of TeV energies has been supported by the detection of non-thermal synchrotron X-ray emission (Koyama et al., 1995). In the $\gamma$-ray band, Fermi-LAT and several ground-based Cherenkov telescopes, like HESS, VERITAS, MAGIC, etc, have detected tens of SNRs. Especially, the $\gamma$-ray emission of SNRs IC 443 and W44 observed by Fermi-LAT provides direct evidence for proton acceleration in SNRs (Ackermann et al., 2013). However, we still lack direct observational evidence for CRs in SNRs up to the spectral knee. Pulsar wind nebulae (PWNe) driven by the central pulsars are also important particle accelerators and a large amount of PWNe have been detected in the $\gamma$-ray band (Abramowski et al., 2012; Xin et al., 2018). High-energy particles can be accelerated at the termination shocks of PWNe. In particular, GeV flares from the Crab nebula show that it can accelerate electrons to PeV energy (Balbo et al., 2011), implying that PWNe could be important contributors to VHE Galactic CRs.

SNR G106.3+2.7 was first discovered by the northern Galactic plane survey of the Dominion Radio Astrophysical Observatory (DRAO) at 408 MHz(Joncas & Higgs, 1990). Pineault & Joncas (2000) later found that it has a comet-shaped radio morphology with the spectral index $\alpha\approx 0.57\pm 0.04$, where the flux density S${}_{\nu}\propto\nu^{-\alpha}$. G106.3+2.7 consists of two distinct regions, a compact “head” region produced by the SNR shock interacting with the ambient dense material and an elongated “tail” region which corresponds to an outbreak to the interior of an neutral hydrogen (HI) bubble created by the stellar wind or supernova explosion (Kothes et al., 2001, 2006b). The “head” region also contains an off-center PWN in the north named “Boomerang” powered by the pulsar PSR J2229+6114, whose characteristic age and luminosity are 10 kyr and 2.2 $\times$ $10^{37}$erg s*-1*, respectively (Halpern et al., 2001b). For the PWN, a spectral break at about 4.3 GHz is detected and has been attributed to the synchrotron cooling (Kothes et al., 2006a). The “tail” region with a lower surface brightness also has a steeper spectrum with $\alpha\approx 0.70\pm 0.07$ than the “head” region ($\alpha\approx 0.49\pm 0.05$; Pineault & Joncas, 2000). Kothes et al. (2001) detected HI and molecular materials associated with SNR G106.3+2.7 and derived a distance of 800 pc, which is much smaller than the value of 3 kpc estimated from X-ray absorption by Halpern et al. (2001b). Moreover, strong linearly polarized radio emission was detected in G106.3+2.7 (Kothes et al., 2001; Gao et al., 2011), which supports a close distance. At a distance of 800 pc, SNR G106.3+2.7 is 14 pc long and 6 pc wide.

In the $\gamma$-ray band, Fermi-LAT has detected GeV emission from PSR J2229+6114 (Abdo et al., 2009a), which is associated with the previously unidentified EGRET source 3EG J2227+6122 (Hartman et al., 1999). Milagro has reported extended VHE $\gamma$-ray emission from G106.3+2.7 (Abdo et al., 2007, 2009b). However, it was not clear whether the emission origins from the SNR or pulsar complex. In Acciari et al. (2009), VERITAS collaboration detected significant TeV $\gamma$-ray emission from the elongated radio extension of SNR G106.3+2.7, $\sim 0.4^{\circ}$ away from the position of PSR J2229+6114 and named it VER J2227+608. The VHE energy spectrum can be well fitted by a power-law of the form $dN/dE=N_{0}(E/3\rm TeV)^{-\Gamma}$, with an index of $\Gamma$ = 2.29 $\pm$ 0.33stat $\pm$ 0.30sys and the Milagro flux is coincide with the high-energy extrapolation of this spectrum (Acciari et al., 2009). Meanwhile, the extended $\gamma$-ray emission is also coincident with molecular clouds traced by 12CO ($J$=1$-$0) emission (Heyer et al., 1998; Kothes et al., 2001), favoring a hadronic origin of the $\gamma$-ray emission.

In this work, we carry out a detailed analysis of GeV emission near VER J2227+608 with ten years of Fermi-LAT Pass 8 data. In Section 2, the data analysis routines and results are presented. The origin of the $\gamma$-ray emission and the expected detection in the TeV band are discussed in Section 3, followed by conclusions in Section 4.

2 Data Analysis

2.1 Data Reduction

In the following analysis, we select the latest Pass 8 version of the Fermi-LAT data with “Source” event class (evclass=128 & evtype=3), collected from August 4, 2008 (Mission Elapsed Time 239557418) to August 4, 2018 (Mission Elapsed Time 555033605). Only events with energies above 3 GeV are selected to avoid a too large point spread function (PSF) in the lower energy band. To suppress the contamination of the Earth Limb, the maximum zenith angle is adopted to be 90*∘. The region of interest (ROI) is a square area of $14^{\circ}\times 14^{\circ}$ centered at the position of VER J2227+608 (R.A., decl.=$337.00^{\circ}$, $60.88^{\circ}$; Acciari et al., 2009). We used the standard LAT analysis software Fermitools, available from the Fermi Science Support Center 111http://fermi.gsfc.nasa.gov/ssc/data/analysis/software/. Meanwhile, the binned likelihood analysis method with gtlike and the instrument response function (IRF) of “P8R3_SOURCE_V2” are adopted to fit the data. The diffuse backgrounds including the Galactic emission and the isotropic component, are modeled according to gll_iem_v07.fits and iso_P8R3_SOURCE_V2_v1.txt222http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html. All sources in the fourth Fermi catalog (4FGL; The Fermi-LAT collaboration, 2019) within a radius of 20∘* from the ROI center and the two diffuse backgrounds are included in the source model. During the fitting procedure, the spectral parameters and the normalizations of sources within distance smaller than $5^{\circ}$ to VER J2227+608 are set to be free, as well as the normalizations of the two diffuse backgrounds.

2.2 Results

In the direction of SNR G106.3+2.7, there is only one 4FGL point source, 4FGL J2229.0+6114, which is the counterpart of PSR J2229+6114 (The Fermi-LAT collaboration, 2019). We first created a $2.0^{\circ}\times 2.0^{\circ}$ Test Statistic (TS) map using gttsmap by subtracting the emission from the diffuse backgrounds and all 4FGL sources in the best-fit model, which is shown in the Figure 1. The TS map shows significant $\gamma$-ray emission toward the “tail” region of SNR G106.3+2.7, which is marked as SrcX afterwards. Then we added SrcX in the model file as a point source with power-law spectrum and redo the likelihood fitting. The best-fit position of SrcX optimized by the gtfindsrc command is R.A., decl.=$336.68^{\circ}$, $60.88^{\circ}$ with 1$\sigma$ error radius of $0.07^{\circ}$. The TS value of SrcX as a point source is 35.5, corresponding to a significance level of 5.1$\sigma$. The spectral index and integral photon flux are 1.60$\pm$0.17 and (9.18$\pm$2.68) $\times\rm 10^{-11}$ photon cm*-2* s*-1*, respectively.

We carried out a spatial extension test for the $\gamma$-ray emission of SrcX using four uniform disks with different radii. These uniform disks are centered at the best-fit position of SrcX and the radii of them are adopted to be $0.15^{\circ}$, $0.20^{\circ}$, $0.25^{\circ}$ and $0.30^{\circ}$. The TS values for the different spatial templates are 46.2, 49.8, 51.5 and 49.7, respectively. The Akaike information criterion (AIC; Akaike, 1974) is adopted to compare the different spatial models, which shows $\Delta$AIC = AICpoint - AICext = 14. Here AICpoint and AICext represent the AIC values for point source and extended uniform disk with radius of $0.25^{\circ}$, respectively. This result shows an obvious improvement for using the extended spatial template than a point-source hypothesis. Meanwhile, the extended $\gamma$-ray emission of SrcX is consistent with the extended TeV $\gamma$-ray emission of VER J2227+608 as shown in Figure 1 (Acciari et al., 2009), which supports SrcX to be the counterpart of VER J2227+608. Hereafter, we adopted the uniform disk with radius of $0.25^{\circ}$ to describe the $\gamma$-ray emission of SrcX in the spectral analysis.

With the extended spatial template for SrcX, the spectral index and integrated photon flux of SrcX in the energy band of 3-500GeV are fitted to be 1.81$\pm$0.16 and (2.04$\pm$0.48) $\times\rm 10^{-10}$ photon cm*-2* s*-1*, respectively. To obtain the spectral energy distribution (SED) of SrcX, the data are further divided into five equal logarithmic energy bins between 3 GeV and 500 GeV. And the likelihood analysis is repeated for each energy bin with only the spectral normalizations of all sources within $5.0^{\circ}$ from SrcX and the two diffuse backgrounds left free. For any bin with TS value of SrcX smaller than 5.0, a 95% confidence level upper limit is calculated. The SED of SrcX is shown in Figure 2, together with the SED of VER J2227+608 (Acciari et al., 2009). Interestingly, the GeV spectrum of SrcX is well consistent with the TeV spectrum of VER J2227+608 and the global spectral index is fitted to be 1.90$\pm$0.04 shown as the blue lines in Figure 2.

3 Discussion

The above data analysis shows that SrcX with the extended $\gamma$-ray emission is positionally coincident with VER J2227+608 and the GeV spectrum of SrcX is also in agreement with the TeV spectrum of VER J2227+608, which all support SrcX as the counterpart of VER J2227+608. In Figure 1, the radio emission of SNR G106.3+2.7 at 1420 MHz by the Canadian Galactic Plane Survey (CGPS; Kothes et al., 2001; Taylor et al., 2003) shows that the $\gamma$-ray emission of VER J2227+608 is spatially consistent with the “tail” region of SNR G106.3+2.7. And in the right panel of Figure 1, the 12CO ($J$=1$-$0) emission integrated between -5.59 and -7.23 km $\rm s^{-1}$ from the Five College Radio Astronomy Observatory (FCRAO; Heyer et al., 1998; Kothes et al., 2001) shows the molecular clouds around VER J2227+608, whose distance of 800 pc is derived from the center velocity of -6.4 km $\rm s^{-1}$ (Kothes et al., 2001). This suggests that the $\gamma$-ray emission of VER J2227+608 may be associated with the molecular clouds.

For the $\gamma$-ray origin of VER J2227+608, several radiation mechanisms are considered here. For the hadronic model, the $\gamma$-ray emission is from the decay of neutral pion mesons produced by inelastic pp collisions. And for the leptonic model, the $\gamma$-ray emission can be produced by inverse Compton scattering (ICS) or a bremsstrahlung process of relativistic electrons. The radio and non-thermal X-ray emission is produced by high-energy electrons via the synchrotron process.

Pineault & Joncas (2000) reported the radio flux densities and spectral indices of the “tail” and “head” regions of G106.3+2.7. And here we treat the “tail” region as the radio counterpart of VER J2227+608 considering the spatial correlation. In the X-ray band, only PSR J2229+6114 and its PWN in the “head” region have been detected and the “tail” region has no significant X-ray emission (Halpern et al., 2001a, b, 2002). Here, we adopted the intrinsic 2-10 keV flux of PSR J2229+6114, 1.3$\times$10*-12* erg cm*-2* s*-1* (Halpern et al., 2001b), as an upper limit for VER J2227+608 in the X-ray band, which can be lower than that obtained via proper data analyses for the expended $\gamma$-ray source region.

Here we quantitatively discuss a pure leptonic model and a hadronic-leptonic hybrid one for the multi-wavelength spectrum of VER J2227+608. And the spectrum of electrons and/or protons is simply assumed to be a power-law spectrum with an exponential cutoff: $dN/dE_{i}\propto E_{i}^{-\alpha_{i}}\exp(-E_{i}/E_{i,\rm cut})$, where $\alpha_{i}$ and $E_{i,\rm cut}$ are the spectral index and the cutoff energy, respectively, for $i=e$ or $p$. The distance of VER J2227+608 is adopted to be 800 pc (Kothes et al., 2001), and the radius of the radiation region is $\sim$ 3.5 pc considering the best-fit uniform disk with radius of $0.25^{\circ}$ for the $\gamma$-ray emission in the data analysis above. Besides the cosmic microwave background (CMB), an infrared (IR) radiation field ($T$ = 30 K, $u$ = 1.0 eV cm*-3*; Porter et al., 2006) is also considered in the ICS process. The gas density is assumed to be $n_{\rm gas}$ = 1.0 cm*-3* for the bremsstrahlung and hadronic processes.

In the pure leptonic model, the spectral index of electrons is about 2.3. And the total energy of electrons above 1 GeV and the magnetic field strength are 4$\times$1046 erg and 8 $\mu$G, respectively. Without considering the TeV data measured by Milagro (Abdo et al., 2007, 2009b), the cutoff energy of electrons need to be less than about 20 TeV to avoid exceeding the upper limit in the X-ray band. With $B$ = 8 $\mu$G and $E_{\rm e,cut}$ = 20 TeV, the synchrotron cooling timescale of electrons is estimated to be about 9400 yr, which is close to the age of G106.3+2.7, $\sim 10^{4}$yr, inferred from the pulsar data (Halpern et al., 2001b). And such a scenario is likely consistent with the expectation about the evolution of high-energy particle distribution in SNRs by Zeng et al. (2017, 2019) and Zhang & Liu (2019a), which suggests that the acceleration of high-energy electrons has stopped at present. However, this leptonic model can not reproduce the flux in the two highest energy bins and it requires that high-energy electrons accelerated by shocks in the early stage of the SNR evolution are still trapped inside the SNR.

For the hadronic-leptonic hybrid model, the spectral indices of electrons and protons are adopted to be 2.0 to reduce the number of parameters. Meanwhile, the ratio of the normalization of electron distribution to that of protons, $K_{\rm ep}$, is fixed to be 0.01, which is in accord with the locally measured CR electron to proton flux ratio (Yuan et al., 2012). With the gas density of $n_{\rm gas}$ = 1.0 cm*-3*, the total energy of protons with the energy above 1 GeV is $\sim 6.0\times 10^{48}(n/1.0\,\mathrm{cm}^{-3})^{-1}\ \mathrm{erg}$. To explain the flat $\gamma$-ray spectra, the cutoff energy of protons should be greater than 400 TeV. The cutoff energy of electrons should be less than $\sim$ 0.1 TeV to avoid excessive the flux in the GeV $\gamma$-ray band. And the magnetic field strength of $\sim 11\mu$G is needed to explain the radio flux by synchrotron process.

Considering the presence of molecular clouds in this region as shown in the right panel of Figure 1, we get an approximate density of about 10.0 cm*-3* based on the values of the column density, the CO brightness temperature, the velocity width and the approximate size of the molecular cloud (Kothes et al., 2001). And we also present the hadronic-leptonic hybrid model with the gas density of $n_{\rm gas}$ = 10.0 cm*-3*. Here neutral gases and dense molecular clumps can serve as targets for hadronic processes. With the cutoff energy of protons $\sim$ 400 TeV, the total energy content of protons above 1 GeV is reduced to be $\sim 6.0\times 10^{47}(n/10.0\,\mathrm{cm}^{-3})^{-1}\ \mathrm{erg}$. Here the cutoff energy of electrons $E_{\rm e,cut}$ is poorly constrained and we fix $E_{\rm e,cut}$ by making the synchrotron radiation loss timescale equal to the age of SNR ($\sim$ 104 yr). And the radio flux gives a magnetic field strength of $\sim 50\mu$G and the corresponding $E_{\rm e,cut}$ is about 0.5 TeV. Based on the observations of the HI and CO emission, Kothes et al. (2001) suggested that the “tail” of G106.3+2.7 was created by a breakout of the SNR shock wave into the HI bubble with a lower density. However, the relatively lower energy content in CRs, the hard $\gamma$-ray spectrum, and its center bright morphology are distinct from typical SNRs (Zeng et al., 2019), implying that this source may be driven by the PWN instead of shocks of the SNR. Nevertheless, if high-energy ions accelerated in the early stage of SNR evolution (Zhang et al., 2017) can be trapped within the SNR, the $\gamma$-ray emission may be attributed to the SNR.

If the tail is due to a breakout of the PWN into the HI bubble, we can get a value of the total ejected mass Mej of the supernova using equation (2) in Chevalier (2004), which describes the dynamical evolution of a PWN. And for PSR J2229+6114 with $\dot{E}$ = 2.2 $\times$ $10^{37}$erg s*-1* and the age of 10 kyr (Halpern et al., 2001b), the value of Mej is about 2.5 M⊙, assuming a typical explosion energy of 1051 erg and the rough radius of 14 pc (the long axis of G106.3+2.7). And this value of Mej seems to be reasonable for typical SNRs Chevalier (2004). The $\gamma$-ray emission from PWN is usually attributed to the ICS process with electrons accelerated by its termination shock. However, in dense environments, relativistic protons can also be produced and contribute to the $\gamma$-ray emission of PWN with the hadronic process (Cheng et al., 1990; Bartko & Bednarek, 2008). Although evidence for the presence of relativistic hadrons in PWNe has been elusive so far, the broadband spectra of several sources have been explained by the hadronic model, e.g. Vela X (Zhang & Yang, 2009), PWN G54.1+0.3 (Li et al., 2010) and DA 495(Coerver et al., 2019). Moreover, Ohira et al. (2018) suggested that PWNe inside SNRs could further re-accelerate the relativistic protons up to the energy of 1 PeV and such PWNe may therefore be important PeVatrons. The hard spectrum of protons and high proton cutoff energy of $>$400 TeV in the hybird model make VER J2227+608 a potential PeVatron. Other TeV sources with very hard spectra may have a similar origin (Guo et al., 2017).

Upcoming $\gamma$-ray experiments LHAASO (Bai et al., 2019; Liu et al., 2016) and CTA-North (CTA Consortium, 2019) can detect $\gamma$-ray emission from VER J2227+608. As shown in Figure 2, the global spectral index of the $\gamma$-ray emission from VER J2227+608 is hard with spectral index of 1.90$\pm$0.04. And the differential sensitivities of LHAASO (1 year) with different sizes of photomultiplier tube (PMT; Bai et al., 2019) and CTA-North (50 hrs; CTA Consortium, 2019) are also overplotted. Especially in the higher energy band (e.g. above 20 TeV), LHAASO has a great advantage than CTA-North. The different radiation models discussed above predict different spectra of VER J2227+608 for the energy above several tens of TeV, and the future accurate measurements of the TeV-PeV spectrum can distinguish these different radiation models and clarify whether it is a PeVatron.

4 Conclusion

In this work, we analyze the $\gamma$-ray emission toward VER J2227+608, which is associated with the “tail” region of SNR G106.3+2.7, using ten years Pass 8 data of Fermi-LAT. The GeV $\gamma$-ray emission of VER J2227+608 is extended and the global GeV to TeV spectra can be described by a hard power-law form with an index of 1.90$\pm$0.04. Both the pure leptonic and the hadronic-leptonic hybrid models are considered for the multiwavelength data of VER J2227+608. The pure leptonic model shows that the synchrotron cooling timescale of electrons is approximately equal to the age of SNR G106.3+2.7, implying quenching of VHE electron acceleration. However, the cutoff energy of electrons derived from the upper limit in the X-ray band seems to contradict with the TeV flux measured by Milagro. The hybrid model seems to be favored considering the spatial coincident between the $\gamma$-ray emission and the likely dense molecular clouds traced by CO emission. The cutoff energy of protons is constrained to be higher than $\sim$ 400 TeV. Such high-energy protons makes VER J2227+608 a promising PeVatron candidate. The low energy content of CRs in this source, its hard $\gamma$-ray spectrum and center bright morphology suggest that it may be powered by the PWN instead of shocks of the SNR. Future observations by LHAASO and CTA-North in the TeV band, especially accurate measurements of the TeV-PeV spectrum by LHAASO, would distinguish the different $\gamma$-ray origins of VER J2227+608 and test its PeVatron nature.

We would like to thank the anonymous referee for very helpful comments, which help to improve the paper significantly. This work is partially supported by the National Key R&D Program of China 2018YFA0404203 and 2018YFA0404204, NSFC grants U1738122, U1931204 and 11761131007, the Natural Science Foundation for Young Scholars of Jiangsu Province, China (No. BK20191109), and by the International Partnership Program of Chinese Academy of Sciences, grant No.114332KYSB20170008.