© The Authors 2024

1. Introduction

Cosmic-ray (CR) accelerators have been probed indirectly across a wide range of energies, from the radio to ultra-high-energy gamma rays. Galactic CRs are expected to be accelerated up to the knee of the CR spectrum, meaning they can reach up to PeV energies (Cristofari 2021). It is suggested that these PeV CR accelerators (PeVatrons) are most commonly supernova remnants (SNRs), where charged particles are accelerated in the shock fronts (Cristofari et al. 2018; Reynolds 2010; Bell et al. 2013).

SNR G106.3+2.7 is a PeVatron candidate that has been under study for many years now. It was first detected in 1990 by Joncas & Higgs (1990) at 408 MHz with the Dominion Radio Astrophysical Observatory (DRAO). The head of the SNR is home to the pulsar J2229+6114 and its boomerang-like pulsar wind nebula (PWN) (Kothes et al. 2001), known as the Boomerang region. An extended shell then reaches into the tail region of the SNR (Pineault & Joncas 2000). There have been many proposed distances for this region. Pineault & Joncas (2000) adopted a distance of 12 kpc for the SNR based on a flat rotation curve and Halpern et al. (2001) found the distance to be 3 kpc using X-ray absorption column density. HI observations from Kothes et al. (2001) place the system at 800 PC.

X-ray observations show similar structures coincident with the head and tail regions. Fujita et al. (2021) finds that the X-ray emission from the entire SNR is from nonthermal synchrotron radiation. However, the X-ray brightness in the head region contrasts with the gamma ray observations, which instead see a brightness in the tail region.

Ground-based gamma ray particle detector arrays have so far only been able to resolve one emission spot in the Boomerang region. Milagro, VERITAS, Tibet, and LHAASO report extended emission coincident with the tail of the SNR. Fermi-LAT (Fang et al. 2022) and MAGIC (MAGIC Collaboration 2023) report gamma ray emission from the head and tail regions of the SNR separately, but no other gamma ray observations have been used to disentangle the sources thus far. Fang et al. (2022) separates the pulsar emission from the SNR and performs a detailed analysis of the tail region using Fermi-LAT data. They find that pion decay (PD) is the primary mechanism for producing the high-energy emission seen in the tail region.

MAGIC Collaboration (2023) reports 0.16° extended gamma ray emission from the head and tail regions. They also show that the emission shifted in the SNR depending on the energy threshold. Their highest energy emission (6.0–30 TeV) is consistent with the tail of the SNR. Their multi-wavelength analysis reveals that the head region can have a leptonic or hadronic explanation, while the tail region is most likely of hadronic origin.

The previous HAWC detection analyzes 1350 days of data using older reconstruction algorithms (Albert et al. 2020). Albert et al. (2020) finds point-like emission consistent with the other gamma ray observations, which see extended emission in the tail region. The work presented in this paper uses more data and improved reconstruction algorithms (Albert et al. 2024), which further extends the previous HAWC measurement past 110 TeV with flux points. We also explore a molecular cloud template fit using data from the Planck 353 GHz dust opacity map and estimate the CR flux in the region using distance information from the Dame CO Survey (Dame et al. 2001). Section 2 describes the HAWC observatory and data used for this analysis. Section 3 describes the analysis procedure, as well as the results from the procedure. Gamma ray production mechanisms and their probable CR sources are described in Section 4. Conclusions and further outlooks are presented in Section 5.

2. HAWC Observatory and data

The High-Altitude Water Cherenkov (HAWC) Observatory is located in Puebla, Mexico at an altitude of 4100 m. HAWC utilizes the water Cherenkov detection technique and consists of 300 water Cherenkov detectors (WCDs). Each WCD contains four photo-multiplier tubes (PMTs) and the entire array covers ∼22 000 m² (Abeysekara et al. 2023). For this analysis, we used 2565 days of data, which is about 1000 days more data than the previous publication Albert et al. (2020). A neural net (NN) algorithm was used to reconstruct the energy of the primary gamma ray that initiated the extensive air shower. From here, the data was binned in a 2D binning scheme based on the fraction of PMTs triggered and reconstructed energy (Abeysekara et al. 2019). In addition to having more data, the new data reconstruction process used algorithms that better optimize angular resolution.

3. Analysis and results

3.1. Systematic source search

