HAWC Study of the Very-high-energy γ-Ray Spectrum of HAWC J1844−034

The HAWC observatory is a ground-based particle sampling array located on Sierra Negra, Mexico, at an altitude of 4100 m. It is designed to observe very-high-energy (VHE) γ-ray photons by detecting Cherenkov radiation produced from electromagnetic air showers as they propagate through the water Cherenkov detectors (WCDs). The main array of HAWC, composed of 300 WCDs, covers a geometrical area of 22,000 m² (Abeysekara et al. 2023). Each WCD consists of a bladder containing 180,000 L of purified water that is protected by a metal casing with a roof on top. A WCD has a height of 5.4 m and a diameter of 7.3 m. Four photomultiplier tubes (PMTs) are located at the bottom of each WCD to collect the Cherenkov light emitted by air shower particles. For more detailed information on the design and operation of HAWC, refer to Abeysekara et al. (2017a, 2023).

The γ-ray events collected by HAWC are categorized into analysis bins based on the fraction of the PMTs used in the event reconstruction (f_hit) and the estimated energies of the particular event (Abeysekara et al. 2017a, 2019). We call these the "ground parameter" (GP) and "neural network" (NN) energy estimators.

The HAWC High-Energy Catalog released in 2020 used the GP data set that contains 1039 days of cumulative data (Abeysekara et al. 2020). In this work, we have used 1910 days of the GP energy estimator data above 17.8 TeV and f_hit above 24.5% that are within our region of interest (ROI) to compute the most sensitive results on the high-energy spectrum of HAWC J1844−034. The energy resolution for the GP energy estimator is 50% at 20 TeV and 34% at 100 TeV (Abeysekara et al. 2019). The angular resolution range for the GP energy estimator bins used in the analysis is 0.1°–0.6° (Abeysekara et al. 2019). Our ROI is a rectangular region defined as 28° < l < 30° and −5° < b < 5°, where l and b are the Galactic longitude and latitude, respectively. The ROI is presented as the brighter region in Figure 1, which has been chosen to include the extended emission.

Zoom In Zoom Out Reset image size

The multimission maximum-likelihood ³⁰ framework and HAWC accelerated-likelihood ³¹ plug-in (Vianello et al. 2015; Abeysekara et al. 2021) are used to simultaneously fit multiple source models and their relevant free parameters based on maximum-likelihood estimations. Spectral and morphological models mentioned in this work can be found in astromodels. ³²

In our analysis, a test statistic (TS) is used to evaluate the pretrial statistical significance of a fitted model in a region with a predetermined number of free parameters. The TS provides a statistical measure of how well an alternative hypothesis performs over a null hypothesis. The TS values can be used to compare between different nested models. The TS is defined as

$$\begin{eqnarray}&&\mathrm{TS}=2\mathrm{ln}\left(\displaystyle \frac{{L}_{1}}{{L}_{0}}\right),\end{eqnarray} \tag{ 1 }$$

àwhere L₁ is the likelihood for an alternative hypothesis, and L₀ is the likelihood for a null hypothesis. To compare two different alternative hypotheses, ${\rm{\Delta }}\mathrm{TS}={\mathrm{TS}}_{2}-{\mathrm{TS}}_{1}\,=2\mathrm{ln}\left({L}_{2}/{L}_{1}\right)$ can be used, where L₂ is the likelihood for a different alternative hypothesis (Abeysekara et al. 2017a). Using Wilks's theorem (Wilks 1938), a pretrial significance can be computed with $\sigma \simeq \sqrt{\mathrm{TS}}$, which is used to generate HAWC significance maps (Abeysekara et al. 2017b).

In order to find a statistically optimal list of sources in the ROI, a systematic approach developed by the Fermi Large Area Telescope (Fermi-LAT) Collaboration is used (Ackermann et al. 2017). For the source search, we recursively add a point source to the model in our ROI until the improvement of the overall TS becomes less than 16. Then, we test an extended source model for each of the found point sources. Starting from the brightest source, we accept an extended source model if ΔTS > 16.

During this process, we assumed a power-law spectral model for all of the sources involved,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\gamma }\left({E}_{\gamma }\right)={{\rm{\Phi }}}_{0}{\left(\displaystyle \frac{{E}_{\gamma }}{{E}_{\mathrm{piv}}}\right)}^{-{\rm{\Gamma }}},\end{eqnarray} \tag{ 2 }$$

