Study of the IC 443 region with the HAWC observatory

We use data collected from 2015 to 2023. The livetime is 2966 days using the $f_{\rm hit}$ bin scheme, while with the reconstructed energy scheme, it is 2769 days since we remove data from the time period where HAWC was under construction and not all tanks were operational.

We apply a gamma-hadron selection cut found in Albert et al. (2024), and estimate the background using the method described in section 4.3 of (Abeysekara et al., 2019). Once we have the estimated background, we proceed with the analysis, based on the maximum likelihood estimation method. For the analysis we use the Multi-Mission Maximum Likelihood (threeML) framework with the HAWC Accelerated Likelihood (HAL plugin) as described in Vianello et al. (2015); Abeysekara et al. (2022). ThreeML uses the likelihood ratio test to determine whether one model is better than another. The test statistic is defined as

where $L_{i,j}$ are the maximum likelihoods of models $i,j$. It is important to remark that we also use this statistic when generating sky maps, where the calculation is performed for each individual pixel.

The algorithm we use to determine the morphology of the region is described in Albert et al. (2023), with the difference that we now use the High-Energy Radiative Messengers (HERMES) diffuse emission template as the galactic diffuse emission (GDE) model (Dundovic et al., 2021). HERMES is a spatial propagation template based on the DRAGON cosmic-ray propagation code (Evoli et al., 2017). The model includes pionic gamma rays, produced by the interaction of the sea of cosmic rays with neutral and molecular hydrogen. It also includes inverse Compton from the cosmic-ray leptons scattering on low-energy photons in the Galaxy such as UV, optical, infrared or microwave photons. This new GDE model has the advantage that it considers the morphology of the gas in the region, whereas Albert et al. (2023) used a Gaussian shape as a function of galactic latitude. When fitting the GDE model, we free a unit-less factor parameter that scales the normalization of the baseline diffuse flux calculated by HERMES.

To summarize the algorithm used in Albert et al. (2023), we first search for all possible point sources in the region and assume a simple power-law spectrum for each one as shown in Algorithm 1. Then, each point source is tested as an extended source as described in Algorithm 2. Once we have found the final morphology of each source, we test for curvature in the spectrum, with the process described in Algorithm 3.

The final step is to check whether the spectrum of the sources has a cutoff or curvature. We have tested for a logarithmic parabola and a power law with exponential cutoff, whose functional forms are:

where $\Phi_{SPL}$ is the same as in Equation 2 and $\beta$ is the curvature parameter;

where $\Phi_{SPL}$ is the same as in Equation 2 and $E_{\rm cutoff}$ is the cutoff energy.

Figure 1 shows the HAWC $f_{\rm hit}$ significance map of the region of IC443. The bright region shows the region of interest used for the analysis. The best-fit model describing the region together with the GDE model is a point source (PS) and an extended source (ES). We found no significant preference for a spectral curvature for either source. We associate the point source with IC 443. For the extended source, which we name HAWC J0615+2213, we did not find a possible counterpart in the literature. Therefore, we investigate two possibilities: emission originating from the interaction of freshly injected cosmic rays with the interstellar medium, or a new TeV Halo object. We defer this discussion to Sect. 4.

The position of the point source is located 0.26^∘ away from the Green catalog position of IC 443 and 0.14^∘ from the pulsar wind nebula associated with the neutron star CXOU J061705.3+222212. The centroid of the extended source is 0.63^∘ from IC 443, 0.57^∘ from CXOU J061705.3+222212, and 0.29^∘ from the pulsar B0611+22. We observe no significant excess in the direction of the supernova remnant G189.6+03.3. The source parameter results are summarized in Table 1. The pivot energy, $E_{\rm piv}$, that reduces the correlation between the parameters was found to be 2.3 TeV. Table 2 shows the comparison of the fitted spectral parameters between the $f_{\rm hit}$ binning scheme and the energy binning scheme.

