HAWC J2227+610 and Its Association with G106.3+2.7, a New Potential Galactic PeVatron

A. Albert; R. Alfaro; C. Alvarez; J. R. Angeles Camacho; J. C. Arteaga-Velázquez; K. P. Arunbabu; D. Avila Rojas; H. A. Ayala Solares; V. Baghmanyan; E. Belmont-Moreno; S. Y. BenZvi; C. Brisbois; K. S. Caballero-Mora; T. Capistrán; A. Carramiñana; S. Casanova; U. Cotti; J. Cotzomi; S. Coutiño de León; E. De la Fuente; L. Diaz-Cruz; B. L. Dingus; M. A. DuVernois; J. C. Díaz-Vélez; R. W. Ellsworth; K. Engel; C. Espinoza; K. L. Fan; K. Fang; M. Fernández Alonso; H. Fleischhack; N. Fraija; A. Galván-Gámez; D. Garcia; J. A. García-González; F. Garfias; G. Giacinti; M. M. González; J. A. Goodman; J. P. Harding; S. Hernandez; J. Hinton; B. Hona; D. Huang; F. Hueyotl-Zahuantitla; P. Hüntemeyer; A. Iriarte; A. Jardin-Blicq; V. Joshi; W. H. Lee; H. León Vargas; J. T. Linnemann; A. L. Longinotti; G. Luis-Raya; J. Lundeen; K. Malone; S. S. Marinelli; O. Martinez; I. Martinez-Castellanos; J. Martínez-Castro; J. A. Matthews; P. Miranda-Romagnoli; J. A. Morales-Soto; E. Moreno; M. Mostafá; A. Nayerhoda; L. Nellen; M. Newbold; M. U. Nisa; R. Noriega-Papaqui; N. Omodei; A. Peisker; Y. Pérez Araujo; E. G. Pérez-Pérez; C. D. Rho; D. Rosa-González; E. Ruiz-Velasco; H. Salazar; F. Salesa Greus; A. Sandoval; M. Schneider; H. Schoorlemmer; J. Serna Franco; G. Sinnis; A. J. Smith; R. W. Springer; P. Surajbali; E. Tabachnick; M. Tanner; O. Tibolla; K. Tollefson; I. Torres; R. Torres-Escobedo; F. Ureña-Mena; L. Villaseñor; T. Weisgarber; A. Zepeda; H. Zhou; C. de León; J. D. Álvarez; (HAWC Collaboration)

doi:10.3847/2041-8213/ab96cc

1. Introduction

Although their existence has been known for more than 100 years (Hess 1912), the origins of cosmic rays are still not fully understood. These charged particles, mostly protons and fully ionized nuclei, of extraterrestrial origin have been detected over many orders of magnitude in energy (Tanabashi et al. 2018). According to our current understanding, cosmic rays with energies up to a few PeV (the "knee," a softening in the measured energy spectrum) are thought to originate from sources (accelerators) within our own Galaxy. But what are those sources? Supernova remnants (SNRs) have been proposed as potential sources of Galactic cosmic rays, mainly for two reasons (see e.g., Bell 2013 and references therein). First, the diffusive shock acceleration mechanism provides an efficient way to transfer energy from the exploding supernova shell into cosmic rays and accelerate them to relativistic energies. Second, with a rate of a few supernovae per century, SNRs have a sufficient energy budget to accelerate the bulk of the Galaxy's cosmic rays.

As cosmic rays are deflected by Galactic magnetic fields, their arrival directions at Earth generally do not point back to their sources. However, relativistic protons interacting with gas and dust near their sources can produce GeV–TeV gamma-ray emission. And indeed, at least two SNRs have gamma-ray spectra with a characteristic "pion bump" feature, indicating that they accelerate protons to relativistic energies (Ackermann et al. 2013). But so far, no SNR has been shown to emit gamma rays to hundreds of TeV as would be indicative of a PeVatron (a source capable of accelerating protons to at least PeV energies or higher).

Nonetheless, the search for PeVatrons has not been fruitless. H.E.S.S. has detected evidence for the existence of a PeVatron near the Galactic center (Abramowski et al. 2016), although the source of the high-energy protons there is now yet known. HAWC recently performed a blind search for gamma-ray emission above 100 TeV (Abeysekara et al. 2020). All three >100 TeV sources identified in that study show the characteristic spectral curvature indicative of gamma-ray emission from relativistic electrons rather than protons, and they all have energetic pulsars nearby, which are potential sources of relativistic electrons. Still, more in-depth studies regarding the nature of these sources are in progress.

The SNR G106.3+2.7 is a comet-shaped radio source, with a brighter "head" and an extended "tail" region (Joncas & Higgs 1990; Pineault & Joncas 2000). The pulsar PSR J2229+6114, seen in radio, X-rays, and gamma-rays (Hartman et al. 1999; Halpern et al. 2001; Abdo et al. 2009a), and its pulsar wind nebula (PWN) G106.65+2.96, the Boomerang Nebula (Halpern et al. 2001; Kothes et al. 2006), appear to be contained inside the remnant.