This first part of this analysis utilized the Multi-Mission Maximum Likelihood (threeML) framework (Vianello et al. 2015), along with the HAWC Accelerated Likelihood (HAL) plugin (Abeysekara et al. 2022) to work through a step-by-step multi-source search method inspired by the Fermi-LAT Extended Source Catalog (Ackermann et al. 2017; Albert et al. 2023). ThreeML uses a likelihood ratio as a test statistic (TS) to reach a parameter set that maximizes the likelihood of the hypothesized model over the background-only model. The value of the TS is defined as follows:

$$\begin{array}{c}\hfill \mathrm{TS}=2ln\frac{L_{\mathrm{alt}}}{L_{\mathrm{null}}},\end{array}$$(1)

where L_alt is the model hypothesis and L_null is the background only hypothesis. In order to arrive at the most statistically preferred model in our multi-source search method, we looked at the difference in TS between each nested model that is tested. There are two main components to this analysis described below: point source search and extension testing.

Point Source Search. In this section, source TS is calculated assuming L_alt to be the current model and L_null to be the previous model.

1. The first step of our multi-source search is to add one point source with a power-law (PL) spectrum to the location of the highest TS value in the significance map, leaving the position and spectral values of the newly added source free while the position of the previous sources remains fixed.

2. After fitting the model in step 1, we check the ΔTS between the previous model and the most recently tested model:

   1. If the ΔTS > 25, we keep the source in the model and return to step 1.

   2. If the ΔTS < 25, we do not add the additional source to the model and stop adding point sources. We then start our extended source testing.

Extended Source Testing. In this section, source TS is calculated assuming L_alt to be the entire model and L_null to be the entire model with the source in question subtracted out.

1. We start by replacing the highest TS point source with a 2D extended Gaussian and leave all other sources as point sources. We keep all source locations fixed.

2. After fitting the model in step 1, we check the ΔTS between the previous model and the most recently tested model:

   1. If the ΔTS > 16, we accept the extended source model and go to step 3.

   2. If the ΔTS < 16, we reject the extended source model and move to the next highest TS value source. We then start back at step 1.

3. Next, we check the TS values of the point sources in the model:

   1. If the TS values of all point sources are > 25, we refit the model and float all parameters. If untested point sources remain, we go back to step 1. We end the study if there are no more point sources remaining for extension testing.

   2. If the TS of any point source is < 25, we remove the source(s) and refit the new model while floating all parameters. If untested point sources remain, we go back to step 1. We end the study if there are no more point sources remaining for extension testing.

Once all point sources in the final point-source model are tested for extension, the source search ends. Figure 1 visualizes the point and extended source systematic steps. After completing the source search, we then test each source for a curvature in its spectrum using a log-parabola (LP), which is defined by:

Fig. 1. Diagram of the multi-source search method described in Section 3. Yellow boxes (first two columns) describe how point sources (PSs) are added, where PL refers to a power-law spectrum. Orange/red boxes (last four columns) describe how extended sources are added.

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}E}=N_0(\frac{E}{30\mathrm{TeV}}{)}^{\alpha -\beta log\left(\frac{E}{30\mathrm{TeV}}\right)},\end{array}$$(2)

where N₀ is the normalization parameter, the pivot energy is 30 TeV, α is the spectral index, and β dictates the curve in the spectrum. Since the simple power-law is nested within the LP, we can use the same TS comparison to decide if a source’s spectrum is curved or not. We start the curvature testing with the brightest source. If the new model (with the curved spectrum) has a ΔTS > 16, we keep it and move on to test the next source. If the new model (with the curved spectrum) has a ΔTS < 16, it is rejected and we test the next source. Once all sources have been tested for curvature, we refit the final model to ensure the best-fit parameters are further optimized.

After completing this procedure in the Boomerang region, the best-fitting model was found to be a single extended source with an LP spectrum. Table 1 lists the best-fit position and spectral parameters with their statistical and systematic uncertainties. The HAWC significance maps using the new reconstruction algorithms and best-fit values are shown in Figure 2. The systematic uncertainties were determined using various detector response files that test possible effects on the data resulting from differences in the modeling of the instrument. These differences are calculated following Abeysekara et al. (2019).

