Evidence of 100 TeV $\gamma$-ray emission from HESS J1702-420: A new PeVatron candidatePreprint: H.E.S.S. collaboration,
;
* Corresponding authors

Three-dimensional likelihood analysis, routinely used for high energy (HE; $0.2\lesssim E_{\gamma}\lesssim 100$ GeV) $\gamma$-ray data processing (Mattox et al. 1996; Abdo et al. 2009), has been recently introduced in the VHE $\gamma$-ray astronomy domain (Mohrmann et al. 2019; Donath et al. 2015; Vovk et al. 2018). In its binned version, this technique allows the adjustment of a spectro-morphological model to a data cube, which carries information on the number of reconstructed events within 3D bins. The term 3D refers to the fact that the data are distributed along 2 spatial dimensions (e.g., Galactic longitude and latitude) plus 1 energy dimension. The model can be seen as a collection of spectral and spatial parametric shapes that are assumed to describe all $\gamma$-ray sources in the (source model), plus the residual hadronic background of $\gamma$-like events (background model). The model is convolved with the IRFs to predict the number of photons that would be detected by the telescope array within each spatial and spectral bin, based on the assumed model and its given parameter values. The 3D analysis allows the fine-tuning of all free parameters of the model, in such a way that the cube of model-predicted counts mimics as closely as possible the measured data cube. This approach is known as “forward-folding” (Piron et al. 2001).

We performed a 3D binned likelihood analysis of the region surrounding HESS J1702-420, in order to determine the best spectro-morphological model to describe the observed TeV emission. The next sections discuss the analysis setup (Section 2.1.1) and results (Section 2.1.2), while Section 2.1.3 focuses on the most relevant components of the model, that is those describing HESS J1702-420.

As a complementary study, we performed a 2D analysis of the energy-integrated morphology of HESS J1702-420 in different energy bands. This technique is useful to assess the overall source morphology and verify the persistence of the TeV emission up to the highest energies, even if it does not allow to disentangle HESS J1702-420A from HESS J1702-420B. The level of cosmic ray background in the region was estimated using the adaptive ring background estimation method (Berge et al. 2007; Carrigan et al. 2013). We also verified that consistent flux and significance distributions can be obtained with the FoV background estimation method (Section 2.1). After subtracting the $\gamma$-like hadronic background, we measured $\gamma$-ray flux integrated above 2, 5, 15 and $40\,\text{TeV}$ inside a $1.6^{\,\text{o}}\times\,1.6^{\,\text{o}}$ region encompassing HESS J1702-420. The result is shown in Figure 3. The figure suggests a shrinking of the VHE emission at high energy, with a shift of the $\gamma$-ray peak toward the position of the unidentified source Suzaku src B. Based on the 3D analysis results (Section 2.1), this effect is understood as the transition between a low energy regime — dominated by the steep spectrum of HESS J1702-420B — to a high energy one, in which HESS J1702-420A stands out with its exceptionally hard power law spectrum. Quantitatively, the distance between the low and high energy emission peaks — estimated from the distance between the centroids of HESS J1702-420A and HESS J1702-420B — amounts to $(0.14\pm 0.04_{\text{stat}}\pm 0.02_{\text{sys}})^{\,\text{o}}$.

Launched in 2008, the Fermi-LAT is a pair-conversion instrument sensitive to the HE $\gamma$-ray domain (Atwood et al. 2009). We analyzed $\approx 12\,\text{yr}$ of events with energies in the $10-900\,\text{GeV}$ interval, within a $10^{\,\text{o}}\times\,10^{\,\text{o}}$ RoI encompassing HESS J1702-420. Event selection and binning criteria are detailed in Appendix E. The analysis, performed with fermipy (Wood et al. 2017), made use of Pass 8 IRFs (Atwood et al. 2013).

To build the source model, we selected all sources from the Fourth Fermi General Catalog (4FGL), Thompson 2019, and second Fermi-LAT Catalog of High Energy Sources (2FHL), Ackermann et al. 2016, within $20^{\,\text{o}}$ from the RoI center. In the model, we also included recent diffuse $\gamma$-ray emission templates, both Galactic and extra-Galactic. More details on the source modeling can be found in Appendix E. After the maximum likelihood fit, we produced a TS map to investigate the presence of statistically significant excesses. For each spatial bin, the algorithm compared the maximum log-likelihood obtained by fitting the model, with the addition of a point source ($\Gamma=2$ frozen, amplitude free) at that position, with that of the starting model alone (null hypothesis). We verified that the TS map does not significantly depend on the spectral index or spatial morphology chosen for the test source. The TS map displayed in Figure 4 shows that, within the source region, there is no evidence for a significant excess, but some low-significance fluctuations are present. For comparison, Figure 10 shows a TS map computed before removing the contribution from point sources.