Measurements of the gas structure in the region put PSR J2229+6114 and G106.3+2.7 at a distance of about 800 pc from Earth (Kothes et al. 2001), closer than originally thought. These measurements suggest that PSR J2229+6114 and G106.3+2.7 were created from the same progenitor event, with the asymmetric structure resulting from inhomogeneities in the medium surrounding the progenitor star.

Very-high-energy (VHE; E > 100 GeV) emission from this region has been reported by the Milagro collaboration at 20 TeV (Abdo et al. 2007) and 35 TeV (Abdo et al. 2009b, 2009c; Goodman & Sinnis 2009), and by the VERITAS collaboration in the energy range from 900 GeV to 16 TeV (Acciari et al. 2009). The Milagro source (MGRO J2228+61) is spatially extended and consistent with both the pulsar position and the radio "tail" within uncertainties. The VERITAS source is spatially asymmetric. The centroid of the emission is offset from the pulsar, consistent with a molecular cloud in the radio "tail," and consistent with the Milagro source within uncertainties. The emission detected by VERITAS follows a hard power-law energy spectrum (∝E^−2.3). Milagro reported measurements of the differential flux at 20 and 35 TeV, which are consistent with the extrapolation of the VERITAS spectrum to higher energies within uncertainties.

While PSR J2229+6114 is a well-established GeV gamma-ray source, there is no clear GeV counterpart of G106.3+2.7 (Acero et al. 2016; Ackermann et al. 2017; Abdollahi et al. 2020). Recently, an independent group of authors analyzed 10 yr of Fermi-LAT from 3 to 500 GeV and reported a weak, but significantly detected, GeV gamma-ray source with a hard energy spectrum. This GeV source is best described by a disk morphology with a radius of 0 fdg 25, overlapping with the VHE emission region (Xin et al. 2019).

Xin et al. (2019) also show models of the radio and (V)HE gamma-ray emission from G106.3+2.7. Their study prefers a lepto-hadronic model—with the gamma-ray emission dominated by hadronic processes and the proton energy spectrum extending to at least 400 TeV—over a purely leptonic model.

The previous detections of gamma-ray emission up to tens of TeV make G106.3+2.7 a potential Galactic PeVatron: a source that would be able to accelerate cosmic rays up to PeV energies. It would only be the second such source in our Galaxy.

In this Letter, we report on the HAWC detection of multi-TeV emission coincident with G106.3+2.7. Section 2 describes the main results, including constraints derived on the underlying particle population, for two models, assuming either a hadronic or a leptonic origin of the VHE emission. The potential sources of the particles producing the gamma-ray emission are discussed in Section 3. Future prospects are summarized in Section 4. Details on the data analysis can be found in Appendix A.

2. Results

2.1. Source Search

The blind point-source search described in Appendix A found a local maximum at R.A. = 336 fdg 96, decl. = 61 fdg 05, with a statistical uncertainty of 0 fdg 16 and a systematic uncertainty of 0 fdg 1. The significance map for the search is shown in Figure 1. The excess is well isolated and is inconsistent with background fluctuations at the 6.2σ level (pre-trials), or about 4.3σ post-trials considering HAWC's entire field of view. It is designated as a new source, HAWC J2227+610.

2.2. Morphology

The best-fit position of HAWC J2227+610 is consistent with the VHE detections by VERITAS and Milagro within uncertainties (see Figure 1). It is also consistent with the pulsar position.

In addition to the point-source hypothesis, HAWC J2227+610 is fit with a symmetric Gaussian morphology, with the width left free. The extended model is not preferred over the point-source hypothesis. (The faint "tail" toward the southeast in Figure 1 is not statistically significant.) At 90% confidence level (CL), the Gaussian extent is constrained to be less than $\left({0.232}_{-0.004}^{+0.024}\left(\mathrm{syst}.\right)\right)$ °. The HAWC morphology is consistent with the asymmetric Gaussian morphology measured by VERITAS.

As the morphology of the VERITAS source and the HAWC source are consistent with each other, we assume that the VHE emission seen by the two instruments is due to the same particle population/underlying emission mechanism.

2.3. Gamma-Ray Energy Spectrum

Considering only HAWC data, the energy spectrum of HAWC J2227+610 is well fit by a power law. Including spectral curvature or a cutoff does not significantly improve the likelihood. The fit parameters and lower limits on the cutoff energy are given in Table 1. The best-fit spectral index agrees well with the VERITAS measurement. After scaling the VERITAS data points to account for the emission outside the VERITAS integration region (see Appendix A.2), the energy spectrum measured by HAWC lines up very well with the extrapolation of the VERITAS measurement to higher energies, as well as with the Milagro flux points (see Figure 2), without any indication of a break or a cutoff.