The next step is to determine the energy range in which the source is confidently detected. The energy range was calculated using the method described in Abeysekara et al. (2017). The energy range at a 1σ confidence level, which corresponds to a 68% confidence limit, is 4–145 TeV. The flux points were calculated by first binning the data into energy bins and determining the median energy in a given bin. We then used threeML to fit the normalization parameter N₀ at the median energy (E) in Equation (2) and all other parameters are fixed at their best-fit values from Table 1 (Abeysekara et al. 2019).

The position (Figure 2b) and spectrum (Figure 3) of the new HAWC observation are consistent with the other gamma ray observations in this region. Amenomori et al. (2021) reports a source extension of 0.24${}_{\pm 0.1(\mathrm{sys})}^{\pm 0.1(\mathrm{stat})}$ degrees (shown in Figure 2b) and LHAASO uses an extension of 0.3 deg in their high-energy source catalog (shown in Figure 2b), as well as reports an extension of 0.35 ± 0.01 (WCDA) and 0.25 ± 0.02 (KM2A) degrees in their first catalog (Cao et al. 2021, 2024). All of these measurements are within the systematic uncertainties of the extension found in this work (Table 1). The main differences between the various gamma ray observations are most likely attributed to each observatory detecting different contributions from the head and tail regions across varying energy ranges.

Fig. 2. Significance maps for HAWC’s all-energy emission and its comparison to the SNR in the region. a) The HAWC all-energy significance map of the region with the 4, 6, and 8σ contours overlaid in green. b) The HAWC significance map of the region with labels showing the best-fit positions from this analysis, the previous HAWC publication (Albert et al. 2020), and VERITAS (Acciari et al. 2009). The green dashed, white dotted, cyan dot-dashed, and blue solid contour indicate LHAASO’s analysis region (Cao et al. 2021), Tibet ASγ’s 1σ extension (Amenomori et al. 2021), Fermi-LAT’s 1σ extension (Fang et al. 2022), and MAGIC’s analysis regions (MAGIC Collaboration 2023), respectively. c) The brightness temperature image (Taylor et al. 2003) of the SNR overlaid with the HAWC 4, 6, and 8σ contours in white. The bright emission on the left shows the head of the SNR, while the faint emission extending to the right shows the tail.

Fig. 3. Spectral energy distribution (SED) of the ground-based VHE gamma ray observations. The green SED with circle flux points are the all-energy HAWC results from this work. The magenta dotted line is the HAWC result for the > 56 TeV molecular cloud template analysis. The red line (KM2A) and blue line (WCDA) are from the first LHAASO catalog (Cao et al. 2024). All SEDs show only statistical uncertainties. The open stars, right-facing triangles, left-facing triangles, squares, and x-marks correspond to the VERITAS (Acciari et al. 2009), MAGIC tail (MAGIC Collaboration 2023), MAGIC head (MAGIC Collaboration 2023), Tibet ASγ (Amenomori et al. 2021), and LHAASO (Cao et al. 2021) measurements, respectively.

3.2. HAWC data > 56 TeV

We will now explore the possible production mechanisms that would contribute to the PeVatron’s UHE emission in this region. Most PeVatron theory assumes that these systems are hadronic, which is confirmed by correlation with molecular clouds in the area. However, recent studies suggest that the leptonic scenario is more than capable of producing UHE emission as well (Vannoni et al. 2009; Breuhaus et al. 2021). Both scenarios will be considered for this region.

3.2.1. Molecular cloud template fitting

We will first look at the UHE gamma ray emission in the region using the pion decay production mechanism. The head of the SNR houses the pulsar and its PWN, which would mean that the protons are most likely being accelerated in the tail region where there is supernova (SN) ejecta and the shock front, which is coincident with a nearby molecular cloud. The molecular cloud is an ideal target for the accelerated protons to interact with and undergo PD. This means that UHE gamma ray morphology should take the shape of the molecular clouds that it is being produced in. This analysis explores fitting HAWC data above 56 TeV using a molecular cloud template from the Planck Collaboration (Ade et al. 2011) to explain a possible hadronic mechanism.

