Observation of gamma rays up to 320 TeV from the middle-aged TeV pulsar wind nebula HESS J1849-000

Events are taken by opening a circular ON-source window with an angular radius of $0.6^{\circ}$ centered at HESS J1849$-$000, $({\alpha},\,{\delta})=(282{\fdg}24,-0{\fdg}04)$ in the J2000 coordinates (H.E.S.S. collaboration, 2018a). On the other hand, the background is estimated by counting events in 20 OFF windows, each of which has the same size as the ON-source window, opened with the equi-zenith angle method (Amenomori et al., 2003). The events in the ON-source and OFF windows are then binned by energy into five logarithmically equal bins per decade. After the event selection described in Appendix A, the detection significance of gamma rays from HESS J1849$-$000 reaches $4.0\,\sigma$ and $4.4\,\sigma$ levels above $25\,{\rm TeV}$ and $100\,{\rm TeV}$, respectively, in units of Gaussian standard deviation $\sigma$ (Li & Ma, 1983).

Figure 1 shows the significance maps of the HESS J1849$-$000 region above $25\,{\rm TeV}$ (top) and $100\,{\rm TeV}$ (bottom) pixelized by $0.05^{\circ}\times 0.05^{\circ}$ pixels. The map is smoothed with the point-spread function (PSF) of the experiment. In the map above $100\,{\rm TeV}$, the pixel with the maximum significance (let us call it “the brightest pixel”) is found at $({\alpha},\,{\delta})=(282{\fdg}33,0{\fdg}08)$, deviating by $0{\fdg}15$ from HESS J1849$-$000. A toy Monte Carlo simulation shows that the $68\%$ statistical uncertainty in the brightest pixel’s orientation is $0.18^{\circ}$. The pointing systematics of the experiment along the right ascension and the declination are also estimated at $0{\fdg}058$ and $0{\fdg}055$, respectively, from the data analysis of the gamma rays coming from the Crab Nebula; see Appendix C. From the above statistical and systematic uncertainties, the total uncertainty in the center position of the observed gamma-ray emission with the $68\%$ confidence level is estimated at $0{\fdg}22$ following the methodology used in the previous studies (H.E.S.S. collaboration, 2018a; Amenomori et al., 2022). Therefore, the center position of the sub-PeV gamma-ray emission observed in this study is consistent with that of HESS J1849$-$000. On the other hand, HESS J1852$-$000 (H.E.S.S. collaboration, 2018a) deviates by $0.74^{\circ}$ in angular distance from the brightest pixel, and the deviation corresponds to $3.0\,\sigma$ significance taking into account both the uncertainty in the center position of the observed sub-PeV gamma-ray emission ($0{\fdg}22$) and that in the position of HESS J1852$-$000. The result thus disfavors HESS J1852$-$000 as the origin of the sub-PeV gamma rays.

The energy spectrum of gamma rays from HESS J1849$-$000 is measured for the first time in an energy range between $40\,{\rm TeV}<E<320\,{\rm TeV}$ as shown with the red points in Figure 2. The differential energy flux in each energy bin is calculated only if the detection significance of gamma rays exceeds $2\,\sigma$, otherwise, the $99\%$ upper limit on the flux is calculated. Table 1 summarizes the result of the calculation. A simple power-law (PL) function can be fitted to the measured spectrum and the best-fit result is ${\rm d}N/{\rm d}E=(2.86\pm 1.44)\times 10^{-16}(E/40\,{\rm TeV})^{-2.24\pm 0.41}\,{\rm TeV}^{-1}\,{\rm cm}^{-2}\,{\rm s}^{-1}$ with $\chi^{2}/{\rm d.o.f.}=0.47/3$. The experiment’s absolute energy scale uncertainty of $12\%$ (Amenomori et al., 2009) dominates the systematic uncertainty in the flux normalization of $\simeq 27\%$. The fraction of contamination to the number of events above $100\,{\rm TeV}$ from the lower energy range due to the finite energy resolution is estimated at $\simeq 20\%$. Furthermore, there may be a spillover of gamma rays from the nearby gamma-ray source HESS J1852$-$000 into the ON-source window (see also Figure 1). Assuming a two-dimensional Gaussian with a $0.28^{\circ}$ extension for the morphology of HESSJ1852$-$000 (H.E.S.S. collaboration, 2018a), the 95% upper limit on its integral gamma-ray flux above $100\,{\rm TeV}$ is calculated as $2.08\times 10^{-15}\,{\rm cm}^{-2}\,{\rm s}^{-1}$ if the spectrum extends beyond the sub-PeV range. The spillover into the ON-source window is thus estimated at $<3.65\times 10^{-16}\,{\rm cm}^{-2}\,{\rm s}^{-1}$, which is less than $20\%$ of the integral flux of HESS J1849$-$000 above $100\,{\rm TeV}$, $1.93^{+0.85}_{-0.65}\times 10^{-15}\,{\rm cm}^{-2}\,{\rm s}^{-1}$. According to Vernetto & Lipari (2016), the attenuation of gamma rays due to the interactions with the interstellar radiations on the way to Earth is $\simeq 12\%$ at $150\,{\rm TeV}$ if the distance to HESS J1849$-$000 is assumed as $7\,{\rm kpc}$ (H.E.S.S. collaboration, 2018b).