where ${{\rm{\Phi }}}_{\gamma }\left({E}_{\gamma }\right)$ is the γ-ray flux at the γ-ray energy E_γ , Φ₀ is the flux normalization at a pivot energy E_piv, and Γ is the spectral index. Here Φ₀ and Γ are set as free parameters to be determined by the likelihood fit, and E_piv is fixed at 20 TeV. For the pivot energy, we have scanned through the entire energy range from which we have obtained a minimum correlation between the free parameters at 20 TeV.

Once the optimal number of point and extended sources is determined from the systematic source search, more realistic spectral models are applied and tested. A power law with an exponential cutoff and log parabola models have been used:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\gamma }\left({E}_{\gamma }\right)={{\rm{\Phi }}}_{0}{\left(\displaystyle \frac{{E}_{\gamma }}{{E}_{\mathrm{piv}}}\right)}^{-{\rm{\Gamma }}}\exp \left(-\displaystyle \frac{{E}_{\gamma }}{{E}_{c}}\right),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\gamma }\left({E}_{\gamma }\right)={{\rm{\Phi }}}_{0}{\left(\displaystyle \frac{{E}_{\gamma }}{{E}_{\mathrm{piv}}}\right)}^{-{\rm{\Gamma }}-\beta \mathrm{ln}\left({E}_{\gamma }/{E}_{\mathrm{piv}}\right)},\end{eqnarray} \tag{ 4 }$$

where E_c of Equation (3) is the cutoff energy, and β of Equation (4) describes the curvature of the power-law spectrum.

The radially symmetric Gaussian morphology (Abeysekara et al. 2019) is formulated as

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{d{\rm{\Omega }}}={\left(\displaystyle \frac{180}{\pi }\right)}^{2}\displaystyle \frac{1}{2\pi {\sigma }^{2}}\exp \left(-\displaystyle \frac{{\theta }^{2}}{2{\sigma }^{2}}\right),\end{eqnarray} \tag{ 5 }$$

where θ and σ are the angular distance and Gaussian width as a measure of the size of the extended source model, respectively. Here σ is set as a free parameter during likelihood fits.

The Naima software ³³ (Zabalza 2015) is a nonthermal γ-ray modeling framework based on Markov Chain Monte Carlo methods. It is used to model and fit the spectra of accelerated particles in hadronic and leptonic scenarios based on the observed electromagnetic data. The Naima software provides numerical models for a hadronic scenario based on neutral pion decays from inelastic proton–proton collisions (Kafexhiu et al. 2014) and a leptonic scenario based on inverse Compton scattering of relativistic electrons with low-energy photons and synchrotron emission from relativistic electrons (Aharonian et al. 2010; Khangulyan et al. 2014). We have used the Naima framework to generate particle spectra for the hadronic and leptonic scenarios of HAWC J1844−034. The results are presented in Sections 4.4 and 5.2.

Figure 4 shows the spectrum of HAWC J1844−034. The gray and red bands represent the statistical and systematic uncertainties, respectively. The various systematic uncertainty components have been summed in quadrature (later discussed in Section 4.3). We have calculated the flux points of HAWC J1844−034 (black) by refitting Φ₀ at each energy bin while maintaining the best-fit parameters for the spectral shape, Γ and β (Abeysekara et al. 2019).

Zoom In Zoom Out Reset image size

In the figure, the TeV γ-ray flux data points from Tibet AS-Gamma (green) and LHAASO (purple) at 100 TeV have also been added for comparison (Cao et al. 2021; Amenomori et al. 2022). The fitted HAWC spectrum is in good agreement with the results from both LHAASO and Tibet AS-Gamma. According to the HAWC results, the falloff in the spectrum is gradual, which may be consistent with the 260 TeV photons observed by LHAASO (purple dashed line; Cao et al. 2021). The best-fit spectrum of TASG J1844−038 has been observed with a power-law function in 25 TeV < E_γ < 130 TeV (Amenomori et al. 2022). In this energy range, TASG J1844−038 has four flux points exceeding 2σ, which correspond well with the best-fit spectrum of HAWC J1844−034. The best-fit power-law spectrum for TASG J1844−038 has a spectral index of Γ = 3.26 ± 0.30, which falls more sharply than the index of HAWC J1844−034, ${\rm{\Gamma }}={2.72}_{-0.22}^{+0.21}$, at 20 TeV (Amenomori et al. 2022). This may be due to the different conditions we are performing the fits in, such as the energy range, pivot energy, and spectral model. However, there is sufficient overlap despite having different conditions, as shown in Figure 4. Also, at the pivot energy of TASG J1844−038 (E_piv = 40 TeV), the spectral index of HAWC J1844−034 is found to be ${\rm{\Gamma }}={3.07}_{-0.37}^{+0.36}$, which is consistent with the results from Tibet AS-Gamma. Furthermore, the sources found by the three γ-ray observatories are spatially consistent given their angular resolutions: ∼0.1°–0.6° for ≥10 TeV (HAWC), 0.24°–0.3° for 100 TeV (LHAASO), and ∼0.28° for 100 TeV (Tibet AS-Gamma; Abeysekara et al. 2019; Wang et al. 2021; Amenomori et al. 2022).