**Figure 2.** VHE gamma-ray energy spectrum of HAWC J2227+610, statistical uncertainties only. VERITAS data points and fit results from Acciari et al. (2009) have been scaled to account for the difference between the size of the emission region and the region over which the spectrum was extracted. Two joint fits to HAWC and VERITAS data are shown: one assuming a pion decay spectrum from a proton population following a power-law energy spectrum (hadronic model), and one assuming the gamma-ray emission is due to inverse Compton emission from electron (leptonic model). Milagro data points from Abdo et al. (2007, 2009c) are shown for reference only.
Download figure:
Standard image High-resolution image

Table 1. Spectral Fit Parameters for Each Data Set as Well as for the Joint Fit

	VERITAS	HAWC	Joint Fit
E₀ (TeV)	3	80	20
K (TeV⁻¹ cm⁻² s⁻¹)	$\left(1.15\pm 0.27\pm 0.35\right)\times {10}^{-13}$	$\left(1.02\pm {0.17}_{-0.22}^{+0.17}\right)\times {10}^{-16}$	$\left(2.46\pm {0.35}_{-0.47}^{+0.33}\right)\times {10}^{-15}$
γ	−2.29 ± 0.33 ± 0.30	$-2.25\pm {0.23}_{-0.19}^{+0.03}$	−2.29 ± 0.08 ± 0.09
E_min (TeV)	0.9	40	0.9
E_max (TeV)	17	110	180

E_C (90% CL) (TeV)		$\gt {35.7}_{-16.9}^{+0.1}$	$\gt {120}_{-76}^{+81}$

Note. VERITAS results are taken from Acciari et al. (2009) (not adjusted emission outside the VERITAS integration region) and shown here for convenience only. For the power-law flux normalization and spectral index, both statistical and systematic uncertainties are given. For the lower limit on the cutoff energy, only systematic uncertainties are given.

Download table as: ASCII Typeset image

The joint VERITAS–HAWC spectrum is well fit by a power law from 900 to 180 TeV. At 90% CL, a lower limit on the cutoff energy is placed at ${120}_{-60}^{+81}\,\left(\mathrm{syst}.\right)$ TeV (assuming an exponential cutoff).

2.4. Primary Proton Energy Spectrum

The VHE gamma-ray emission from HAWC J2227+610 is modeled as a pion decay spectrum, assuming that the underlying proton spectrum follows a power law. The best-fit spectral index of the relativistic proton population is given by ${\gamma }_{p}=-2.35\pm 0.07{\left(\mathrm{stat}.\right)}_{-0.10}^{+0.09}\,\left(\mathrm{syst}.\right)$ . The best-fit normalization (total proton energy above 2 TeV times gas density) is given by $n\,{W}_{p}=(4.4\pm 0.9\,{\left(\mathrm{stat}.\right)}_{-1.2}^{+1.3}$ $(\mathrm{syst}.))\times {10}^{48}$ erg cm⁻³ · (d/(800 pc))², where d is the distance to the source.

An exponential cutoff in the proton spectrum does not significantly improve the likelihood. The 90% lower limit on the proton cutoff energy is given by ${E}_{C,p}\gt {800}_{-640}^{+990}\left(\mathrm{syst}.\right)\,\mathrm{TeV}$ .

2.5. Primary Electron Energy Spectrum

The model was optimized using the same maximum-likelihood procedure as for the hadronic case. The best-fit values of the parameters are ${W}_{e}=\left(4.2\pm 1.2{\left(\mathrm{stat}.\right)}_{-1.5}^{+1.6}\left(\mathrm{syst}.\right)\right)\times {10}^{45}\,\mathrm{erg}\cdot {\left(d/\left(800\mathrm{pc}\right)\right)}^{2}$ and ${\gamma }_{e}=-2.75\pm 0.11{\left(\mathrm{stat}.\right)}_{-0.14}^{+0.15}\left(\mathrm{syst}.\right)$ , where d is the distance to the source. The resulting predicted gamma-ray energy spectrum can be seen in Figure 2. No evidence for curvature or a cutoff in the electron spectrum was found. The 90% lower limit on the electron cutoff energy is given by ${E}_{C,e}\gt {270}_{-170}^{+140}\left(\mathrm{syst}.\right)$ TeV.

The GHz radio emission seen from G106.3+2.7 (see, e.g., Pineault & Joncas 2000) could be explained by synchrotron emission from electrons with energies in the GeV range. Without precise knowledge of the B-field in the region, the radio and TeV data alone are not sufficient to constrain the electron spectrum in the GeV–TeV range. Better constraints on the synchrotron peak, e.g., from X-ray measurements, would be needed. Due to the relatively large size of the emission region, multiwavelength coverage is rather sparse. For this study, no attempt was made to model the emission at radio wavelengths or any other energy range other than TeV gamma-rays.

3. Origin of the VHE Gamma-Ray Emission

