Context. Supernova remnants SNRs are believed to be capable of accelerating cosmic rays CRs to PeV energies. SNR G106.32.7 is a prime PeVatron candidate. It is formed by a head region where the pulsar J22296114 and its boomerangshaped pulsar wind nebula are located and a tail region containing SN ejecta. The lack of observed gamma ray emission from the two regions of this SNR has made it difficult to assess which region would be responsible for the PeV CRs. Aims. We aim to characterize the veryhighenergy VHE 0.1100 TeV gamma ray emission from SNR G106.32.7 by determining the morphology and spectral energy distribution of the region. This is accomplished using 2565 days of data and improved reconstruction algorithms from the High Altitude Water Cherenkov HAWC Observatory. We also explore possible gamma ray production mechanisms for different energy ranges. Methods. Using a multisource fitting procedure based on a maximumlikelihood estimation method we evaluate the complex nature of this region. We determine the morphology spectrum and energy range for the source found in the region. Molecular cloud information is also used to create a template and evaluate the HAWC gamma ray spectral properties at ultrahighenergies UHE gt 56 TeV. This will help probe the hadronic nature of the highestenergy emission from the region. Results. We resolve one extended source coincident with all other gamma ray observations of the region. The emission reaches above 100 TeV and its preferred logparabola shape in the spectrum shows a flux peak in the TeV range. The molecular cloud template fit on the higher energy data reveals that the SNRs energy budget is fully capable of producing a purely hadronic source for UHE gamma rays. Conclusions. The HAWC observatory resolves one extended source between the head and the tail of SNR G106.32.7 in the VHE gamma ray regime. The template fit suggests the highest energy gamma rays could come from a hadronic origin. However the leptonic scenario or a combination of the two cannot be excluded at this time.
['astroparticle physics', 'pulsars: general', 'ISM: supernova remnants', 'gamma rays: general', 'Astrophysics - High Energy Astrophysical Phenomena']
Testing the molecular cloud paradigm for ultrahighenergy gamma ray emission from the direction of SNR G106.32.7
{'Testing the molecular cloud paradigm for ultra-high-energy gamma ray emission from the direction of SNR G106.3+2.7': 'R. Alfaro 1 , C. Alvarez 2 , J.C. Arteaga-Velázquez 3 , D. Avila Rojas 1 , H.A. Ayala Solares 4 , R. Babu 5 , E. Belmont-Moreno 1 , A. Bernal 6 , K.S. Caballero-Mora 2 , T. Capistrán 6 , A. Carramiñana 7 , S. Casanova 8 , U. Cotti 3 , J. Cotzomi 9 , S. Coutiño de León 10 , E. De la Fuente 11 , C. de León 3 , D. Depaoli 12 , P. Desiati 10 , N. Di Lalla 13 , R. Diaz Hernandez 7 , B.L. Dingus 14 , M.A. DuVernois 10 , K. Engel 15 , T. Ergin 5 , C. Espinoza 1 , K.L. Fan 15 , K. Fang 10 , N. Fraija 6 , S. Fraija 6 , J.A. García-González 16 , F. Garfias 6 , M.M. González 6 , J.A. Goodman 15 , S. Groetsch 17 , J.P. Harding 14 , S. Hernández-Cadena 26 , I. Herzog 5 , J. Hinton 12 , D. Huang 15 , F. Hueyotl-Zahuantitla 2 , T.B. Humensky 28 , P. Hüntemeyer 17 , S. Kaufmann 18 , D. Kieda 25 , W.H. Lee 6 , J. Lee 19 , H. León Vargas 1 , J.T. Linnemann 5 , A.L. Longinotti 6 , G. Luis-Raya 18 , K. Malone 14 , O. Martinez 9 , J. Martínez-Castro 20 , J.A. Matthews 21 , P. Miranda-Romagnoli 22 , J.A. Montes 6 , E. Moreno 9 , M. Mostafá 27 , M. Najafi 17 , L. Nellen 23 , M.U. Nisa 5 , L. Olivera-Nieto 12 , N. Omodei 13 , Y. Pérez Araujo 1 , E.G. Pérez-Pérez 18 , C.D. Rho 24 , D. Rosa-González 7 , H. Salazar 9 , D. Salazar-Gallegos 5 , A. Sandoval 1 , M. Schneider 15 , J. Serna-Franco 1 , A.J. Smith 15 , Y. Son 19 , R.W. Springer 25 , O. Tibolla 18 , K. Tollefson 5 , I. Torres 7 , R. Torres-Escobedo 26 , R. Turner 17 *, F. Ureña-Mena 7 , E. Varela 9 , L. Villaseñor 9 , X. Wang 17 , Z. Wang 15 , I.J. Watson 19 , E. Willox 15 , S. Yu 4 , S. Yun-Cárcamo 15 , and H. Zhou 26  \n(A ffi liations can be found after the references)  \nReceived July 15, 2024; Accepted Sept. 9, 2024', 'ABSTRACT': "Context. Supernova remnants (SNRs) are believed to be capable of accelerating cosmic rays (CRs) to PeV energies. SNR G106.3 + 2.7 is a prime PeVatron candidate. It is formed by a head region, where the pulsar J2229 + 6114 and its boomerang-shaped pulsar wind nebula are located, and a tail region containing SN ejecta. The lack of observed gamma ray emission from the two regions of this SNR has made it di ffi cult to assess which region would be responsible for the PeV CRs.  \nAims. Weaim to characterize the very-high-energy (VHE, 0.1-100 TeV) gamma ray emission from SNR G106.3 + 2.7 by determining the morphology and spectral energy distribution of the region. This is accomplished using 2565 days of data and improved reconstruction algorithms from the High Altitude Water Cherenkov (HAWC) Observatory. We also explore possible gamma ray production mechanisms for di ff erent energy ranges.  \nMethods. Using a multi-source fitting procedure based on a maximum-likelihood estimation method, we evaluate the complex nature of this region. We determine the morphology, spectrum, and energy range for the source found in the region. Molecular cloud information is also used to create a template and evaluate the HAWC gamma ray spectral properties at ultra-high-energies (UHE, > 56 TeV). This will help probe the hadronic nature of the highest-energy emission from the region.  \nResults. We resolve one extended source coincident with all other gamma ray observations of the region. The emission reaches above 100 TeV and its preferred log-parabola shape in the spectrum shows a flux peak in the TeV range. The molecular cloud template fit on the higher energy data reveals that the SNR's energy budget is fully capable of producing a purely hadronic source for UHE gamma rays.  \nConclusions. The HAWC observatory resolves one extended source between the head and the tail of SNR G106.3 + 2.7 in the VHE gamma ray regime. The template fit suggests the highest energy gamma rays could come from a hadronic origin. However, the leptonic scenario, or a combination of the two, cannot be excluded at this time.  \nKey words. gamma ray - Cosmic-ray - Particle Acceleration - Supernova Remnant - Molecular Clouds - Pulsar Wind Nebula", '1. Introduction': 'Cosmic-ray (CR) accelerators have been probed indirectly across a wide range of energies, from the radio to ultra-highenergy 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).  \nSNR 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 pul-  \nsar 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.  \nX-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.  \nGround-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 (Abe et al. 2022) report gamma ray  \nFig. 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.  \n<!-- image -->  \nemission 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.  \nAbe et al. (2022) 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.  \nThe 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 ∼ 22000 m 2 (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  \nthe 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.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 multisource 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:  \nTS = 2 ln L alt L null , (1)  \nwhere 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 di ff erence 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.  \nPoint 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.  \n- 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.\n- 2. After fitting the model in step 1, we check the ∆ TS between the previous model and the most recently tested model:\n- (a) If the ∆ TS > 25, we keep the source in the model and return to step 1.\n- (b) 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.  \n<!-- image -->  \n(a)  \n<!-- image -->  \nFig. 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 (Abe et al. 2022), 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.  \n<!-- image -->  \nExtended 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.  \n- 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.\n- 2. After fitting the model in step 1, we check the ∆ TS between the previous model and the most recently tested model:\n- (a) If the ∆ TS > 16, we accept the extended source model and go to step 3.\n- (b) 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.\n- 3. Next, we check the TS values of the point sources in the model:\n- (a) 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.\n- (b) 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.  \nOnce 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:  \ndN dE = N 0 GLYPH<18> E 30TeV GLYPH<19> α -β log ( E 30TeV ) , (2)  \nwhere N 0 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.  \nAfter 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 e ff ects on the data resulting from di ff erences in the modeling of the instrument. These di ff erences are calculated following Abeysekara et al. (2019).  \nThe 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  \nFig. 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 (Abe et al. 2022), MAGIC head (Abe et al. 2022), Tibet AS γ (Amenomori et al. 2021), and LHAASO (Cao et al. 2021) measurements, respectively.  \n<!-- image -->  \nnormalization parameter N 0 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).  \nThe position (figure 2.b) 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 ± 0 . 1( stat ) ± 0 . 1( sys ) degrees (shown in figure 2.b) and LHAASO uses an extension of 0.3 degrees in their highenergy source catalog (shown in figure 2.b), 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 di ff erences between the various gamma ray observations are most likely attributed to each observatory detecting di ff erent contributions from the head and tail regions across varying energy ranges.", '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  \nthere 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.  \nWe used two di ff erent 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 4.c 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.  \nAgain, 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:  \n- 1. We start by calculating the column density- this calculation is taken from Albert et al. (2021):  \nNH = τ D / GLYPH<18> τ D NH GLYPH<19> re f , (3)  \nwhere ( τ D / NH ) re f = (1.18 ± 0.17) × 10 -26 cm 2 for the 353 GHz Planck map (Ade et al. 2011) and τ D is the dust opacity.  \n- 2. We then crop the data to the size of the region of interest.\n- 3. We calculate the mass of the region using the following equation:  \nMdust = NH Ω d 2 mh , (4)  \nwhere Ω is the angular area of the cloud and d = 800 pc as the distance to the cloud (Kothes et al. 2001)  \n- 4. Finally, we normalize the data to 1 / sr. This is required for a 3ML 2D template.  \n.  \nSince 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 4.a. 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 dropo ff at higher energies (figure 3). The best-fit normalization at 50 TeV for the template fit is 2.8 + 0 . 9 -0 7 × 10 -16 (1 / TeV s cm 2 ).", '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  \nTable 1. Best-fit results from the source search. The first uncertainties listed are statistical and the second uncertainties are systematic. The spectral values are quoted assuming the spectral parameters as defined in equation 2.  \nFig. 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.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nbeing 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.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:  \nJ = Z Ef Ei E dN dE dE , (5)  \nwhere  \ndN dE = k GLYPH<18> E E piv GLYPH<19> -α , (6)  \nand Ei = 56 TeV, Ef = 316 TeV, and k , Epiv 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:  \nL ( ≥ E γ ) = 4 π d 2 J , (7)  \nwhere d = 800 pc. Finally, we used the luminosity to calculate the CR proton energy density using the same approach as Abramowski et al. (2016):  \nω CR ( ≥ 10 E γ ) ≈ 1 . 8 × 10 -2 GLYPH<18> η N 1 . 5 GLYPH<19> -1 GLYPH<18> L ( ≥ E γ ) 10 34 erg / s GLYPH<19>GLYPH<18> M 10 6 M ⊙ GLYPH<19> -1 , (8)  \nArticle number, page 5 of 11  \nTable 2. Best-fit results from the simple modeling above 56 TeV, which assumes a point source with a PL spectrum. All errors are statistical.  \n.  \n-  \n0  \nwhere η N = 1.5 for heavier nuclei, M = 0 . 23 × 10 5 M ⊙ and ω CR is in eV / cm 3 , 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 -3 eV / cm 3 for protons > 560 TeV.  \nThe 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 44 erg.  \nThe 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 44 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 49 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 di ff usion 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 multiwavelength 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 (Abe et al. 2022), 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 2 arcmin 2 . 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 E ff elsberg λ 11 and λ 21 cm surveys (Gao et al. 2011). Since the radio surveys for the Sino-German and E ff elsberg surveys are reported for the entire SNR, their fluxes were scaled using the DRAO 408 MHz tail-towhole flux ratio.  \nThe multi-wavelength modeling was done using the Naima framework (Zabalza 2015). All parent particle spectral models used a power-law with an exponential cuto ff ,  \n4  \nFig. 5. Top: Multi-wavelength SED. The data shown includes the DRAO flux points in red squares (Pineault & Joncas 2000), the SinoGerman and E ff elsberg 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 (Abe et al. 2022), 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: O ff set 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.  \n<!-- image -->  \nf ( E ) = A ( E / E 0) -α exp( -( E / E cuto ff ) β ) (9)  \nwhere A is the normalization parameter, E 0 is the pivot energy, α is the spectral index, Ecutof f is the cuto ff energy, and β is the curvature parameter, set to the default value, 1.0, for all modeling.  \nThe 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 Xray 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.  \nThe 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.  \nTable 3. Best-fit results from Naima modeling, where W e is the electron energy budget. All errors are statistical.  \nThe 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 45 erg for electrons > 1 TeV, which is below the energy budget estimation of the SNR.  \nBreuhaus 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 LSD > 10 36 erg / s. While the PWN in this region fits that requirement ( LSD ≈ 2 × 10 37 erg / s) and would be an ideal electron accelerator, it is o ff set 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.  \nAs mentioned in section 1, MAGIC also saw a shift in emission towards the tail of the SNR as energy thresholds increased (Abe et al. 2022). 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 4.b). Since we are only using energy as our di ff erentiator 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 di ff erent energies in the gamma ray regime.  \nWe 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.  \nCurrently, 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.  \nAcknowledgements. We acknowledge the support from: the US National Science Foundation (NSF); the US Department of Energy O ffi ce of High-Energy  \n-  \n115  \nPhysics; 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-202300280210). We also acknowledge the significant contributions over many years of Stefan Westerho ff , 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 Matwiyo ff & Carl Hogberg Endowed Graduate Fellowship at Michigan Technological University.", 'References': 'Abe, H., Abe, S., Acciari, V., et al. 2022, Astronomy and Astrophysics  \n- Abeysekara, A., Albert, A., Alfaro, R., et al. 2023, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 1052, 168253\n- Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017, The Astrophysical Journal, 843, 39\n- Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2019, American Astronomical Society, 881, 134\n- Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2022, in 37th International Cosmic Ray Conference.\n- Abramowski, A., Aharonian, F., Benkhali, F. A., et al. 2016, Nature, 531, 476-479\n- Acciari, V. A., Aliu, E., Arlen, T., et al. 2009, The Astrophysical Journal, 703, L6\n- Ackermann, M., Ajello, M., Baldini, L., et al. 2017, The Astrophysical Journal, 843, 139\n- Ade, P. A. R., Aghanim, N., Arnaud, M., et al. 2011, Astronomy and Astrophysics, 536, A1\n- Albert, A. ., Alfaro, R., Alvarez, C., et al. 2024, Performance of the HAWC Observatory and TeV Gamma-Ray Measurements of the Crab Nebula with Improved Extensive Air Shower Reconstruction Algorithms\n- Albert, A., Alfaro, R., Alvarez, C., et al. 2023, in 38th International Cosmic Ray Conference.\n- Albert, A., Alfaro, R., Alvarez, C., et al. 2021, The Astrophysical Journal, 914, 106\n- Albert, A., Alfaro, R., Alvarez, C., et al. 2020, The Astrophysical Journal Letters, 896  \nAmenomori, M., Bao, Y. W., Bi, X. J., et al. 2021, Nature Astronomy, 5, 460-464 Bell, A. R., Schure, K. M., Reville, B., & Giacinti, G. 2013, Monthly Notices of the Royal Astronomical Society, 431, 415-429  \n- Breuhaus, M., Hahn, J., Romoli, C., et al. 2021, The Astrophysical Journal Letters, 908, L49\n- Cao, Z., Aharonian, F., An, Q., et al. 2024, The Astrophysical Journal Supplement Series, 271, 25\n- Cao, Z., Aharonian, F., An, Q., et al. 2021, Nature, 594, 33-36  \nCristofari, P. 2021, Universe, 7, 324  \n- Cristofari, P., Gabici, S., Terrier, R., & Humensky, T. B. 2018, Monthly Notices of the Royal Astronomical Society, 479, 3415\n- Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, The Astrophysical Journal, 547, 792\n- Fang, K., Kerr, M., Blandford, R., Fleischhack, H., & Charles, E. 2022, Physical Review Letters, 129  \nFujita, Y., Bamba, A., Nobukawa, K. K., & Matsumoto, H. 2021, The Astrophysical Journal Letters, 912, 133  \n- Gao, X. Y., Han, J. L., Reich, W., et al. 2011, Astronomy and Astrophysics, 529 Ge, C., Liu, R.-Y., Niu, S., Chen, Y., & Wang, X.-Y. 2021, The Innovation, 2, 100118\n- Halpern, J. P., Gotthelf, E. V., Leighly, K. M., & Helfand, D. J. 2001, The Astrophysical Journal, 547, 323-333\n- Joncas, G. & Higgs, L. A. 1990, Astron. Astrophys. Suppl. Ser., 82, 113 Kothes, R., Uyaniker, B., & Pineault, S. 2001, The Astrophysical Journal Supplement Series, 560  \nPineault, S. & Joncas, G. 2000, The Astronomical Journal, 120, 3218 Reynolds, S. P. 2010, Astrophysics and Space Science, 336, 257-262 Taylor, A. R., Gibson, S. J., Peracaula, M., et al. 2003, AJ, 125, 3145  \n- Vannoni, G., Gabici, S., & Aharonian, F. A. 2009, Astronomy & Astrophysics, 497, 17-26\n- Vianello, G., Lauer, R. J., Younk, P., et al. 2015, The Multi-Mission Maximum Likelihood framework (3ML)  \nZabalza, V. 2015, Proceedings of Science  \n- 1 Instituto de Física, Universidad Nacional Autónoma de México, Ciudad de Mexico, Mexico\n- 2 Universidad Autónoma de Chiapas, Tuxtla Gutiérrez, Chiapas, México\n- 3 Universidad Michoacana de San Nicolás de Hidalgo, Morelia, Mexico\n- 4 Department of Physics, Pennsylvania State University, University Park, PA, USA\n- 5 Department of Physics and Astronomy, Michigan State University, East Lansing, MI, USA\n- 6 Instituto de Astronomía, Universidad Nacional Autónoma de México, Ciudad de Mexico, Mexico\n- 7 Instituto Nacional de Astrofísica, Óptica y Electrónica, Puebla, Mexico\n- 8 Institute of Nuclear Physics Polish Academy of Sciences, PL-31342 IFJ-PAN, Krakow, Poland\n- 9 Facultad de Ciencias Físico Matemáticas, Benemérita Universidad Autónoma de Puebla, Puebla, Mexico\n- 10 Department of Physics, University of Wisconsin-Madison, Madison, WI, USA\n- 11 Departamento de Física, Centro Universitario de Ciencias Exactase Ingenierias, Universidad de Guadalajara, Guadalajara, Mexico\n- 12 Max-Planck Institute for Nuclear Physics, 69117 Heidelberg, Germany\n- 13 Department of Physics, Stanford University: Stanford, CA 94305-4060, USA\n- 14 Los Alamos National Laboratory, Los Alamos, NM, USA\n- 15 Department of Physics, University of Maryland, College Park, MD, USA\n- 16 Tecnologico de Monterrey, Escuela de Ingeniería y Ciencias, Ave. Eugenio Garza Sada 2501, Monterrey, N.L., Mexico, 64849\n- 17 Department of Physics, Michigan Technological University, Houghton, MI, USA\n- 18 Universidad Politecnica de Pachuca, Pachuca, Hgo, Mexico\n- 19 University of Seoul, Seoul, Rep. of Korea\n- 20 Centro de Investigación en Computación, Instituto Politécnico Nacional, México City, México.\n- 21 Dept of Physics and Astronomy, University of New Mexico, Albuquerque, NM, USA\n- 22 Universidad Autónoma del Estado de Hidalgo, Pachuca, Mexico\n- 23 Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de Mexico, Ciudad de Mexico, Mexico\n- 24 Department of Physics, Sungkyunkwan University, Suwon 16419, South Korea\n- 25 Department of Physics and Astronomy, University of Utah, Salt Lake City, UT, USA\n- 26 Tsung-Dao Lee Institute & School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai, China\n- 27 Department of Physics, Temple University, Philadelphia, Pennsylvania, USA\n- 28 NASA Goddard Space Flight Center, Greenbelt, MD, USA', 'Appendix A: Validating the Planck 353 GHz dust opacity map': "Two di ff erent 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 A.1 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 1 . 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.  \nFrom 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 A.2 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.  \nThe 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 A.1 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.", '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 B.1).  \nThe 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 di ff use 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 di ff usion. 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.  \nFig. A.1. Dame CO survey for the Boomerang region (Dame et al. 2001). Top: Velocityintegrated CO map. The black box indicates the Boomerang region. Bottom: Velocity range over a given longitude. The black box indicates the Boomerang longitude.  \n<!-- image -->  \nFig. A.2. Velocity distribution of molecular hydrogen for a longitude and latitude of 106.4 · ± 0.5 and 2.9 · ± 0.5, respectively  \n<!-- image -->  \n/sr  \n<!-- image -->  \nFig. 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.  \n<!-- image -->"}