Figure 5 displays the radial distribution of the γ-ray excess profile from the best-fit location of HAWC J1844−034 (magenta points). The brown solid line represents the overall best-fit model that includes all three HAWC sources of the region. The black dashed, blue dotted–dashed, and red dotted lines represent the best-fit models for each source, and the dotted vertical lines represent the distances of each labeled source from HAWC J1844−034. Note that TASG J1844−038 (green) and LHAASO J1843−0338 (purple) are not included in the inclusive model. The inclusive model shows a good match with the excess. Also, TASG J1844−038 presents a good spatial agreement with HAWC J1844−034. Conversely, LHAASO J1843−0338 is positioned between the peaks of HAWC J1844−034 and HAWC J1843−032, which corresponds to a plateau of the overall excess peak.

Zoom In Zoom Out Reset image size

Factors that may contribute to the systematics are listed in Table 2. The first four components in the table arise from various detector and simulation effects (see Abeysekara et al. 2017a, 2019). The diffuse background emission (DBE)  is a systematic effect arising from any unresolved sources and the true Galactic diffuse emission. We have calculated this effect by adding a Gaussian background model along the Galactic equator with a fixed width of 1° (Rho et al. 2018; Albert et al. 2021b) in the nested fit. We have compared the best-fit parameters of the three sources and added the differences as a systematic component. A power-law spectral model (Equation (2)) with an index of 2.75 is assumed for the DBE based on Abdo et al. (2008) and Albert et al. (2022).

Table 2. Systematic Uncertainties

	HAWC J1844−034				HAWC J1843− 032		HAWC J1846−025
Uncertainty Origin	Φ0	Γ	β	Size	Φ0	Γ	Φ0	Γ
	[(TeV cm2 s)−1]			(deg)	[(TeV cm2 s)−1]		[(TeV cm2 s)−1]
Late light	${}_{-2.21}^{+1.08}\times {10}^{-15}$	${}_{-0.08}^{+0.22}$	${}_{-0.27}^{+0.10}$	${}_{-0.00}^{+0.00}$	${}_{-6.73}^{+3.87}\times {10}^{-17}$	${}_{-0.02}^{+0.01}$	${}_{-7.01}^{+3.75}\times {10}^{-17}$	${}_{-0.03}^{+0.08}$
Charge uncertainty	−7.48 × 10−16	+0.08	−0.04	+0.00	−4.75 × 10−17	+0.03	−2.26 × 10−17	+0.03
PMT efficiency	+4.67 × 10−16	${}_{-0.04}^{+0.02}$	−0.04	${}_{-0.00}^{+0.00}$	+6.13 × 10−17	${}_{-0.02}^{+0.01}$	+2.67 × 10−17	${}_{-0.03}^{+0.02}$
PMT threshold	${}_{-1.11}^{+1.36}\times {10}^{-16}$	${}_{-0.05}^{+0.09}$	${}_{-0.07}^{+0.06}$	${}_{-0.00}^{+0.00}$	${}_{-1.00}^{+1.86}\times {10}^{-18}$	${}_{-0.02}^{+0.04}$	+9.06 × 10−18	+0.03
DBE	−2.70 × 10−15	−0.11	+0.37	−0.08	−1.52 × 10−17	−0.02	−8.06 × 10−19	+0.06
ROI	+7.17 × 10−16	+0.00	−0.02	+0.03	−2.78 × 10−18	−0.02	+1.61 × 10−16	−0.34
Quadratic summed	${}_{-3.57}^{+1.38}\times {10}^{-15}$	${}_{-0.15}^{+0.25}$	${}_{-0.29}^{+0.39}$	${}_{-0.08}^{+0.03}$	${}_{-8.38}^{+7.25}\times {10}^{-17}$	${}_{-0.04}^{+0.05}$	${}_{-0.74}^{+1.68}\times {10}^{-16}$	${}_{-0.34}^{+0.11}$

Download table as:  ASCIITypeset image