Using VHE gamma-ray data alone, we cannot distinguish whether the gamma-ray emission is due to leptonic or hadronic processes. Both the best-fit leptonic and the best-fit hadronic models for G106.3+2.7 have similar values for the optimized likelihood. We use the Bayesian information criterion (BIC; Schwarz 1978) to compare the two models. We find ${\rm{\Delta }}\mathrm{BIC}={\mathrm{BIC}}_{\mathrm{leptonic}}-{\mathrm{BIC}}_{\mathrm{hadronic}}=0.34$ (in favor of the hadronic model), which does not indicate that one of the models is significantly preferred over the other.

Observations at other wavelengths will be needed to ascertain the true origin of the gamma-ray emission. If the observed radiation is of hadronic origin, one might be able to observe a pion-bump signature at hundreds of MeV (Ackermann et al. 2013). On the other hand, if the emission is of leptonic origin, the same electrons are expected to also emit synchrotron emission up to X-ray energies (depending on the ambient magnetic field).

3.1. Hadronic Origin

Assuming a hadronic origin of the observed gamma-ray emission, the underlying population of relativistic protons is described by a hard power-law energy spectrum, with spectral index −2.35 above 2 TeV and no cutoff below at least 800 TeV. What could be the source of these protons?

Assuming a gas density of n ≈ 50 cm⁻³ (see Appendix B), the best-fit normalization of the proton spectrum corresponds to a total energy in protons above 2 TeV of ${W}_{p}=\left(9\pm 2\left(\mathrm{stat}.\right)\pm 3\left(\mathrm{syst}.\right)\right)\times {10}^{46}\,\mathrm{erg}$ . Assuming that the spectrum extends down to 1 GeV without a break, this corresponds to a total energy in protons above 1 GeV of ${W}_{p,\gt 1\mathrm{GeV}}\approx 1\times {10}^{48}\,\mathrm{erg}$ , within a factor of 2 (large uncertainties as no GeV data were considered for this measurement).

Xin et al. (2019) found a total energy in protons above 1 GeV of 6 × 10⁴⁸ erg for n = 1 cm⁻³, which would correspond to ${W}_{p,\gt 1\mathrm{GeV}}\approx 1.2\times {10}^{47}\,\mathrm{erg}$ assuming n ≈ 50 cm⁻³; almost an order of magnitude lower than the HAWC–VERITAS result. This discrepancy can be explained by the fact that they included a new GeV gamma-ray source in the region that they ascribe to the same proton population. Including the GeV emission, their model prefers a harder index of −2 for the proton spectrum and thus contains fewer protons at lower energies.

The most likely source of these protons is diffusive shock acceleration in the SNR G106.3+2.7, which has a sufficient energy budget assuming even a relatively small acceleration efficiency of about 1% (Kothes et al. 2001 derived a lower limit to the total kinetic energy of the SNR shock of about 7 × 10⁴⁹ erg).

However, SNRs are generally only expected to act as PeVatrons for the first few hundred years of their evolution (Aharonian 2013), while G106.3+2.7 could be as old as 10,000 years. We offer two possible scenarios that could explain the observed VHE emission.

3.1.1. Old SNR Scenario

Assuming that the characteristic age of PSR J2229+6114 (10.4 kyr) is close to the real age of the system, G106.3+2.7 is not expected to be an active proton accelerator up to hundreds of TeV at this stage of it evolution. The observed VHE emission would be due to protons previously accelerated to hundreds of TeV by the SNR, which are now diffusing freely in the region and interacting with a molecular cloud outside the accelerating region.

The energy dependence of the diffusion coefficient D(E_p) of protons in the interstellar medium can be approximated by a power law: $D({E}_{p})={D}_{0}{\left({E}_{p}/\mathrm{GeV}\right)}^{\delta }$ , with ${D}_{0}\approx 3\times {10}^{24}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$ and δ ≈ 1 for Kolmogorov diffusion (Aharonian et al. 2012).

The characteristic length scale that a diffusing particle of energy E travels in time t is given by $l=2\,\sqrt{t\,D(E)}$ (Atoyan et al. 1995). Assuming that the source was mainly active for the first few hundred years of its evolution, then t = 10 kyr, which yields l ≈ 600 pc for an 800 TeV proton (and 200 pc for a 2 TeV proton). This distance is much larger than the few parsec size of the emission region, indicating that most of the protons accelerated by the SNR should already have diffused much further out into the Galaxy. However, it is expected that diffusion is suppressed by a factor of 100 or more near SNRs (Fujita et al. 2009), corresponding to a diffusion length of 60 pc or less for an 800 TeV proton, and enough protons could remain to produce the observed emission.

In the case of an "old" SNR, the total energy in relativistic particles would have to be higher than derived above, as only a fraction of the accelerated protons would have diffused into the emission region. The population of freshly accelerated protons must have had a harder spectral index than the one measured here, since higher-energy protons would have had a larger chance of diffusing away from the emission zone.

3.1.2. Young SNR Scenario