Finally, we included an additional model component defined by a power law spectrum and an elliptical Gaussian morphology identical to the spatial model of HESS J1702-420B (Section 2.1.3). Its $1\sigma$ contour is indicated by the red ellipse in Figure 4. We left free to vary the normalization and index of source spectrum, performed a maximum likelihood fit and compared the resulting model likelihood with the null hypothesis (no source). We found only marginal significance ($4.3\sigma$) for a positive excess corresponding to the chosen Gaussian template. In the absence of a clear detection, we estimated the 99% confidence-level upper limit for the HE emission, associated with the HESS J1702-420B template shape, at the level of

where $E_{\text{ref}}\approx 95\,\text{GeV}$ is the geometric mean of the Fermi-LAT energy range. This value was used to constrain the low-energy extrapolation of the HESS J1702-420B spectrum to the Fermi-LAT energy range (Section 4.2).

Observations of the southern Galactic plane in the 109–115 GHz radio band with Mopra (Braiding et al. 2018), together with the Southern Galactic Plane Survey of the $\lambda=21\,$cm line emission with the ATCA and Parkes telescopes (McClure-Griffiths et al. 2005), allow the study of the molecular and atomic gas distribution in the direction of HESS J1702-420.

The gas densities along the line of sight were measured by Lau et al. 2018, integrating the velocity peaks within a $0.30^{\,\text{o}}\times\,0.15^{\,\text{o}}$ ellipse centered at $l\,$=$\,344.30^{\,\text{o}}$ and $b\,$=$\,-0.18^{\,\text{o}}$. This choice of integration region reflected the approximate shape of the TeV source, from Aharonian et al. 2008. Based on the new H.E.S.S. observations and improved analyses presented in this paper (Section 2), we repeated the ISM analysis, with the same radio dataset and approach as in Lau et al. 2018 but adopting a smaller extraction window to focus on HESS J1702-420A. Our conclusions agree with Lau et al. 2018, in that dense target material, although present at various distances along the line of sight (see Figure 11), does not exhibit any obvious correlation with the VHE $\gamma$-ray maps (see Figure 12). In particular, no hydrogen cloud clearly correlates with HESS J1702-420A or HESS J1702-420B.

In the X-ray domain, deep Suzaku observations of the HESS J1702-420 region revealed the presence of two faint point-like objects (src A and src B, indicated in Figure 3) and the absence of diffuse X-ray emission in the Suzaku FoV, whose dimensions were however insufficient to fully cover the whole TeV source (Fujinaga et al. 2011). Suzaku src B, in particular, is positionally close to the newly discovered component HESS J1702-420A (see Figure 3), which might hint at the first multi wavelength association for HESS J1702-420. For Suzaku src B Fujinaga et al. 2011 estimated a very low flux (in the 2-$10\,\mathrm{keV}$ band) of $(1.9\pm 0.7)\times 10^{-14}\,\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}$, and did not report any evidence of source extension linked with a compact PWN. However, we point out that the Suzaku measurement likely suffered from strong systematics at the position of src B. Indeed, referring to Figure 2 in Fujinaga et al. 2011, it appears that src B was probably not fully contained in the Suzaku FoV, thus leading to an underestimated flux. The higher level of X-ray fluctuations in the corner surrounding src B suggests that the actual level of diffuse emission at that position might be larger than elsewhere in the FoV. Therefore, an association between src B and HESS J1702-420A cannot be ruled out at this stage.

HESS J1702-420A has one of the hardest $\gamma$-ray spectra ever detected in a VHE $\gamma$-ray source. This means that the spectral indices of the underlying particle distributions, responsible for the $\gamma$-ray flux via hadronic or leptonic processes, have to be extremely hard themselves. A pure power law distribution of protons (electrons) with slope $\Gamma_{p}=1.58\pm 0.14_{\text{stat}}$ ($\Gamma_{e}=1.61\pm 0.15_{\text{stat}}$) is well suited to produce the $\gamma$-ray emission of HESS J1702-420A, via hadronic (leptonic) radiative processes. The two spectra, with their $1\sigma$ butterfly envelopes, are shown in Figure 5 (left panel), where the H.E.S.S. spectral points — see Table 9 — are shown for reference purpose only as they were not used for the fit.