Finally, our choice of ROI has been considered as a systematic uncertainty. Our ROI has been chosen as the best compromise between considering a large enough region to include the extended γ-ray excess of HAWC J1844−034 and minimizing the number of sources we would have to simultaneously fit. During the analysis, the ROI is predetermined before the systematic source search; hence, we have used a larger ROI of 28° < l < 30.3° to see how much of an effect this would have on the source of our main interest, HAWC J1844−034. Although extending the ROI has shown a notable effect on HAWC J1846−025, there is a minimal effect on HAWC J1844−034.

Table 3. A List of Known γ-Ray Sources near HAWC J1844−034

Source	α (deg)	δ (deg)	Δ (deg)	Size (deg)
HAWC J1844−034	${281.02}_{-0.05}^{+0.05}$	$-{3.64}_{-0.04}^{+0.05}$	⋯	${0.48}_{-0.02}^{+0.02}$
LHAASO J1843−0338	280.75	−3.65	0.27	0.3
TASG J1844−038	${281.09}_{-0.10}^{+0.10}$	$-{3.76}_{-0.09}^{+0.09}$	0.14	0.34 ± 0.12
HESS J1843−033	280.95	−3.55	0.11	0.24 ± 0.06
HESS J1844−030	281.17	−3.10	0.56	0.02 ± 0.013
4FGL J1843.7−0326	280.93	−3.44	0.22	⋯
4FGL J1844.4−0345	281.12	−3.75	0.15	⋯

Notes. Here Δ is the angular distance from the center of HAWC J1844−034. The size of LHAASO J1843−0338 is referred from the extension template of 0.3° (Cao et al. 2021). The size of TASG J1844−038 does not include the PSF contribution (Amenomori et al. 2022, see Equation (1)). HESS J1843−033 consists of two Gaussian peaks (H.E.S.S. Collaboration et al. 2018). The positional error of LHAASO J1843−0338 was not reported.

Download table as:  ASCIITypeset image

The nature of the γ-ray emission from HAWC J1844−034 has been modeled and studied to test if this source is a potential hadronic accelerator. In the hadronic scenario, the emission of γ-rays is dominated by the process π⁰ → 2γ. If the neutral pions are produced from the collisions of ambient protons that are accelerated by HAWC J1844−034 (Kafexhiu et al. 2014), then the observed γ-ray spectrum can be used to constrain the proton spectrum of HAWC J1844−034 (Celli et al. 2020). The hadronic scenario is powered by the proton–proton scatterings of accelerated protons and protons of the ambient gas. During the fit, we have assumed a proton density of 1 cm⁻³.

In the case of the leptonic scenario, the TeV γ-ray emission is produced mostly from the inverse Compton scattering between low-energy ambient photons and high-energy electrons accelerated by HAWC J1844−034 (Khangulyan et al. 2014). Due to the Klein–Nishina effects, the γ-ray spectrum originating from the inverse Compton scattering has a spectral cutoff at tens of TeV (Aliu et al. 2014). For the inverse Compton scattering of the leptonic scenario, we used cosmic microwave background (CMB) photons as seed photons that are upscattered by the accelerated electrons.

Figure 6 presents the locations of the sources inside our ROI. The labels at the top of the figure are the VHE γ-ray sources found by H.E.S.S., Tibet AS-Gamma, and LHAASO and previously detected sources from HAWC catalogs (H.E.S.S. Collaboration et al. 2018; Cao et al. 2021; Amenomori et al. 2022). The purple labels indicate the sources specifically found by this work. The lower red labels indicate GeV γ-ray sources observed by Fermi-LAT (Abdollahi et al. 2020, denoted as 4FGL), the green labels are for pulsars from the ATNF pulsar catalog (Manchester et al. 2005), and the blue label marks a known supernova remnant, SNR G28.6−0.1 (Condon et al. 1998; Bamba et al. 2001); Ebisawa et al. 2005. Furthermore, Table 3 contains the positional properties of the sources observed by other experiments labeled in Figure 6, namely, LHAASO, Tibet AS-Gamma, H.E.S.S., and Fermi-LAT.

Zoom In Zoom Out Reset image size