We used two different surveys to assess the molecular clouds in the region: the Dame CO survey (Dame et al. 2001) and the 353 GHz Planck dust opacity map (Ade et al. 2011). Figure 4c shows significant CO emission from the Dame CO survey coincident with HAWC’s > 56 TeV 3, 4, and 5σ contours. Ultimately, we used the Planck dust map for the template fitting because it is more recent than the Dame survey (Dame 2001 vs. Planck 2011) and provides information for HI gas that is optically thick in CO. However, the Planck dust map has no distance information, so we use the Dame CO survey to confirm that there is no emission in front of or behind the Boomerang region that is being included in the Planck dust map. Appendix A gives a more in-depth explanation of how the Planck dust map was validated for use.

Fig. 4. Significance maps for HAWC’s > 56 TeV emission and its comparison to the SNR and molecular clouds in the region. a) The HAWC significance map of the region above 56 TeV. b) The brightness temperature image (Taylor et al. 2003) of the SNR overlaid with the HAWC 3, 4, and 5σ contours in white. c) Molecular hydrogen column density integrated over a velocity range of −20 km/s to 0 km/s (Dame et al. 2001). d) The Planck 353 GHz template normalized to 1/sr that is used in the UHE HAWC fit.

Again, we utilized the threeML/HAL framework to do a maximum likelihood fitting of the molecular cloud spatial template. A similar procedure to that of Albert et al. (2021) was used to create the spatial template that is fed to threeML for the template fitting:

1. We start by calculating the column density– this calculation is taken from Albert et al. (2021):

   $$\begin{array}{c}\hfill N_{\mathrm{H}}={\tau }_{\mathrm{D}}/(\frac{{\tau }_{\mathrm{D}}}{N_{\mathrm{H}}}{)}^{\mathrm{ref}},\end{array}$$(3)

   where (τ_D/N_H)^ref = (1.18 ± 0.17)×10⁻²⁶ cm² for the 353 GHz Planck map (Ade et al. 2011) and τ_D is the dust opacity.

2. We then crop the data to the size of the region of interest.

3. We calculate the mass of the region using the following equation:

   $$\begin{array}{c}\hfill M_{\mathrm{dust}}=N_{\mathrm{H}}\mathrm{\Omega }{d}^2{m}_{\mathrm{h}},\end{array}$$(4)

   where Ω is the angular area of the cloud and d = 800 pc as the distance to the cloud (Kothes et al. 2001)

4. Finally, we normalize the data to 1/sr. This is required for a 3ML 2D template.

Since Kothes et al. (2001) used HI observations to determine the distance to the region, we adopted the same distance in this paper for our calculations. Figure 4.d shows the Planck dust map template with HAWC’s > 56 TeV 3, 4, and 5σ contours. Appendix B gives additional information and important caveats for the template itself. This template, along with a power-law spectrum and a fixed index of −3.0, was used to fit HAWC’s UHE emission shown in Figure 4a. The spectral index for this fit was chosen to be slightly softer than the all-energy measured index (−2.76) as a way to capture the curvature that causes a sharp dropoff at higher energies (Figure 3). The best-fit normalization at 50 TeV for the template fit is 2.8${}_{-0.7}^{+0.9}\times {10}^{-16}$ (1/TeV s cm²).

3.2.2. Simple leptonic model fitting

We also explore the possibility of a leptonic mechanism producing the UHE HAWC data. This would mean that electrons are being accelerated by the SN ejecta or the PWN winds. For this scenario, we modeled the region with a point source morphology and a power-law spectrum, which has a fixed index of −3.0. Table 2 lists the best-fit values for the position and flux normalization after fitting the emission shown in Figure 4.a with the simple model.

4. Scenarios for gamma ray production

4.1. Cosmic ray energy density from protons

In Section 3.2.1, we explored a hadronic-related fit to the HAWC UHE data with a molecular cloud template. We estimated the CR energy density using purely hadronic interactions and the measured gamma ray flux. The CR energy density can then be used to roughly determine the amount of energy that would be needed for the CR population in the region. We started by calculating the gamma ray energy flux using the best-fit gamma ray spectrum from our molecular cloud template fit:

$$\begin{array}{c}\hfill J={\displaystyle {\int }_{\mathit{Ei}}^{\mathit{Ef}}}E\frac{\mathrm{d}N}{\mathrm{d}E}\mathrm{d}E,\end{array}$$(5)

where

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}E}=k(\frac{E}{E_{\mathrm{piv}}}{)}^{-\alpha },\end{array}$$(6)