The gamma-ray energy spectrum from the sub-TeV ($E<1\,{\rm TeV}$) to the sub-PeV ranges observed by Fermi-LAT (Acero et al., 2013), H.E.S.S. (H.E.S.S. collaboration, 2018a), LHAASO (Cao et al., 2021a), and this work is modeled with the leptonic scenario using Naima (Zabalza, 2015). In this scenario, gamma rays are generated from inverse Compton scattering (ICS) off the interstellar radiations by high-energy electrons (Khangulyan et al., 2014). The target radiation field is assumed to consist of the cosmic microwave background and near- (temperature of $20\,{\rm K}$ and energy density of $0.75\,{\rm eV}\,{\rm cm}^{-3}$) and far-infrared blackbody components ($3,000\,{\rm K}$ and $1.26\,{\rm eV}\,{\rm cm}^{-3}$) whose energy densities are determined using GALPROP (Porter et al., 2017). The distance to HESS J1849$-$000 is assumed as $7\,{\rm kpc}$, the same as that assumed for PSR J1849$-$0001 in previous studies (Gotthelf et al., 2011; H.E.S.S. collaboration, 2018b). The energy spectrum of parent electrons is assumed to follow a simple PL function

where the normalization constant $A_{\rm e}$ and the spectral index $\alpha_{\rm e}$ are free parameters. Under the above model assumptions, Naima computes the posterior distributions of the parameters and gives the median value and its uncertainty corresponding to the central $68\%$ of each distribution. The best-fit result of the modeling is shown with the black curve in the top panel of Figure 2 and the spectral parameters are estimated at ${\rm log}_{10}A_{\rm e}=31.98^{+0.06}_{-0.07}$ and $\alpha_{\rm e}=2.46^{+0.08}_{-0.07}$. Assuming a PL function with an exponential cutoff (ECPL) for the electron spectrum does not improve the goodness of fit, and the $95\%$ lower limit on the cutoff energy is estimated at $740\,{\rm TeV}$. The obtained limit is extremely high but is not implausible because the acceleration of PeV electrons takes place in the Crab Nebula (Cao et al., 2021b). A broken PL function is not significantly preferred as the electron spectrum either and the current gamma-ray observations can be well explained with ICS by electrons following a simple PL spectrum.

In the leptonic scenario considered above, the gamma rays observed can be interpreted as ICS radiations from high-energy electrons accelerated by the PWN of PSR J1849$-$0001. On the other hand, the total energy of electrons above $100\,{\rm GeV}$ is calculated as $2.8^{+1.0}_{-0.7}\times 10^{47}\,{\rm erg}$, which occupies only $\simeq 2\%$ of the total spin-down energy of PSR J1849$-$0001 ($\simeq 1.3\times 10^{49}\,{\rm erg}$). If the rest of the spin-down energy is assigned to the magnetic field inside the diffuse X-ray nebula around the pulsar ($\simeq 75^{{}^{\prime\prime}}$ in radius or $\simeq 2.5\,{\rm pc}$ at the distance of $7\,{\rm kpc}$ according to Gotthelf et al. (2011)), then the magnetic field strength $B$ in the nebula would be $\sim 400\,\mu{\rm G}$. Such a strong field, however, leads to the synchrotron X-ray nebula far outshining that observed with the energy flux of $\sim 10^{-12}\,{\rm erg\,cm^{-2}\,s^{-1}}$ in $2-10\,{\rm keV}$ (Gotthelf et al., 2011; Kuiper & Hermsen, 2015; Vleeschower Calas et al., 2018), which can be reproduced with $B\sim 2\,\mu{\rm G}$ under a simple one-zone model in the spectral fit assuming the aforementioned PL electron population. Furthermore, the one-zone model would not be adequate to interpret the observations since the extensions of the X-ray nebula and TeV gamma-ray emission ($\simeq 0.09^{\circ}$, H.E.S.S. collaboration (2018b)) are much different, and more detailed theoretical modeling is thus required when one considers the leptonic scenario.