Figure 7 shows the spectra of the GeV to TeV γ-ray sources within the extension of HAWC J1844−034 as shown in Figure 6. HESS J1843−033 has been considered as the counterpart of TASG J1844−038 (Amenomori et al. 2022). The observed spectrum of HESS J1843−033 (brown) has a similar amplitude to that of HAWC J1844−034 (black) and TASG J1844−038 (green), so it has also been assumed as the sub-TeV counterpart in our analysis. The flux points of HESS J1844−030 (orange) and 4FGL J1843.7−0326 (blue) are smoothly connected, so it is likely that these two observations come from the same source. However, it is difficult to find an association of these two sources with HAWC J1844−034 (as well as TASG J1844−038 and LHAASO J1843−0338) because HESS J1844−030 has not been observed with significant emission at energies above 10 TeV (H.E.S.S. Collaboration et al. 2018). The spectrum of 4FGL J1844.4−0345 (red) has been found to follow a log parabolic shape in the GeV range and cuts off sharply; it is not detected above 100 GeV from Fermi-LAT (Abdollahi et al. 2020). Therefore, it is difficult to identify a source in the GeV energies that shows a clear association with HAWC J1844−034 at this time.

Zoom In Zoom Out Reset image size

The two scenarios have been fitted with the Naima framework using only the TeV γ-ray flux data from HESS J1843−033 (H.E.S.S. Collaboration et al. 2018), TASG J1844−038 (Amenomori et al. 2022), LHAASO J1843−0338 (Cao et al. 2021), and HAWC J1844−034 (this work). The particle spectra of each scenario have been assumed to be a power law with an exponential cutoff with a pivot energy of 7 TeV,

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{{\rm{P}}}({E}_{{\rm{P}}}) & = & {\rho }_{0,{\rm{P}}\left(1\mathrm{kpc}\right)}{\left(\displaystyle \frac{{E}_{{\rm{P}}}}{7\,\mathrm{TeV}}\right)}^{-{{\rm{\Gamma }}}_{{\rm{P}}}}\\ & & \times \exp \left(-\displaystyle \frac{{E}_{{\rm{P}}}}{{E}_{c,{\rm{P}}}}\right),\end{array}\end{eqnarray} \tag{ 6 }$$

where ρ_P(E_P) is the spectrum of the accelerated particles, either electrons or protons; E_P is energy of the accelerated particles; and ${\rho }_{0,{\rm{P}}\left(1\mathrm{kpc}\right)}$ is the value of ρ_P(E_P) at a distance of 1 kpc and pivot energy of 7 TeV.

Table 4 shows the best-fit values for the proton spectrum of the hadronic scenario and the electron spectrum of the leptonic scenario. The three free parameters of Equation (6) are listed for each scenario in the table. The obtained spectral index of 1.56 for the hadronic scenario may seem hard considering other hadronic spectra in general. However, hadronic scenarios with such extreme spectral indices are discussed in some cases (Abdalla et al. 2021). The best-fit cutoff energy for the hadronic scenario at 169 TeV is too low for HAWC J1844−034 to be defined as a PeVatron. Note that the particle spectra have been fitted while assuming that HAWC J1844−034 is located at a unit distance of 1 kpc, since the exact counterpart source has not been identified. The source distance is later accounted for with a known source distance in Section 5 when we try to identify a counterpart source.

Table 4. Best-fit Parameters for the Hadronic and Leptonic Scenarios

	Hadronic Scenario	Leptonic Scenario
${\rho }_{0,{\rm{P}}\left(1\mathrm{kpc}\right)}$ [1/eV]	${1.85}_{-0.37}^{+0.34}\times {10}^{34}$	${7.04}_{-1.13}^{+1.17}\times {10}^{31}$
ΓP	${1.56}_{-0.17}^{+0.11}$	${2.17}_{-0.20}^{+0.14}$
Ec,P [TeV]	${169}_{-35}^{+32}$	${72}_{-13}^{+13}$
BIC	85.9	84.4

Note. Here ${\rho }_{0,{\rm{P}}\left(1\mathrm{kpc}\right)}$ is the normalization at a source distance of 1 kpc, Γ_P is the spectral index, and E_c,P is the cutoff energy.

Download table as:  ASCIITypeset image

In Figure 8, we have the data points from LHAASO J1843−0338 (purple), HESS J1843−033 (brown), TASG J1844−038 (green), and HAWC J1844−034 (black). The dashed purple vertical line is the maximum photon energy as observed by LHAASO on this source. The γ-ray spectra obtained from fitting the data points with the Naima framework are presented as blue and red lines for the leptonic and hadronic scenarios, respectively, with their statistical uncertainties shown as lighter bands. The overall trends of the spectra are very similar at energies below 100 TeV. At higher energies, the two spectra start to diverge as the leptonic spectrum cuts off at lower energies. From the Naima fit, the leptonic model is slightly more favorable than the hadronic model with a difference in BIC of 1.5.

Zoom In Zoom Out Reset image size