and E_i = 56 TeV, E_f = 316 TeV, and k, E_piv and α are the same best-fit values from the fitting in Section 3.2.1. Next, we used the integrated flux to calculate the luminosity:

$$\begin{array}{c}\hfill L(\ge E_{\gamma })=4\pi d^{2}J,\end{array}$$(7)

where d = 800 pc. Finally, we used the luminosity to calculate the CR proton energy density using the same approach as Abramowski et al. (2016):

$$\begin{array}{c}\hfill {\omega }_{\mathrm{CR}}(\ge 10E_{\gamma })\approx 1.8\times {10}^{-2}(\frac{{\eta }_N}{1.5}{)}^{-1}(\frac{L(\ge E_{\gamma })}{{10}^{34}\mathrm{erg}/\mathrm{s}})(\frac{M}{{10}^6{M}_{\odot }}{)}^{-1},\end{array}$$(8)

where η_N = 1.5 for heavier nuclei, M = 0.23 × 10⁵ M_⊙ and ω_CR is in eV/cm³, which is calculated the same way as in Albert et al. (2021). Here, we assumed that the proton energy scales as 10× that of the gamma ray energy (Cristofari et al. 2018), so we are probing protons that are > 560 TeV. The CR density was found to be 10.3 × 10⁻³ eV/cm³ for protons > 560 TeV.

The SNR has a length of 14 pc and a width of 6 pc (Kothes et al. 2001). Since the tail is more elongated than the head, we assumed that the tail is two-thirds of the length of the SNR. We also assumed that the SNR is also capable of expanding 14 pc in all other directions as well, so we use a cube of 6 pc × 9.33 pc × 28 pc as the volume of the tail region as a higher estimate of the energy budget. Using these dimensions and the CR energy density, we found that the energy budget for these CR protons that are > 560 TeV is 7.62 × 10⁴⁴ erg.

The hadronic-type modeling in this section assumes that protons are producing the UHE HAWC data. After being accelerated in the shock fronts of the SNR ejecta, the protons travel to the nearby molecular cloud and produce gamma rays through PD. We get a proton energy budget of 7.62 × 10⁴⁴ erg for protons > 560 TeV using the results from the molecular cloud template fit. This is well below the SNR’s lower-end energy budget estimation of 7 × 10⁴⁹ ergs (Kothes et al. 2001), so the SNR is fully capable of producing the UHE gamma rays we are seeing from HAWC. Following the logic from Amenomori et al. (2021), this SNR, which is capable of accelerating protons up to at least ≈0.5 PeV, should also be capable of producing protons up to ≈1.6 PeV during its free expansion phase at 1 kyr. Since the SNR age is believed to be 10 kyr, the diffusion of protons out to the molecular cloud may have been suppressed until now.

4.2. Multi-wavelength modeling for electrons

We used the results from Section 3.2.2 to explore the possible leptonic nature of the UHE HAWC emission through multi-wavelength modeling. Since we assumed that the hadronic nature of the UHE gamma ray emission stems from the shock fronts in the tail of the SNR, we assumed the same scenario for the leptonic emission we model here. Therefore, we used other gamma ray observations that also modeled the tail region for additional information during our modeling. We used flux points from HAWC, MAGIC’s tail source (MAGIC Collaboration 2023), and Fermi-LAT (Fang et al. 2022) to model the gamma ray production using Naima (Zabalza 2015). The X-ray observation is from the Suzaku telescope (Ge et al. 2021) and its flux was calculated using the “East” region and a total solid angle of 17.8² arcmin². The radio observations come from the DRAO at 408 MHz and 1240 MHz (Kothes et al. 2001), the Sino-German λ6 cm polarization survey (Gao et al. 2011), and the Effelsberg λ11 and λ21 cm surveys (Gao et al. 2011). Since the radio surveys for the Sino-German and Effelsberg surveys are reported for the entire SNR, their fluxes were scaled using the DRAO 408 MHz tail-to-whole flux ratio.

The multi-wavelength modeling was done using the Naima framework (Zabalza 2015). All parent particle spectral models used a power-law with an exponential cutoff,

$$\begin{array}{c}\hfill f(E)=A{(E/E_0)}^{-\alpha }exp(-{(E/E_{\mathrm{cutoff}})}^{\beta })\end{array}$$(9)