Based on the currently available H.E.S.S. data, any attempt of fitting an additional parameter for the cut-off energy of the particle spectra led either to a non-covergence of the fit or to an unphysically high cutoff energy value. We therefore computed lower limits on the particle cut-off energy, following the procedure described in Appendix F. To estimate the lower limits, we modified the model likelihood adding a Gaussian prior on the particle spectral index, to prevent it from floating toward nonphysical regions (i.e., very small or even negative values), due to the trial of low cut-off energies and the reduced lever arm for this spectral modeling. In the case of the hadronic model, we assumed as a prior a Gaussian distribution centered at $\Gamma_{p}=2$ and with $\sigma=0.5$, based on standard diffusive shock acceleration theory (DSA; Bell 1978). We estimated the impact of this prior choice by varying the Gaussian central values to $\Gamma_{p}=1.7$ and $\Gamma_{p}=2.3$. We found that for a prior centered at $\Gamma_{p}=2$ (1.7, 2.3) the 95% confidence-level lower limit on the proton cut-off energy is 0.82 (0.55, 1.16) PeV. The fact that --- independenly of the chosen prior --- the cut-off energy lower limit is found at $E_{p}>0.5\,\text{PeV}$ means that in a hadronic scenario the source likely harbors PeV cosmic rays. In the leptonic case, we tested three different Gaussian priors, all having width $\sigma=0.5$. Based on Sironi & Spitkovsky 2011 and Werner et al. 2015, we chose a prior centered at $\Gamma_{e}=1.5$ to probe shock-driven magnetic reconnection in conditions of moderate wind magnetization, $\Gamma_{e}=2.5$ to account for Fermi-like acceleration at the termination shock in conditions of low upstream magnetization, and finally $\Gamma_{e}=2.0$ as an intermediate scenario. Our results showed that assuming $\Gamma_{e}=2.0$ (1.5, 2.5) the 95% confidence-level lower limit on the electron cut-off energy is 106 (64, 152) TeV.

The energy contents in protons and electrons, necessary to sustain the $\gamma$-ray emission of HESS J1702-420A, were computed integrating the particle spectra above $1\,\text{TeV}$:

Given the best-fit proton and electron distributions found for HESS J1702-420A, with power law indices $\Gamma_{p/e}\approx 1.6$, Eq. 4 would diverge unless the presence of a high energy cut-off is assumed. We therefore adopted the 95% confidence level lower limits on the cut-off energies, thus obtaining lower limits on the integrated particle energetics. We verified that the results are not strongly influenced by the choice of spectral index prior. They are:

In a leptonic scenario, HESS J1702-420A would be powered by an electron popultion with unusually hard spectral index, $\Gamma_{e}\approx 1.6\,$, and the electron energy required to power the $\gamma$-ray emission (see Eq. 6) would be high compared to the typical values for TeV detected PWNe (H.E.S.S. Collaboration et al. 2018). A simple one-zone leptonic model is therefore challenged, also because it would imply the unlikely presence of inverse-Compton emitting electrons with $E_{e}\approx 100\,\mathrm{TeV}$. Indeed, given the $\propto 1/E_{e}$ dependence of the synchrotron loss timescale in the Thomson regime (see Eq. 15), such energetic electrons would cool down extremely fast creating a high energy spectral curvature or break, which is not observed for HESS J1702-420A. To further understand whether a pulsar-PWN association between Suzaku src B (see Section 3.3) and HESS J1702-420A is plausible, we made use of the simple one-zone leptonic model derived in this Section to match the synchrotron emission of HESS J1702-420A with the measured X-ray flux of src B, thus estimating the magnetic field value in the vicinity of the source. This turns out to be unrealistically low: $B\approx 0.3\,\mu G$ (see Figure 13). In other words, if the Suzaku measurement is reliable (see Section 3.3) an association between HESS J1702-420A and Suzaku src B is very unlikely in a simple one-zone leptonic scenario. For all these reasons, a standard PWN model is disfavored, but cannot be definitively ruled out mainly due to the uncertainties on the X-ray measurement. We notice that an alternative interpretation is possible, in which the observed $\gamma$-ray emission is due to electrons that are accelerated by the reconnection electric field at X-points in the current sheets of a pulsar striped wind, where the magnetic field value is expected to be low (Sironi & Spitkovsky 2011; Werner et al. 2015; Guo et al. 2015). In this case, a Doppler boost of the VHE emission due to relativistic plasma motions might be invoked to explain the apparent presence of $100\,\text{TeV}$ inverse-Compton emitting electrons (Cerutti et al. 2020). If true, this would be the first time that a TeV measurement probes the reconnection spectrum immediately downstream of the termination shock of a pulsar wind. However, the lack of a clear multi wavelength detection of the compact object providing the necessary electron population renders this hypothesis unlikely.