Figure 2 shows the map of the best-fit model, as well as the residual map of the region, and the significance distribution of the sky. The significance distribution is relatively close to a Gaussian distribution, showing that the model adequately describes most of the signal and the remaining is only background. Figure 3 shows each source after subtracting either the extended source or the point source from the original map.

The energy ranges in which the sources are significantly detected are calculated using the method described in Abeysekara et al. (2017). At the 1$\sigma$, or 68% confidence level, these are [0.3,30] TeV for the point source and [0.7,47] TeV for the extended source.

As mentioned before, the 1D binning scheme is used to obtain the best spectral shape of a source with threeML. To obtain the spectral points, we use the 2D binning scheme. We fix all the parameters of the model except for the normalization of the source in each energy bin. Tables 4 and 4 present the differential particle flux measured by HAWC and the E²-weighted differential particle flux is shown in Fig. 4. The figure also compares the spectral model fitted with the $f_{\rm hit}$ scheme and the energy bin scheme, showing that they are consistent with each other.

Several researchers have performed extensive multi-wavelength studies of IC 443 over the years. Early works by Chevalier (1999) and Baring et al. (1999) laid the groundwork of multi-wavelength observations, focusing mainly on the radio and X-ray components; while more recent studies as Ackermann et al. (2013); Zhang et al. (2018) and Sharma et al. (2023), have furthered our understanding at higher energies up to 1 TeV. We use the new HAWC observations to extend the multi-wavelength analysis up to $30$ TeV. As mentioned earlier, we associate the point source we found in Sect. 3 with IC 443.

In the very-high-energy ($>$100 MeV) part of the electromagnetic spectrum, Ackermann et al. (2013) demonstrated that gamma-ray emission is primarily due to the decay of neutral pions produced in inelastic collisions between cosmic-ray protons and dense molecular clouds in the region (Seta et al., 1998; Ustamujic et al., 2021). They also explored leptonic models, showing that these models could not account for the observed gamma-ray flux due to energetic constraints (see Abdo et al., 2009). These constraints arise from the maximum energy that electrons can achieve in an environment like IC 443. Electron acceleration is restricted by the underlying mechanism, such as diffusive shock acceleration (DSA), and by substantial energy losses due to synchrotron radiation, inverse Compton scattering, and bremsstrahlung. These losses are further amplified by the dense molecular cloud surrounding IC 443.

Under typical magnetic fields ($\sim$100 $\mu$G) and accounting for both synchrotron and inverse Compton losses, the maximum energy an electron can gain is approximately 5 TeV (Petrosian, 2012; Yamaguchi et al., 2009). When bremsstrahlung losses are included in regions with higher gas densities (up to 100 cm^-3) (Troja et al., 2008), this maximum decreases further to around 3 TeV. Even with non-linear shock effects that amplify the magnetic field up to  300 $\mu$G, the maximum achievable electron energy is only $\sim$10 TeV. Thus, while electrons can contribute to high-energy emission, they cannot explain gamma rays up to 30 TeV. This limitation makes hadronic processes a more plausible explanation for IC 443’s observed gamma-ray spectrum (Guo et al., 2014; Zhang & Yan, 2011).

Protons can achieve significantly higher energies due to their relatively low energy losses. The maximum energy for protons, achievable via DSA, depends on factors such as shock speed, magnetic field strength, and remnant size. For IC 443, typical values include a shock speed of 300 km s^-1 (Lee et al., 2008; Troja et al., 2008), a magnetic field strength of 10 $\mu$G (Yamaguchi et al., 2009), and a remnant size of 10.9 pc (Troja et al., 2008). These parameters suggest that protons could theoretically reach energies up to $\sim$65 TeV (Drury, 1983; Lagage & Cesarsky, 1983; Bell et al., 2013). With the inclusion of non-linear diffusive shock acceleration (NLDSA) effects, which amplify the magnetic field, this maximum energy may extend to 100 TeV or even higher, potentially approaching 1 PeV (Bell, 2004; Bell et al., 2013; Caprioli et al., 2009).