There are three pulsars that are within the best-fit size of HAWC J1844−034: PSR J1843−0355, PSR J1844−0310, and PSR J1844−0346 (Manchester et al. 2005). Also, PSR J1841−0345 is suggested as a possible counterpart source to LHAASO J1843−0338 in de Oña Wilhelmi et al. (2022). Therefore, we have examined each of these four pulsars as a candidate counterpart source to HAWC J1844−034. The properties of the pulsars used in the analysis have been referenced from the ATNF catalog (Manchester et al. 2005) and listed in Table 5.

Table 5. Various Properties of the Four Pulsars Found Near HAWC J1844−034 Used for Calculations (Manchester et al. 2005; Devin et al. 2021)

			Spin-down		Angular	γ-Ray
Pulsar	R.A.	Decl.	Power	Distance	Distance	Luminosity
	α (deg)	δ (deg)	$\dot{E}$ (erg s−1)	dsrc (kpc)	Δ (deg)	${L}_{\gamma }^{\mathrm{HAWC}}$ (erg s−1)
PSR J1841−0345	280.4112	−3.8118	2.7 × 1035	3.776	0.63	${2.14}_{-0.27}^{+0.44}\times \,{10}^{34}$
PSR J1844−0346	281.13704	−3.7752	4.2 × 1036	4.3	0.18	${2.78}_{-0.35}^{+0.53}\times \,{10}^{34}$
PSR J1843−0355	280.77776	−3.93239	1.77 × 1034	5.798	0.38	${5.05}_{-0.64}^{+0.97}\times \,{10}^{34}$
PSR J1844−0310	281.1895	−3.1770	2.79 × 1033	5.973	0.49	${5.36}_{-0.67}^{+1.03}\times \,{10}^{34}$

Note. Angular distance is the distance from the center of HAWC J1844−034.

Download table as:  ASCIITypeset image

In Table 5, we have computed a number of parameters of the four pulsars for comparison. Here ${L}_{\gamma }^{\mathrm{HAWC}}$ is the γ-ray luminosity estimated from integrating the spectrum of HAWC J1844−034 in Figure 4 as

$$\begin{eqnarray}&&{L}_{\gamma }^{\mathrm{HAWC}}=4\pi {d}_{\mathrm{src}}^{2}{\int }_{\,11\,\mathrm{TeV}}^{202\,\mathrm{TeV}}{{dE}}_{\gamma }\,{E}_{\gamma }{{\rm{\Phi }}}_{\gamma }.\end{eqnarray} \tag{ 7 }$$

The angular distance is the distance to each pulsar from the best-fit location of HAWC J1844−034.

Assuming that the spin-down of a pulsar powers the acceleration of the particles that then produce TeV γ-rays, the γ-ray luminosity is required to be less than the spin-down luminosity of a pulsar to be a possible counterpart source (Cataldo et al. 2020; de Oña Wilhelmi et al. 2022; Albert et al. 2023). Objects PSR J1843−0355 and PSR J1844−0310 both have spin-down powers that are lower than their γ-ray luminosities if they are associated with HAWC J1844−034. Therefore, they are disfavored to be the counterpart source of HAWC J1844−034.

Object PSR J1841−0345 is located 0.63° away from HAWC J1844−034, which is outside of the extension of HAWC J1844−034. Therefore, we rule out PSR J1841−0345 as a counterpart source candidate.

Finally, PSR J1844−0346 has a spin-down power that is approximately two orders of magnitude higher than the γ-ray luminosity. Also, it is located closest to the best-fit positions of HAWC J1844−034, HESS J1843−033, and TASG J1844−038 with angular distances of 0.18°, 0.28°, and 0.05°, respectively. The angular resolutions for HAWC, H.E.S.S., and Tibet AS-Gamma are 0.1°–0.6° for ≥10 TeV, ∼0.08° for ≥0.2 TeV, and ∼0.28° for 100 TeV, respectively (H.E.S.S. Collaboration et al. 2018; Abeysekara et al. 2019; Amenomori et al. 2022). Here the angular resolution of H.E.S.S. is smaller than the angular distance between HESS J1843−033 and PSR J1844−0346. However, HESS J1843−033 is an extended source with a best-fit size of 0.24° (Table 3). Furthermore, although the angular distance of 0.18° is larger than the positional uncertainties given in Table 3, the 0.18° separation is contained within the 0.48° extension of HAWC J1844−034. Therefore, of the four pulsars found near the region of HAWC J1844−034, PSR J1844−0346 is considered to be the most favorable candidate according to this work.