In a hadronic scenario, VHE $\gamma$-ray emission is attributed to the interaction of energetic protons with target material within a source or a nearby molecular cloud. In this case, the $100\,\mathrm{TeV}$ $\gamma$-ray emission from HESS J1702-420A, together with its proton cut-off energy lower limit at $0.55-1.16\,\mathrm{PeV}$, would make it a compelling candidate site for the presence of PeV cosmic ray protons. Therefore HESS J1702-420A becomes one of the most solid PeVatron candidates detected in H.E.S.S. data, also based on the modest value of the total energy in protons that is necessary to power its $\gamma$-ray emission (see Eq. 5) and the excellent agreement of a simple proton power law spectrum with the data. However, we notice that a proton spectrum with a slope of $\Gamma_{p}\approx 1.6$ over two energy decades is hard to achieve in the standard DSA framework (Bell 1978). This fact may suggest that HESS J1702-420A, instead of being a proton accelerator, is in fact a gas cloud that, being illuminated by cosmic rays transported from elsewhere, acts as a passive $\gamma$-ray emitter. In that case, the hard measured proton spectrum could result from the energy-dependent particle escape from a nearby proton PeVatron (Gabici et al. 2009). Alternatively, the $\gamma$-ray emission from HESS J1702-420A might be interpreted as the hard high energy end of a concave spectrum arising from nonlinear DSA effects (Kang et al. 2009), or originate from the interaction of SNR shock waves with a young stellar cluster wind (Bykov et al. 2015). The absence of a clear spatial correlation between the ISM and the observed TeV emission (see Section 3.2) prevents a confirmation of the hadronic emission scenario, unless an extremely powerful hidden PeVatron is present. In the latter case, even a modest gas density would suffice to produce the measured $\gamma$-ray emission of HESS J1702-420A, which would explain the observed nonlinearity between the ISM and TeV maps (see Section 3.2).

The baseline proton and electron spectra, used to model the $\gamma$-ray emission of HESS J1702-420B, are broken power laws of the form

where $\tilde{E}$ and $E_{0}$ are the energy of the spectral break and the reference energy, respectively. The introduction of a spectral break was necessary, because a simple power law extrapolation from the VHE to the HE $\gamma$-ray range would have led to unrealistic energy budgets and an overshoot of the Fermi-LAT upper limit (Section 3.1). The first power law index, $\alpha_{1}$, was adjusted manually with respect to the Fermi-LAT upper limit — its value is therefore not to be interpreted as a fit result, but rather as a working assumption.

In the hadronic (leptonic) scenario, the best-fit proton (electron) spectrum corresponds to a broken power-law with slopes $\alpha_{1}=1.6$ (1.4) and $\alpha_{2}=2.66\,\pm 0.11_{\text{stat}}$ ($3.39\,\pm 0.11_{\text{stat}}$), and with break energy of $\tilde{E}=(6.77\,\pm 3.64_{\text{stat}})\,\text{TeV}$ ($\,(4.19\,\pm 1.25_{\text{stat}})\,\text{TeV}\,$). The 95% confidence-level lower limit on the proton (electron) cut-off energy — computed as described in Appendix F — is 550 (140) TeV.

The values of proton and electron energetics, necessary to power the $\gamma$-ray emission of HESS J1702-420B, were computed integrating the broken power law particle spectra above $1\,\text{GeV}$. They are:

In a leptonic scenario, HESS J1702-420A and HESS J1702-420B could be seen as different emission zones belonging to the same PWN complex. However, we deem this interpretation unlikely for several reasons. First of all, a leptonic scenario for HESS J1702-420A is disfavored by the arguments in Section 4.1. Also, the only known nearby pulsar is PSR J1702-4128, that to power the whole TeV source would require an extremely high conversion efficiency of its spin down luminosity into $1-10\,\text{TeV}$ $\gamma$-rays;