The characteristic age of a pulsar is only a valid approximation for its true age if its current spin period is much slower than its initial period, and its braking index is close to 3. This is not necessarily the case for PSR J2229+6114 with its current period of about 51.6 ms (Halpern et al. 2001), as initial periods of even more than 50 ms are not uncommon (Noutsos et al. 2013). Assuming its initial period was close to the current value, PSR J2229+6114 could be much younger than its characteristic age. On the other hand, PSR J2229+6114 is also thought to be the source of the Boomerang PWN G106.6+2.9, which has been estimated to be at least 3900 yr old (Kothes et al. 2006).

There is, however, still a possibility that the SNR G106.3+2.7 is not connected to PSR J2229+6114 and merely happens to lie on the same line of sight as the Boomerang PWN. In that case, the age of G106.3+2.7 is not well constrained, although is is likely to be at least a few hundred years old as there are no historic records of the supernova explosion.

If the SNR is only hundreds of years old, it could still be an active particle accelerator, and the gamma-ray emission would be due to freshly accelerated protons interacting with molecular gas inside or very nearby the acceleration region.

3.1.3. HAWC J2227+610 as a Potential Neutrino Emitter

Assuming that the VHE gamma-ray emission from HAWC J2227+610 is indeed dominated by the decays of neutral pions from pp collisions, it should be a source of VHE neutrinos from the decay of charged pions, which are also produced in pp interactions. These neutrinos could then be detected by suitable detectors such as IceCube (IceCube Collaboration 2013) or ANTARES (Ageron et al. 2011). However, neither experiment has detected significant steady neutrino point sources yet (Albert et al. 2017; Aartsen et al. 2019; Illuminati 2019).

According to Ahlers & Murase (2014), the flux of muon neutrinos (including antineutrinos) should be proportional to the gamma-ray flux and is approximated by

$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{{\nu }_{\mu }}}{{{dE}}_{{\nu }_{\mu }}}=\displaystyle \frac{{E}_{\gamma }}{{E}_{{\nu }_{\mu }}}\,\displaystyle \frac{{{dN}}_{\gamma }}{{{dE}}_{\gamma }},\end{eqnarray} \tag{ 1 }$

where the neutrino energy E_ν and the gamma-ray energy E_γ are related as E_γ ≈ 2 E_ν. Here, we neglected any gamma-ray absorption effects, and assumed equal flavor ratios at Earth due to flavor mixing. Only muon (anti)neutrinos are considered here as IceCube has published effective area files for track-like, muon-induced events.

Using the best-fit VHE energy spectrum from Table 1 and the IceCube IC-86 effective areas for track-like events for the year 2012 from IceCube Collaboration (2018), the expected detection rate for muon neutrinos from HAWC J2227+610 in the IceCube detector is $0.58\,{\mathrm{yr}}^{-1}$ , or about six per decade, with a median neutrino energy of about 6 TeV. While this is likely not a strong enough signal to be picked up by IceCube, future upgrades to the detector (Aartsen et al. 2014) or next-generation neutrino observatories (e.g., Aiello et al. 2019) might be able to detect neutrino emission from this nearby PeVatron.

3.2. Leptonic Origin

If, on the other hand, the observed VHE gamma-ray emission is of leptonic origin, there must be a source nearby that is capable of accelerating electrons to hundreds of TeV. The observed gamma-ray flux levels require at least 4 × 10⁴⁵ erg in electrons above 2 TeV, which the supernova shell should be able to provide.

Another possible source of electrons could be the pulsar PSR J2229+6114 and/or its PWN G106.65+2.96. PSR J2229+6114 has a spindown luminosity of $2.2\times {10}^{37}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ (Halpern et al. 2001). As discussed earlier, it is not clear that its characteristic age of 10.4 kyr is a good estimate for the actual age of the system. Additionally, the pulsar has shown several timing glitches and its braking index is not well measured (Espinoza et al. 2016). Still, assuming that the pulsar is at least 1 kyr old and its energy output has not been less than the currently measured value, the pulsar would have released at least 7 × 10⁴⁷ erg of kinetic energy over its lifetime; a sufficient energy budget to produce the necessary population of electrons. However, if the pulsar is the source of the relativistic electrons producing the gamma-ray emission, one would expect this emission to be centered around the pulsar.

4. Conclusions and Outlook

The HAWC detection of hard-spectrum gamma-ray emission from the G106.3+2.7 region up to more than 100 TeV is interpreted in the context of both a hadronic and a leptonic emission model. Assuming a hadronic origin of the VHE gamma-ray emission observed by HAWC and VERITAS, the cutoff energy in the underlying proton spectrum is constrained to be above 800 TeV. This would make the source a Galactic PeVatron. The supernova shockwave could have released sufficient energy to account for the observed VHE gamma-ray emission.