where A is the normalization parameter, E₀ is the pivot energy, α is the spectral index, E_cutoff is the cutoff energy, and β is the curvature parameter, set to the default value, 1.0, for all modeling.

The gamma ray observations are assumed to be produced through IC. For IC, the target photon fields considered are the cosmic microwave background (CMB) and near-infrared (NIR) (same as that used in Amenomori et al. (2021)). The radio and X-ray observations are assumed to be produced through nonthermal synchrotron radiation. This modeling also assumes a single population of electrons that are responsible for the radio, X-ray, and gamma ray observations. The best-fit values from the leptonic fit can be found in Table 3 and Figure 5 shows the multi-wavelength model.

The leptonic modeling in this section assumed electrons are producing the UHE HAWC data which, after being accelerated in the shock fronts of the SNR ejecta or the PWN’s wind, interact with the nearby photon fields and undergo IC scattering. The synchrotron and IC emission can explain the emission being seen in this region well (Figure 5). Using this multi-wavelength model, we got a total electron energy budget of 1.28 × 10⁴⁵ erg for electrons > 1 TeV, which is below the energy budget estimation of the SNR.

Fig. 5. Top: multi-wavelength SED. The data shown includes the DRAO flux points in red squares (Pineault & Joncas 2000), the Sino-German and Effelsberg measurements in dark green circles (Gao et al. 2011), the Suzaku spectrum in yellow (Ge et al. 2021), the Fermi-LAT flux points in black dots (Fang et al. 2022), the MAGIC flux points in green diamonds (MAGIC Collaboration 2023), and the HAWC flux points from the VHE work in blue triangles. Synchrotron is shown in dashed light grey, IC scattering is shown in dotted dark grey, and black shows the total SED from synchrotron and IC scattering combined. Bottom: offset between the data points used and the model determined with Naima. The Fermi-LAT upper limits are removed to better see the deviation about zero for the other data points.

Breuhaus et al. (2021) suggests that electron cooling dominated by Inverse Compton (IC) losses could produce UHE emission, resulting in a PeVatron. One requirement for this system is a pulsar with a spin-down energy of L_SD > 10³⁶ erg/s. While the PWN in this region fits that requirement (L_SD ≈ 2 × 10³⁷ erg/s) and would be an ideal electron accelerator, it is offset from the location of the UHE HAWC emission and located in the head of the SNR, which we are not assuming the UHE gamma ray emission comes from. This would mean that electrons are accelerated in the shock fronts of the SN ejecta and produce gamma rays through IC.

As mentioned in Section 1, MAGIC also saw a shift in emission towards the tail of the SNR as energy thresholds increased (MAGIC Collaboration 2023). Even though HAWC’s angular resolution is not as small as MAGIC’s, we do see our UHE emission trend towards the tail of the SNR (Figure 4b). Since we are only using energy as our differentiator for the head and tail regions of the SNR, it is also possible that electrons accelerated from the SN ejecta are producing most of the gamma ray emission, but some “leakage” from the PWN is also being accounted for in the UHE gamma ray emission.

5. Conclusions and outlook

The morphological studies we performed reveal a single extended source detected over 4–145 TeV. This modeling is consistent with other gamma ray observatories. As previously stated, any discrepancies between observations are most likely attributed to the varying parts of the SNR emitting at different energies in the gamma ray regime.

We also looked at HAWC’s > 56 TeV data to explore its correlation with nearby molecular clouds, which would hint at a hadronic gamma ray production mechanism. The molecular cloud template fits the region well and shows that protons could account for the gamma ray flux detected at these energies. We also explored the possibility of a leptonic origin for HAWC’s > 56 TeV data. While the multi-wavelength model fits all the data well, it is still hard to say what part of the region would be responsible for the UHE electrons producing gamma rays. If protons are the sole source of the UHE gamma ray emission, then SNR G106.3+2.7 is indeed a hadronic PeVatron. However, we cannot rule out the possibility of electrons producing the gamma ray emission at this time.

Currently, VERITAS, Tibet ASγ, and LHAASO also report gamma ray emission from this SNR but cannot resolve the head and tail regions of the SNR at this time. More VHE gamma ray observations disentangling the head and tail of the SNR would greatly help in describing the gamma ray emission in this region. Particularly, it would help in narrowing down which, if not both, of the regions of the SNR are capable of being a PeVatron.