where $L_{[1,10]\,\text{TeV}}$ was obtained considering both HESS J1702-420A and HESS J1702-420B, and assuming the same pulsar’s distance from Earth $d=5.2\,\text{kpc}$ (Kramer et al. 2003). The result of Eq. 10 is well above the efficiency of all other PWNe identified by H.E.S.S. in the same energy range (H.E.S.S. Collaboration et al. 2018). Finally, several PWNe detected by H.E.S.S. have an energy-dependent morphology with spectral softening away from the pulsar position (e.g., H.E.S.S. Collaboration et al. 2012; H.E.S.S. Collaboration et al. 2019), which seems not to be the case for HESS J1702-420 (see Appendices D and A). However, we point out that this might be due to insufficient statistics or spatial resolution, and that not all TeV-bright PWNe detected by H.E.S.S. have an energy-dependent morphology (e.g., H. E. S. S. Collaboration 2020). Therefore, leptonic scenarios cannot be definitively ruled out. In particular, as argued in Gallant 2007, the PSR J1702-4128 might power only part of the TeV emission. Indeed, significant VHE $\gamma$-ray emission is detected by H.E.S.S. near the pulsar position — see Figure 3 (upper right panel).

In a hadronic scenario, HESS J1702-420B might be interpreted as a proton accelerator, whose spectral break around $E_{p}\approx 7\,\text{TeV}$ is due to energy-dependent cosmic ray escape from the source. In this case, as argued in Section 4.1, the hard $\gamma$-ray spectrum of HESS J1702-420A could be the signature of delayed emission from the highest energy runaway protons, hitting target material in the ISM. This scenario is challenged however by the absence of clear TeV$-\,n_{\text{H}}$ correlation at the location of HESS J1702-420A (see Section 3.2).

Even if an unequivocal identification of HESS J1702-420 remains elusive, mostly due to the uncertain relationship between HESS J1702-420A and HESS J1702-420B, the new H.E.S.S. observations allow us to constrain the source distance from Earth $d$ and the values of the most relevant environmental parameters in a hadronic or leptonic emission scenario, which are respectively the gas density $n_{\text{H}}$ and magnetic field strength $B$. In this section, we make the assumption that the two components are associated. This means assuming that their distance from Earth is roughly the same, and their TeV emissions are connected.

The constraints we found are shown in Figure 6. The left panel focuses on hadronic scenarios: molecular clouds from Lau et al. 2018 are indicated by red circles, with size proportional to (the logarithm of) the proton energy necessary to power the $\gamma$-ray emission of HESS J1702-420B in each case (see Table 10 for more details). For all clouds, the nearer kinematic distance was assumed. The blue exclusion region in the figure was obtained requiring that the measured proton energetics above $1\,\text{GeV}$ (from Eq. 8) do not exceed $10^{50}\,\text{erg}$, which is the kinetic energy transferred to cosmic rays by a typical SNR. Finally, the gray shaded areas exclude portions of the parameter space in which protons of energy $E_{p}\geq 1\,\text{TeV}$ are cooled down due to $p$-$p$ collisions before having time to diffuse across the whole size of HESS J1702-420B. For this calculation, we assumed a standard ISM magnetic field of $3\,\mu G$, and tested different values for the normalization factor of the diffusion coefficient $\chi$ — defined in Gabici et al. 2007. It is clear that, if the source lies in the diluted ISM where $n_{\text{H}}\lesssim 1\,\text{cm}^{-3}$, it has to be relatively close — $d\lesssim 2\,\text{kpc}$ —, unless its proton energy budget exceeds $10^{50}\,\text{erg}$. If the normalization of the diffusion coefficient is low ($\chi\lesssim 0.001$), only the three nearest molecular clouds would be apt to harbor the source, whose distance would again be $d\lesssim 2\,\text{kpc}$.

In the right panel, which focuses on leptonic scenarios, the gray exclusion areas correspond to portions of the parameter space in which electrons with energy $E_{e}\geq 1\,\text{TeV}$ do not have time to fill the whole component HESS J1702-420B before being cooled down. From the figure, it is clear that if the normalization of the diffusion coefficient is low ($\chi\lesssim 0.01$), then the source has to be relatively close — less than $\approx 3\,\text{kpc}$ away, for realistic values of $B$ field.

Further details on the assumptions that were made to produce the Figure 6 can be found in Appendix G.