The shell-type SNR G28.6−0.1 (Bamba et al. 2001) is another candidate counterpart source to HAWC J1844−034. It is located at an angular distance of 0.28° from HAWC J1844−034, which is within the best-fit size of HAWC J1844−034. It was suggested to be a potential candidate counterpart source for HESS J1843−033 with VHE γ-rays originating from leptonic processes due to the observed synchrotron emission (Bamba et al. 2001; Devin et al. 2021). Assuming that SNR G28.6−0.1 is responsible for the observed VHE γ-rays by using the distance to SNR G28.6−0.1 of 9.6 ± 0.3 kpc (Ranasinghe & Leahy 2018), the γ-ray luminosity and proton density can be estimated as presented in Table 6.

Table 6. Computed Values for SNR G28.6−0.1

	Hadronic		Leptonic
${L}_{\gamma }^{\mathrm{HAWC}}$ (erg s−1)	nh,9.6 kpc (cm−3)	Wp (erg)	${\rho }_{0,e\left(9.6\mathrm{kpc}\right)}$ (1/eV)	Γe	Ec,e (TeV)	B (μG)	We (erg)
${1.38}_{-0.17}^{+0.27}\times \,{10}^{35}$	46	${2.37}_{-0.73}^{+1.59}\times {10}^{49}$	${8.03}_{-0.52}^{+0.54}\times {10}^{33}$	${2.37}_{-0.01}^{+0.01}$	${92.20}_{-6.70}^{+6.23}$	${2.85}_{-0.16}^{+0.17}$	${2.60}_{-0.45}^{+0.54}\times \,{10}^{50}$

Notes. The best-fit parameters for the proton spectrum are the same as in Table 4. The best-fit parameters for the electron spectrum (${\rho }_{0e\left(9.6\mathrm{kpc}\right)}$, Γ_e , E_c,e ) have been refitted at 9.6 kpc to include the synchrotron emission spectrum.

Download table as:  ASCIITypeset image

We have calculated the energy budgets corresponding to each scenario at the known distance of SNR G28.6−0.1 to compare between leptonic and hadronic scenarios for the SNR. The proton density in the ambient gas of HAWC J1844−034, ${n}_{h,{d}_{\mathrm{src}}}$, is required for the calculation. We use $N\left(\mathrm{HI}\right)\,=1.8\,\times \,{10}^{21}\,{\mathrm{cm}}^{-2}$ as the H i column density and $N\left({{\rm{H}}}_{2}\right)\,=6.7\,\times \,{10}^{21}\,{\mathrm{cm}}^{-2}$ as the H₂ column density, collected from the HI4PI survey and FUGIN data, respectively, at the center of HAWC J1844−034 (Bekhti et al. 2016; Umemoto et al. 2017). The integration range of the radial velocity that corresponds to the location of SNR G28.6−0.1 that has been used to calculate the column densities is found to be 77–88 km s⁻¹ based on Reid et al. (2019). Assuming that the density is spherically uniform within the best-fit size of HAWC J1844−034, the proton density of the ambient gas can be estimated as

$$\begin{eqnarray}\begin{array}{rcl}{n}_{h,{d}_{\mathrm{src}}} & = & \displaystyle \frac{{\rm{\Theta }}}{V}=\displaystyle \frac{3(N\left(\mathrm{HI}\right)+2N\left({{\rm{H}}}_{2}\right))}{4}\\ & & \times \displaystyle \frac{180^\circ }{0\buildrel{\circ}\over{.} 48\,\pi {d}_{\mathrm{src}}},\end{array}\end{eqnarray} \tag{ 8 }$$

where Θ is the total number of protons, V is the spherical volume of the emission region, and 048π/180° is the best-fit size of HAWC J1844−034 in radians. The estimated proton density for SNR G28.6−0.1 is ∼46 cm⁻³.

The energy budget for each scenario is computed by integrating the particle spectrum ρ_P(E_P) provided in Table 4,

$$\begin{eqnarray}&&{W}_{{\rm{P}}}={\int }_{\,{E}_{\min ,{\rm{P}}}}^{{E}_{\max ,{\rm{P}}}}{{dE}}_{{\rm{P}}}\,{E}_{{\rm{P}}}{\rho }_{{\rm{P}}}({E}_{{\rm{P}}}),\end{eqnarray} \tag{ 9 }$$

where ${E}_{\min ,{\rm{P}}}$ and ${E}_{\max ,{\rm{P}}}$ are the minimum and maximum energies of the particle spectrum, respectively.