Acknowledgments

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, CF-2023-I-645, cátedras 873, 1563, 341, 323, Red HAWC, México; DGAPA-UNAM grants IG101323, IN111716-3, IN111419, IA102019, IN106521, IN114924, IN110521, IN102223; 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, grant CIDEGENT/2018/034; The Program Management Unit for Human Resources & Institutional Development, Research and Innovation, NXPO (grant number B16F630069); Coordinación General Académica e Innovación (CGAI-UdeG), PRODEP-SEP UDG-CA-499; Institute of Cosmic Ray Research (ICRR), University of Tokyo. H.F. acknowledges support by NASA under award number 80GSFC21M0002; National Research Foundation of Korea (RS-2023-00280210). We also acknowledge the significant contributions over many years of Stefan Westerhoff, Gaurang Yodh and Arnulfo Zepeda Domínguez, all deceased members of the HAWC collaboration. Thanks to Scott Delay, Luciano Díaz and Eduardo Murrieta for technical support. This research was also supported by the Nicholas Matwiyoff & Carl Hogberg Endowed Graduate Fellowship at Michigan Technological University.

References

Appendix A: Validating the Planck 353 GHz dust opacity map

Two different surveys were used to assess the molecular clouds in this region: the Dame CO survey (Dame et al. 2001) and the 353 GHz Planck dust opacity map (Ade et al. 2011). We used the Dame CO survey to ensure that there is a molecular cloud at the distance of the Boomerang region. Figure A1 shows the velocity-integrated CO map on top, with a black box indicating the Boomerang region, and the velocity range for a given latitude on bottom¹. From the top image, we can see that there is emission in the region and the bottom shows that there isn’t any CO emission in front, or behind, the Boomerang region.

Fig. A.1. Dame CO survey for the Boomerang region (Dame et al. 2001). Top: Velocity-integrated CO map. The black box indicates the Boomerang region. Bottom: Velocity range over a given longitude. The black box indicates the Boomerang longitude.

From here, we plotted the velocity distribution at the location of the Boomerang to see what velocity range should be integrated over in the CO map. Figure A2 shows that about −20 km/s to 0 km/s is the ideal range. Kothes et al. (2001) found a peak of CO emission at −6.4 km/s, which corresponds to a distance 800 pc. We see a similar peak of emission around approximately −6 km/s, which tells us we are seeing a similar CO emission and can use a distance of 800 pc for the SNR and molecular cloud distances. The right shows the integrated CO map for the same velocity range. It can be seen that there is CO emission coincident with the HAWC VHE emission.

The 353 GHz Planck dust opacity map is newer than the Dame CO survey and provides information on HI emission, as well as CO emission because it shows gas that is optically thick in CO. Figure A1 shows that there isn’t any CO emission in front, or behind, the Boomerang region. This is important because the Planck data does not take into account distance, so it is looking at everything in front, and behind, the intended target. This means that the emission we are seeing in the Planck data is similar to the emission we are seeing in the Dame CO map at the distance of the Boomerang region and can be used for the template analysis.

Fig. A.2. Velocity distribution of molecular hydrogen for a longitude and latitude of 106.4° ± 0.5 and 2.9° ± 0.5, respectively

Appendix B: Fitting details for the Planck template

SNR G106.3+2.7 is a highly asymmetric source, which is not well captured by the symmetrical source assumptions that 3ML makes. We instead chose to do a template fit for the hadronic scenario to hopefully model the asymmetric source more appropriately. The template is also physically justified by the general assumption that if gamma rays are being produced through nearby molecular clouds, then the gamma ray emission might take the shape of the molecular cloud that is helping produce them. However, the Planck data has finer angular resolutions than HAWC’s point-spread function (PSF), so when the template is convolved with HAWC’s PSF in the 3ML fit, we lose the fine details of the molecular clouds (figure B1).

The size of the molecular cloud template is based on the size of HAWC’s VHE emission, as well as the distance that would be reasonable for protons to travel for interaction. The template is 1° × 1.5°, which corresponds to about 14 pc × 21 pc and appropriately accounts for HAWC’s emission size. Protons in this region should be able to diffuse out to about 60 pc based on Albert et al. (2020), so the chosen size for the molecular cloud template is also still appropriate for proton diffusion. Since the flux fit from this template and the mass used in the CR energy density calculation is based on the size of the template we choose for the fit, the results shown here are dependent on our choices made for the template size.