On the other hand, the gamma-ray energy spectrum can also be explained in the context of the hadronic scenario. This scenario supposes that CR protons generate $\pi^{0}$-decay gamma rays through interactions with the ambient interstellar medium (Kafexhiu et al., 2014). To investigate the environment around the source, archive data of the ${}^{12}{\rm CO}\,(J=1-0)$ line emission published by the FUGIN survey (Umemoto et al., 2017) is analyzed. The analysis finds a $\sim 20\,{\rm pc}$ size molecular cloud in the velocity range of $93-100\,{\rm km}\,{\rm s}^{-1}$ corresponding to the distance around $7\,{\rm kpc}$, which resides at the west side of HESS J1849$-$000 as shown by green contours in Figure 1. The cloud has the ${}^{12}{\rm CO}$ line intensity of $\sim 20\,{\rm K}\,{\rm km\,s^{-1}}$ and can provide a gas density of $\gtrsim 10$ protons ${\rm cm}^{-3}$ assuming the canonical CO-to-${\rm H}_{2}$ factor of $2\times 10^{20}\,{\rm cm}^{-2}\,({\rm K\,km\,s^{-1}})^{-1}$ (Bolatto et al., 2013). A wealthy amount of target material for CRs is also implied by the high hydrogen column density of $\simeq 4.5\times 10^{22}\,{\rm cm}^{-2}$ toward the source direction (Gotthelf et al., 2011). The region bright in gamma rays above $100\,{\rm TeV}$ shown in Figure 1 is positionally consistent with the cloud within the uncertainty discussed in Section 3.1.

The modeling of the gamma-ray energy spectrum from the sub-TeV to the sub-PeV ranges is performed under the hadronic scenario assuming an ECPL function for the spectrum of CR protons:

where the normalization constant $A_{\rm p}$, the spectral index $\alpha_{\rm p}$, and the cutoff energy $E_{\rm p,\,cut}$ are free parameters. The density of the interstellar medium is assumed as $10$ protons ${\rm cm}^{-3}$ based on the aforementioned radio data analysis. The result of the modeling is shown in the bottom panel of Figure 2, and the spectral parameters are determined as ${\rm log}_{10}A_{\rm p}=33.93^{+0.09}_{-0.11}$, $\alpha_{\rm p}=2.01^{+0.12}_{-0.21}$, and ${\rm log}_{10}(E_{\rm p,\,cut}/{\rm TeV})=3.73^{+2.98}_{-0.66}$, respectively; it is suggested that the cutoff energy reaches the PeV range. The total energy of the CR protons between $1\,{\rm TeV}<E<10\,{\rm PeV}$ is calculated as $(1.1\pm 0.2)\times 10^{49}\,{\rm erg}$.

The aforementioned spectral index ($\simeq 2$) implies CR protons are accelerated by shock waves of the PWN or the precursor SNR. However, given the high cutoff energy beyond $1\,{\rm PeV}$, the acceleration of CR protons in a PWN-SNR composite system may provide an adequate explanation. Ohira et al. (2018) pointed out that the PWN adiabatically compressed by the reverse shock of the precursor SNR can produce PeV CR protons. The total energy of $\sim 10^{49}\,{\rm erg}$ given to particles by the PWN-SNR system (Gelfand et al., 2009) is also consistent with our result. Moreover, such a composite system can also accelerate and leptons and would produce an X-ray nebula as observed in the adiabatically compressed PWN with the amplified magnetic field (Gelfand et al., 2009; Ohira et al., 2018). On the other hand, there are still no theories that systematically study the distribution of high-energy particles generated in a PWN-SNR system simulating detailed particle acceleration processes and the resultant spectrum of their non-thermal radiations. Neutrinos inevitably generated in the hadronic interaction are not currently detected and the contribution to gamma-ray flux from hadrons is not well constrained (Huang & Li, 2022). Future deep observation of neutrinos by IceCube-Gen2 (Aartsen et al., 2021) and of sub-PeV gamma rays by ALPACA (Kato et al., 2021), CTA (The Cherenkov Telescope Array Consortium, 2017), and SWGO (Hinton & the SWGO Collaboration, 2021) have great importance to find clear evidence for the acceleration of PeV CRs in this energetic source.