As for the hadronic scenario, we have used 1.22 GeV for the minimum proton energy ${E}_{\min ,p}$, since the inclusive π⁰ cross section of a proton–proton collision starts to become meaningful at energies above 1.22 GeV (see Kafexhiu et al. 2014, Section 2.B.2). For the maximum proton energy ${E}_{\max ,p}$, we have used 2.6 PeV, which is 10 times the maximum photon energy observed by LHAASO (Cao et al. 2021).

For the modeling of SNR G28.6−0.1, radio and infrared data are available, as well as the γ-ray data. Hence, the Naima framework has computed both the synchrotron and inverse Compton scattering spectrum for the leptonic scenario (Lee et al. 2019; Castelletti et al. 2021). Despite the near-infrared observations at Lee et al. (2019), the near-infrared photons are not considered as seed photons for inverse Compton scattering because photons with energies higher than the CMB are suppressed via the Klein–Nishina effect (Moderski et al. 2005). With the additional data points, we are able to model and fit the synchrotron regime, as well as the inverse Compton scattering regime for the leptonic scenario. Table 6 shows the best-fit values of the newly defined leptonic scenario. Note that because we use radio data to fit the leptonic model for SNR G28.6−0.1, we have extended ${E}_{\min ,e}$ to 10 MeV. Here we assume ${E}_{\max ,e}=740$ TeV, which is calculated in Khangulyan et al. (2014) and de Oña Wilhelmi et al. (2022) based on the maximum photon energy of LHAASO J1843−0338.

The obtained magnetic field strength of ${2.85}_{-0.16}^{+0.17}\,\mu {\rm{G}}$ from the Naima fit is comparable to that of the average Galactic magnetic field strength of 2–3 μG (Longair 2011). The hadronic energy budget has been derived from the proton spectrum in Table 4, since the hadronic model has been fitted with only the γ-ray data.

The photon spectra obtained from the best-fit scenarios of SNR G28.6−0.1 are presented in Figure 9. The γ-ray data points are same as those of Figure 8. The pink and orange data points are from radio and infrared observations, respectively (Lee et al. 2019; Castelletti et al. 2021), used for the modeling of synchrotron emission. The final leptonic and hadronic spectra for SNR G28.6−0.1 as a counterpart to HAWC J1844−034 are displayed as blue and red lines, respectively. The leptonic spectrum (solid) consists of two subspectra, synchrotron emission (dashed–dotted) and inverse Compton scattering (dotted). The synchrotron and inverse Compton scattering spectra have been fitted simultaneously while assuming that a single population of electrons is responsible for the production of both emissions.

Zoom In Zoom Out Reset image size

We have not used any X-ray data in the fit because the reported X-ray flux values for SNR G28.6−0.1 from Bamba et al. (2001) and Ueno et al. (2003) do not contain any flux uncertainties. Since the X-ray band is located in the keV cutoff region of the synchrotron emission, having the X-ray data in the fit model can alter the spectral break parameter and directly influence the spectral break of the inverse Compton scattering as well. Therefore, incorporating reliable X-ray data should help us obtain a more complete understanding of SNR G28.6−0.1 in the future.

The typical energy of a Type Ia supernova explosion is expected to reach ∼10⁵¹ erg (Khokhlov et al. 1993). According to our fitted model, the energy budget for the leptonic scenario is ∼30% of the typical explosion energy. As for the hadronic scenario, the energy budget is estimated to be ∼2% of the explosion energy, which is more than sufficient (Tatischeff 2008). Furthermore, from Amenomori et al. (2022), the diffusion timescale for protons can be estimated using

$$\begin{eqnarray}\begin{array}{rcl} & & {\tau }_{\mathrm{diff}}\sim 1.2\times {10}^{4}\,{\chi }^{-1}{\left(\displaystyle \frac{{R}_{\mathrm{sys}}}{20\,\mathrm{pc}}\right)}^{2}\,\,\,\,\,\,\,\,\,\,\\ & & \times {\left(\displaystyle \frac{{E}_{p}}{\mathrm{GeV}}\right)}^{-0.5}{\left(\displaystyle \frac{B}{10\,\mu {\rm{G}}}\right)}^{0.5}\,\mathrm{yr},\end{array}\end{eqnarray} \tag{ 10 }$$

where χ is the factor of diffusion suppression, and R_sys is the radius of the emission region. For the calculation, a χ of 0.1 and B of 10 μG are assumed (Amenomori et al. 2022). Using R_sys = (048 π d_src)/180° ∼ 80 pc and E_p = 169 TeV (cutoff energy from the best-fit proton spectrum), we can obtain a τ_diff of ∼5 kyr, which is indeed smaller than the known age of SNR G28.6−0.1 of 19 kyr (Ranasinghe & Leahy 2018).