So far, no hint for a cutoff has been detected in the VHE gamma-ray energy spectrum of HAWC J2227+610. Due to the finite age and size of the SNR, it must cut off somewhere between 120 TeV (the lower limit to the gamma-ray cutoff energy derived from the joint fit) and the PeV range, depending on the magnetic field and the size of the acceleration region. Detecting this cutoff would be important to improve our understanding of this source and how much it could contribute to the Galactic cosmic-ray population at the "knee" of the cosmic-ray spectrum.

Due to its large celestial latitude, HAWC J2227+610 can only be observed at relatively large zenith angles by HAWC and other observatories at similar latitudes. Future upgrades to HAWC's reconstruction algorithms, such as optimizing the energy estimation at large zenith angles or including the additional "outrigger" tanks installed around the main detector array (Marandon et al. 2019), will improve HAWC's sensitivity—especially at high energies. Future efforts could also benefit from cross-calibration of HAWC and VERITAS, which would improve the uncertainty on the spectral index and hence allow better constraints on the cutoff energy.

Several next-generation gamma-ray observatories are currently under development or construction. HAWC J2227+610 will be outside the field of view of both the southern Cherenkov Telescope Array (CTA) observatory (CTA Consortium et al. 2018) and the Southern Wide-Field Gamma-Ray Observatory (SWGO; Abreu et al. 2019), both optimized for high energies and the study of Galactic cosmic-ray sources. CTA North, located at 28.7622°N, will be able to observe HAWC J2227+610. However, the northern CTA observatory will be optimized for the detection of extragalactic sources and will have worse sensitivity³² than HAWC above 20 TeV. CTA North could improve its sensitivity (especially at the highest energies) either by spending more time on this source, or by employing dedicated observing strategies such as observations of the rising or setting source at large zenith angles (see Peresano et al. 2019).

The Large High-Altitude Air Shower Observatory (LHAASO; Di Sciascio 2016) is located at a latitude of 29 fdg 36 and uses the water Cerenkov technique to detect air showers, similar to HAWC. Once complete, LHAASO is expected to have more than an order of magnitude better sensitivity compared to HAWC for a similar run time. LHAASO will be able to make an important contribution to the understanding of the emission spectrum of HAWC J2227+610.

We acknowledge the support from the US National Science Foundation (NSF); the US Department of Energy Office of High-Energy Physics; the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory; Consejo Nacional de Ciencia y Tecnología (CONACyT), México, grants 271051, 232656, 260378, 179588, 254964, 258865, 243290, 132197, A1-S-46288, A1-S-22784, cátedras 873, 1563, 341, 323, Red HAWC, México; DGAPA-UNAM grants IG101320, IN111315, IN111716-3, IN111419, IA102019, IN112218; VIEP-BUAP; PIFI 2012, 2013, PROFOCIE 2014, 2015; the University of Wisconsin Alumni Research Foundation; the Institute of Geophysics, Planetary Physics, and Signatures at Los Alamos National Laboratory; Polish Science Centre grant DEC-2017/27/B/ST9/02272; Coordinación de la Investigación Científica de la Universidad Michoacana; Royal Society—Newton Advanced Fellowship 180385; Generalitat Valenciana, Spain, grant CIDEGENT/2018/034. Thanks to Scott Delay, Luciano Díaz and Eduardo Murrieta for technical support.

Facilities: HAWC - , VERITAS. -

Software: astropy (Astropy Collaboration et al. 2013), 3ML (Vianello et al. 2015), Naima (Zabalza 2015).

Appendix A: Instrument and Data Analysis

A.1. The HAWC Detector

HAWC is a large-field-of-view (FoV), ground-based air shower detector array located at 18°59'42'' N, 97°18'27'' W, in Mexico. It is optimized for gamma-ray astronomy in the TeV regime. The detector, event reconstruction, binning, and background estimation are described in Abeysekara et al. (2017a, 2017b). HAWC's energy threshold depends on the source's decl. and its energy spectrum, and ranges from hundreds of GeV to tens of TeV. HAWC's angular resolution (68% containment radius) improves with the fraction of the array triggered by an air shower and the source elevation, and varies from ∼1° to ∼0 fdg 2 (Abeysekara et al. 2017b). About two-thirds of the sky (8.4 sr, from −26° to 64° decl.) are visible to HAWC, with its >95% duty cycle and an instantaneous FoV covering about 1.8 sr. HAWC data corresponding to 1347 sidereal days of livetime are used for this analysis.

A.2. Data Analysis

The source search follows the method described in Abeysekara et al. (2017a): a putative point source with a power-law energy spectrum (spectral index: −2.5) is moved across the sky. The grid points correspond to the centroids of HEALPix pixels (Gorski et al. 2005) with NSIDE = 1024. For each source position, the flux normalization is fit and a likelihood ratio of the best-fit source+background model ${\hat{{ \mathcal L }}}_{s+b}$ , compared to the background-only model ${{ \mathcal L }}_{b}$ , is calculated. A test statistic (TS) is derived from this likelihood ratio: $\mathrm{TS}=2\,\mathrm{log}\left({\hat{{ \mathcal L }}}_{s+b}/{{ \mathcal L }}_{b}\right)$ . The local maxima of the significance map obtained with this procedure are candidate gamma-ray sources (see Figure 1).