To test these estimates of the maximum proton energies, we perform an analysis using the Markov Chain Monte Carlo (MCMC) software Naima (Zabalza, 2015). This software fits data to non-thermal emission models by optimizing the parameters of the parent particle spectra. Additionally, it calculates the Bayesian Information Criterion (BIC) for the models, which we use to compare and evaluate their performance.

We incorporate observations from Fermi (Abdo et al., 2010), VERITAS (Acciari et al., 2009), MAGIC (Albert et al., 2007), and the results presented here. Additionally, the recently commissioned LHAASO detector—a gamma-ray observatory with better sensitivity above 10 TeV compared to HAWC—has reported extended emission in the region as part of its catalog search. LHAASO observes a source with a width of $(0.59\pm 0.08)^{\circ}$ located at coordinates 94.35^∘, 22.57^∘. Their analysis using LHAASO-WCDA data, based on a simple power-law model, yielded a differential particle flux at 3 TeV of (1.95$\pm$0.27)$\times 10^{-13}$ TeV^-1 cm^-2 s^-1 with a spectral index of $-2.92\pm 0.14$ in the energy range between 1 TeV to 25 TeV. An upper limit of 0.17 $\times 10^{-15}$ TeV^-1 cm^-2 s^-1 at 50 TeV was set using LHAASO-KM2A (Cao et al., 2024). However, further detailed analysis is required to determine whether this observation aligns with the morphology identified in our study. Therefore, we do not include the LHAASO data in our fitting (See appendix B where we test the LHAASO model with HAWC data). Based on the arguments outlined above, we use a pion decay model to explain the gamma-ray emission, similar to the approach in Ackermann et al. (2013). We set the maximum energy values to 65 TeV and 1 PeV for the proton distribution for the test. The proton spectrum is modeled with both a broken power law and a power law with an exponential cutoff:

where $A$ is the normalization of the spectrum at energy $E_{0}$; $\alpha_{1}$ and $\alpha_{2}$ are the spectral indices before and after the break energy $E_{\rm break}$, while $\alpha$ is the spectral index in the exponential cutoff model; $E_{\rm cutoff}$ is the cutoff energy. The fitted parameters results after running the aforementioned models are shown in Table 5.

The results are consistent when we compare between the different maximum energies in each of the cosmic-ray distributions. When comparing between spectral models, the broken power-law model provides a better fit than the exponential cutoff based on the BIC. We also fitted the density of the environment, which is in the same order of magnitude as shown in Ustamujic et al. (2021) and can impact the total energy of the protons required to produce gamma rays. These models yield a total proton energy content with $E\geq 0.8$ GeV of $W_{p}\sim[6-8]\times 10^{48}(n/200{\rm\,cm}^{-3})^{-1}$ erg, which is less than 1 per cent of the total supernova remnant energy of $\sim 10^{51}$ erg, consistent with Ackermann et al. (2013), even though they assumed a density of $20{\rm\,cm}^{-3}$. Figure 5 shows that the tail of the distribution is better described by extending the proton spectrum to 1 PeV.

We tested a model where we left the maximum energy of the proton distribution to be free. Since the HAWC data do not show any cutoff, the posterior distribution of the maximum energy remains relatively flat, as it is shown in Fig. 6, compared to a prior uniform distribution in logarithmic space in the range of 65 TeV to 10 PeV.

This suggests that IC 443 has the potential to accelerate cosmic rays to energies beyond the  65 TeV limit predicted by DSA, making it a strong candidate for a PeVatron. Given the constraints on leptonic processes, our findings further support IC 443 as a compelling source of PeV cosmic rays. However, higher-energy observations are still needed to determine where the gamma-ray spectrum cuts off, as this will provide a direct measure of the maximum energy achievable by protons.

The extended emission we found in the analysis is surprising. First, it does not spatially correlate with the other the supernova remnant G189.6, which is 1.2^∘ away from the centroid of the extended source. We see no significant emission from this other remnant, which means that it does not emit gamma rays at a level detectable by HAWC. After discarding this possibility, we investigate two other scenarios next.