Fig. B.1. a) The Planck 353 GHz template normalized to 1/sr that is used in the VHE HAWC fit. b) The Planck model template after being convolved with HAWC’s PSF.

All Tables

All Figures

	Fig. 1. Diagram of the multi-source search method described in Section 3. Yellow boxes (first two columns) describe how point sources (PSs) are added, where PL refers to a power-law spectrum. Orange/red boxes (last four columns) describe how extended sources are added.
In the text

Fig. 2. Significance maps for HAWC’s all-energy emission and its comparison to the SNR in the region. a) The HAWC all-energy significance map of the region with the 4, 6, and 8σ contours overlaid in green. b) The HAWC significance map of the region with labels showing the best-fit positions from this analysis, the previous HAWC publication (Albert et al. 2020), and VERITAS (Acciari et al. 2009). The green dashed, white dotted, cyan dot-dashed, and blue solid contour indicate LHAASO’s analysis region (Cao et al. 2021), Tibet ASγ’s 1σ extension (Amenomori et al. 2021), Fermi-LAT’s 1σ extension (Fang et al. 2022), and MAGIC’s analysis regions (MAGIC Collaboration 2023), respectively. c) The brightness temperature image (Taylor et al. 2003) of the SNR overlaid with the HAWC 4, 6, and 8σ contours in white. The bright emission on the left shows the head of the SNR, while the faint emission extending to the right shows the tail.

In the text

Fig. 3. Spectral energy distribution (SED) of the ground-based VHE gamma ray observations. The green SED with circle flux points are the all-energy HAWC results from this work. The magenta dotted line is the HAWC result for the > 56 TeV molecular cloud template analysis. The red line (KM2A) and blue line (WCDA) are from the first LHAASO catalog (Cao et al. 2024). All SEDs show only statistical uncertainties. The open stars, right-facing triangles, left-facing triangles, squares, and x-marks correspond to the VERITAS (Acciari et al. 2009), MAGIC tail (MAGIC Collaboration 2023), MAGIC head (MAGIC Collaboration 2023), Tibet ASγ (Amenomori et al. 2021), and LHAASO (Cao et al. 2021) measurements, respectively.

In the text

Fig. 4. Significance maps for HAWC’s > 56 TeV emission and its comparison to the SNR and molecular clouds in the region. a) The HAWC significance map of the region above 56 TeV. b) The brightness temperature image (Taylor et al. 2003) of the SNR overlaid with the HAWC 3, 4, and 5σ contours in white. c) Molecular hydrogen column density integrated over a velocity range of −20 km/s to 0 km/s (Dame et al. 2001). d) The Planck 353 GHz template normalized to 1/sr that is used in the UHE HAWC fit.

In the text

Fig. 5. Top: multi-wavelength SED. The data shown includes the DRAO flux points in red squares (Pineault & Joncas 2000), the Sino-German and Effelsberg measurements in dark green circles (Gao et al. 2011), the Suzaku spectrum in yellow (Ge et al. 2021), the Fermi-LAT flux points in black dots (Fang et al. 2022), the MAGIC flux points in green diamonds (MAGIC Collaboration 2023), and the HAWC flux points from the VHE work in blue triangles. Synchrotron is shown in dashed light grey, IC scattering is shown in dotted dark grey, and black shows the total SED from synchrotron and IC scattering combined. Bottom: offset between the data points used and the model determined with Naima. The Fermi-LAT upper limits are removed to better see the deviation about zero for the other data points.

In the text

	Fig. A.1. Dame CO survey for the Boomerang region (Dame et al. 2001). Top: Velocity-integrated CO map. The black box indicates the Boomerang region. Bottom: Velocity range over a given longitude. The black box indicates the Boomerang longitude.
In the text

	Fig. A.2. Velocity distribution of molecular hydrogen for a longitude and latitude of 106.4° ± 0.5 and 2.9° ± 0.5, respectively
In the text

	Fig. B.1. a) The Planck 353 GHz template normalized to 1/sr that is used in the VHE HAWC fit. b) The Planck model template after being convolved with HAWC’s PSF.
In the text