Following the initial source search, the multi-mission maximum likelihood (3ML³³ ) framework (Vianello et al. 2015) with the HAWC accelerated likelihood (HAL³⁴ ) plugin is used to perform likelihood fits to determine the morphology and energy spectrum of the source. The likelihood calculation within the HAL plugin proceeds similarly to previous HAWC publications (Younk et al. 2015; Abeysekara et al. 2017b).

For joint fits with VERITAS data, spectral points from Acciari et al. (2009) are added to the likelihood calculation via a ${\chi }^{2}$ -like likelihood. VERITAS quotes a Gaussian half-width of 0 fdg 27 ± 0 fdg 05 by 0 fdg 18 ± 0 fdg 03 for the spatial extent of the source. Accounting for their angular resolution of 0 fdg 1, about 40% of their observed emission is expected to fall outside of their source region (a circular region with a radius of 0 fdg 32). For the joint spectral fits, the VERITAS differential flux measurements are scaled up by a factor of 1.67 to account for this effect.

In this study, the differential photon fluxes ${dN}/\left({dE}\,{dA}\,{dt}\right)$ are modeled as power laws (PLs):

$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}\,{dA}\,{dt}}=K\,{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{\gamma }.\end{eqnarray} \tag{ A1 }$

Power-law spectra with an exponential cutoff (CPL) are also tested:

$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}\,{dA}\,{dt}}=K\,{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{\gamma }\,\exp \,\left(-\displaystyle \frac{E}{{E}_{C}}\right).\end{eqnarray} \tag{ A2 }$

Here, E is the (true) gamma-ray energy, K is the flux normalization, γ is the spectral index, E₀ is the pivot energy, and E_C is the cutoff energy. The gamma-ray emission is assumed to be isotropic and constant in time over the past two decades.

As described in Abeysekara et al. (2017b), HAWC uses the fraction of photon detectors measuring a signal (f_hit) as a proxy for the gamma-ray energy. f_hit is also correlated with the zenith angle of a given shower and its position with respect to the array, making it nontrivial to recover the energy of a given gamma-ray shower. Per-event energy estimators, introduced in Abeysekara et al. (2019), are not used here due to the high decl. of the source. To calculate the energy range over which the spectral fits to HAWC data are valid, the hard-cutoff method presented in Abeysekara et al. (2017b) is used.

Lower limits on the exponential cutoff energy E_C are determined from a likelihood ratio test. The likelihood of a model with an exponential cutoff in the energy spectrum, ${\hat{{ \mathcal L }}}_{\mathrm{CPL}}$ , is compared to the best-fit model with no cutoff, ${\hat{{ \mathcal L }}}_{\mathrm{PL}}$ . A test statistic is derived from this: $\mathrm{TS}=2\left(\mathrm{log}{\hat{{ \mathcal L }}}_{\mathrm{CPL}}-\mathrm{log}{\hat{{ \mathcal L }}}_{\mathrm{PL}}\right)$ . We then estimate the 90% confidence interval on E_C by scanning over a range of cutoff energies, re-optimizing the spectral parameters of the cutoff model for each point in the scan, and identifying the value E_C where $\mathrm{TS}=-1.64$ .

A.3. Modeling of the Underlying Particle Population

A.3.1. Hadronic Modeling

The observed VHE emission can be interpreted in the context of a hadronic emission model: relativistic protons interact with ambient hydrogen nuclei, producing a cascade of particles including neutral pions, which decay into gamma-rays. The Naima framework (Zabalza 2015) is used to predict the resulting gamma-ray emission spectrum from a given proton population. The parameters of the underlying proton spectrum are fit to the data.

Two spectral models are tested for the proton energy spectrum ${{dN}}_{p}/{{dE}}_{p}$ : a simple power law (Equation (A1)) and a power law with a cutoff (Equation (A2)). Instead of the normalization K at a given pivot energy, the total proton energy W_p is used to describe the normalization of the proton spectrum:

$\begin{eqnarray}&&{W}_{p}=\underset{{E}_{1}}{\overset{{E}_{2}}{\int }}{E}_{p}\,\displaystyle \frac{{{dN}}_{p}}{{{dE}}_{p}}{{dE}}_{p}.\end{eqnarray} \tag{ A3 }$

E₁ and E₂ are the minimum and maximum energies of the proton population used to predict the gamma-ray spectrum.

Care should be taken in selecting the energy range; choosing a too-small energy range, meaning E₁ is too high (or E₂ is too low), will cause Naima to underestimate the gamma-ray flux at low (high) energies. On the other hand, choosing a too-large range, meaning E₁ is too low (or E₂ is too high), can increase the overall uncertainty on W_p and the correlation between W_p and the spectral index, as we might be extrapolating the proton spectrum to energies where it is unconstrained by our measurements. We chose E₁ = 2 TeV as a compromise where the flux over the gamma-ray energy range used here (0.9 TeV to 180 TeV) is underestimated by at most 5% and ${E}_{2}=1\,\mathrm{EeV}$ .

For a given spectral shape of the proton spectrum, the predicted gamma-ray emission is proportional to $n\,{W}_{p}/{d}^{2}$ , where n is the gas density in the emission region and d is the distance between the observer and the source. Accordingly, two of these three parameters have to be fixed during the fit procedure. In this study, d is fixed to the measured value of 800 pc (Kothes et al. 2001), and the best-fit normalization is reported in terms of $n\,{W}_{p}$ .

A.3.2. Leptonic Modeling

Relativistic electrons and positrons (from hereon, electrons) emit electromagnetic radiation via three main processes (see, e.g., Blumenthal & Gould 1970 for more details): synchrotron emission due to deflection by ambient magnetic fields, bremsstrahlung due to scattering with ambient ions/nuclei, and inverse Compton (IC) upscattering of lower-energy photons. In most cases, TeV gamma-ray emission from electrons is dominated by the IC process. Using Naima, we model the observed VHE emission from HAWC J2227+610 as IC emission from three seed photon fields: the cosmic microwave background, a galactic near-infrared photon field, and a galactic far-infrared photon field. For all three fields, the default values set by the Naima package are used. As in the hadronic model, the energy spectrum of the electrons is modeled as a power-law spectrum (Equation (A1)) between 2 TeV and 1 EeV, with the index γ_e and the total electron energy W_e being the only two free parameters of the fit. A power-law spectrum with an exponential cutoff (Equation (A2)) was tested as well.

A.4. Systematic Uncertainties

Two classes of systematic uncertainty are considered here: the modeling of the HAWC instrument response and the uncertainty of the VERITAS measurements. Uncertainties related to the HAWC detector model are investigated as in Abeysekara et al. (2019). To incorporate the systematic uncertainties on the VERITAS data, the joint fit is repeated four more times, with the VERITAS data points adjusted according to the quoted systematic uncertainty: all fluxes scaled up/down by the constant factor $\left(1\pm {\rm{\Delta }}K/K\right)$ , and scaled individually by an energy-dependent factor ${\left(E/{E}_{0}\right)}^{\pm }{\rm{\Delta }}\gamma$ . Here, E₀ designates the pivot energy, ΔK is the systematic uncertainty on the flux normalization K, and Δγ is systematic uncertainty on the spectral index. The resulting shifts in the fit parameters were added in quadrature, as in the "shift method" from Heinrich & Lyons (2007).

Appendix B: Molecular Gas in the Region

Here, we briefly outline how the density of molecular gas in the G106.3+2.7 region is determined.

We used data from the DAME CO survey (Dame et al. 2001), available online.³⁵

The available data set contains the brightness temperature T of the ${}^{12}{{\rm{C}}}^{16}{\rm{O}}\,J=1\to 0$ line, with a velocity resolution of $1.3\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and a grid spacing of 0 fdg 125, potentially interpolated from sparser measurements.

Scanning a square of width 0 fdg 6 around the nominal position of HAWC J2227+610, we find a peak in the brightness temperature at −5.2 km s⁻¹, presumably corresponding to the SNR. For the gas maps of the region, we therefore consider the three velocity bins centered at −6.5, −5.2, and −3.9 km s⁻¹.

The column density of hydrogen molecules can be determined from the measured temperature brightness as

$\begin{eqnarray}&&{N}_{{{\rm{H}}}_{2}}={X}_{\mathrm{CO}}\,{\int }_{{v}_{\min }}^{{v}_{\max }}T(v){dv}.\end{eqnarray} \tag{ B1 }$

The CO-to-H₂ conversion factor is given by ${X}_{\mathrm{CO}}\approx 2\times {10}^{20}\,{\mathrm{cm}}^{-2}\,{{\rm{K}}}^{-1}$ according to Bolatto et al. (2013). Here, the integration is replaced by the summation over the previously indicated velocity bins.

Figure 1 shows the column density of hydrogen molecules determined in this way. There is a cloud of molecular gas centered around R.A. = 336 fdg 8, decl. = 60 fdg 85, consistent with the VHE emission region. Integrating the molecular hydrogen in a circle of radius 0 fdg 3 around this point yields a gas content of roughly 400 M_⊙. Assuming a spherical region of radius $800\,\mathrm{pc}\cdot \tan (0\buildrel{\circ}\over{.} 3)\approx 4\,\mathrm{pc}$ , this corresponds to an average density of hydrogen atoms of about 50 cm⁻³, which we adopt in the calculation of the proton energy needed to produce the observed gamma-ray emission. We assume here that the molecular gas content dominates the emission region over atomic or ionized hydrogen.