Evidence of 100 TeV $\gamma$-ray emission from HESS J1702-420: A new PeVatron candidate

The identification of PeVatrons, hadronic particle accelerators reaching the knee of the cosmic ray spectrum (few $10^{15}$ eV), is crucial to understand the origin of cosmic rays in the Galaxy. We provide an update on the unidentified source HESS J1702-420, a promising PeVatron candidate. We present new observations of HESS J1702-420 made with the High Energy Stereoscopic System (H.E.S.S.), and processed using improved analysis techniques. The analysis configuration was optimized to enhance the collection area at the highest energies. We applied a three-dimensional (3D) likelihood analysis to model the source region and adjust non thermal radiative spectral models to the $\gamma$-ray data. We also analyzed archival data from the Fermi Large Area Telescope (LAT) to constrain the source spectrum at $\gamma$-ray energies>10 GeV. We report the detection of a new source component called HESS J1702-420A, that was separated from the bulk of TeV emission at a $5.4\sigma$ confidence level. The power law $\gamma$-ray spectrum of HESS J1702-420A extends with an index of $\Gamma=1.53\pm0.19_\text{stat}\pm0.20_\text{sys}$ and without curvature up to the energy band 64-113 TeV, in which it was detected by H.E.S.S. at a $4.0\sigma$ confidence level. This brings evidence for the source emission up to $100\,\text{TeV}$, which makes HESS J1702-420A a compelling candidate site for the presence of extremely high energy cosmic rays. Remarkably, in a hadronic scenario, the cut-off energy of the proton distribution powering HESS J1702-420A is found to be higher than 0.5 PeV at a 95% confidence level. HESS J1702-420A becomes therefore one of the most solid PeVatron candidates detected so far in H.E.S.S. data, altough a leptonic origin of its emission could not be ruled out either.


Introduction
The acceleration sites of cosmic rays 1 are a century-old unknown in modern astrophysics (Hess 1912). The current understanding is that the bulk of cosmic rays reaching Earth -mostly energetic protons -originate within our Galaxy, outside of the solar system, at unknown sites where they are accelerated up to the energy of the knee feature in the cosmic ray spectrum (Berezinskii et al. 1990;Gaisser et al. 2016). Since the measured knee energy is around 3 − 4 PeV (Apel et al. 2013), the Galactic accelerators responsible for cosmic rays up to the knee are called PeVatrons. Several source populations have been proposed as potential PeVatron candidates: among them, supernova remnants (SNRs) and young massive stellar clusters stand out as well motivated cases (Bell 2014;Aharonian et al. 2019). However, to date no observation has definitively linked any particular source class to the acceleration of PeV protons. The H.E.S.S. Collaboration has already reported evidence for the acceleration of PeV protons in the central molecular zone around Sgr A * (H.E.S.S. Collaboration et al. 2016;H. E. S. S. Collaboration et al. 2018), at a level that is presently insufficient to sustain the flux of PeV cosmic rays observed at Earth. Recent searches for a high-energy cut-off in the spectrum of the diffuse emission around Sgr A * have led to unclear conclusions, with MAGIC reporting a 2σ hint for a spectral turnover around ≈ 20 TeV and VERITAS measuring a straight power law up to 40 TeV (MAGIC Collaboration et al. 2020;Adams et al. 2021).
Observations of the very high energy (VHE; 0.1 E γ 100 TeV) γ-ray sky with ground-based telescope arrays such as the High Energy Stereoscopic System (H.E.S.S.), providing relatively good angular and energy resolution as well as high sensitivity (H.E.S.S. Collaboration et al. 2018), represent a unique tool to improve our understanding of cosmic ray physics. Charged cosmic particles radiate in this energy band due to interactions with the interstellar medium (ISM) or to the up-scattering of diffuse low-energy radiation fields. The H.E.S.S. Galactic Plane Survey (HGPS) catalog (H.E.S.S. Collaboration et al. 2018) lists 78 VHE γ-ray sources, most of which have been identified as of today -or at least likely associated -with multi wavelength counterparts. However, a handful of sources remain completely unidentified. Having no clear counterpart at other wavelength, they are categorized as dark TeV sources (Aharonian et al. 2008). Experience has shown that several such objects, despite being unidentified at first, were later classified as evolved pulsar wind nebulae (PWNe), based on the discovery of energy-dependent morphologies and compact X-ray counterparts (H.E.S.S. Collaboration et al. 2012Collaboration et al. , 2019. However, a leptonic scenario, in which the VHE γ-ray emission is powered by relativistic electrons up-scattering ambient radiation fields, might not necessarily suit all of the remaining dark sources. HESS J1702-420 is a long-known but poorly understood VHE γ-ray source. It was discovered during the first Galactic plane survey campaign with a significance of 4σ, based on a 5.7 hr observation livetime (Aharonian et al. 2006). In Aharonian et al. (2008), a dedicated analysis revealed a hard power law 2 spectral index of Γ = 2.07 ± 0.08 stat ± 0.20 sys , with no sign of cut-off, and a significantly extended morphology which is well described by a 0.30 o × 0.15 o elliptical Gaussian template. With better reconstruction and data selection algorithms, the HGPS catalog confirmed the spectral hardness of the source, Γ = 2.09 ± 0.07 stat ± 0.20 sys , and estimated a source significance 2 Hearafter, the term power law refers to the functional form dN/dE(E) = Φ ref (E/E ref ) −Γ , where Φ ref is the spectral normalization at the reference energy E ref -usually chosen to correspond with the pivot energy which minimizes correlations between the spectral parameters -, and Γ is the spectral index. of 15σ based on 9.5 hr of observations. It simplified however the source morphology to a 0.2 o symmetric Gaussian, due to the noninclusion of elongated shapes in the semi-automated survey analysis chain.
The physical origin of the γ-ray emission from HESS J1702-420 is unknown. The remnant SNR G344.7-0.1 and the pulsar PSR J1702-4128 are within a ≈ 0.5 o aperture from the centroid of the TeV emission. The former is a 3 kyr old (Giacani et al. 2011) and small-sized -8 arcmin in diameter -SNR, whose centrally-peaked radio shell (Dubner et al. 1993;Whiteoak & Green 1996;Giacani et al. 2011) is also emitting thermal X-rays (Yamaguchi et al. 2012;Combi et al. 2010), with the brightest X-ray and radio features close to each other (Giacani et al. 2011). Recently, the Fermi-LAT association 2FHL J1703.4-4145, with the hard spectral index Γ ≈ 1.2, was discovered on the western edge of the SNR (Eagle et al. 2020). The core-collapse origin of the supernova is debated, due to the absence of a compact remnant. Also controversial is the SNR distance (Eagle et al. 2020), for which a limit based on the high absorbing hydrogen column density is d 8 kpc (Yamaguchi et al. 2012). The cosmic ray diffusion time from the SNR to the VHE peak is compatible with the remnant age (Eagle et al. 2020), which suggests that both 2FHL J1703.4-4145 and HESS J1702-420 may be associated with SNR G344.7-0.1. However, the detection of an extended and bright TeV source at d 8 kpc in the Galactic plane is unlikely given H.E.S.S. sensitivity (H.E.S.S. Collaboration et al. 2018). Moreover, the surrounding ISM does not exhibit any clear morphological association with the VHE γ-ray source (Lau et al. 2018), a fact that challenges a hadronic interpretation of the TeV emission. With a spin-down luminosityĖ = 3.4 × 10 35 erg s −1 , the large-offset pulsar PSR J1702-4128 would need a conversion efficiency > 10% in order to power the whole TeV source, higher than all other PWNe identified by H.E.S.S. (Gallant 2007;H.E.S.S. Collaboration et al. 2018). This fact, together with the unconclusive searches for an asymmetric X-ray PWN around the pulsar (Chang et al. 2008), tend to disfavor an association of PSR J1702-4128 with HESS J1702-420. Finally, deep X-ray observations of the VHE source with Suzaku revealed the presence of two faint point-like sources close to the line of sight of HESS J1702-420, and the absence of extended emission, with an X-ray flux at least 12 times lower than the TeV flux in the Suzaku field of view (FoV) (Fujinaga et al. 2011).
This paper reports on new H.E.S.S. observations of HESS J1702-420 that have been processed with improved techniques. Additionally, archival Fermi-LAT data were analyzed, to perform a broadband modeling of the TeV source. The paper is structured as follows: first, the H.E.S.S. data analysis and results are presented (Section 2), with focus on the three-dimensional (3D) likelihood anaysis (Section 2.1) and a morphological study made with a more classical background estimation technique (Section 2.2). Then, Section 3 reports on the analysis of archival multi wavelength observations of the region, while Section 4 describes the adjustment of physically-motivated non thermal radiative models to the data and discusses possible interpretations of the new H.E.S.S. results in a broadband context. Finally, Section 5 summarizes the main conclusions of the paper. eter telescopes -CT1-4, whose cameras were upgraded in 2017 (Ashton et al. 2020) -at the corners of a 120 m × 120 m square. In 2012, a second phase began (H.E.S.S. II) with the addition of a 28 m diameter telescope (CT5) at the center of the array.
All results presented in this paper make use of data collected from 2004 to 2019, using CT1-4 observations that have been carried out in multiple contexts: dedicated pointings on HESS J1702-420, observations of other nearby objects -mainly RX J1713-3946 and PSR B1706-44 -and Galactic plane scan observations for the HGPS campaign. This led to the accumulation of an acceptance-corrected livetime of 44.9 hr, obtained after selecting all runs with pointing direction within 3 o from the HGPS position of HESS J1702-420 (l = 344.30 o , b = −0.18 o ) and averaging over a 0.5 o circle centered at the same location.
Observations were processed with the H.E.S.S. analysis package (HAP), applying a Hillas-type shower reconstruction (Hillas 1985) and the multi-variate analysis (MVA) technique (Becherini et al. 2012) for efficient γ/hadron discrimination. Preselection and MVA discrimination cuts were optimized to improve the collection area at high energies (E > 1 TeV), assuming a γ-ray spectral index of 2.3 . The reduced data and instrument response functions 3 (IRFs) were exported to FITS files complying with the standard format developed by Deil et al. (2017a). Finally, all high-level analysis results -that is sky maps and spectra -were obtained using gammapy (version 0.17), an open source python library for γ-ray astronomy (Deil et al. 2017b;Nigro et al. 2019;Deil et al. 2020). The analysis cross-check was also performed with gammapy, by applying the same high-level analysis pipeline to data that were reduced according to an alternative low-level chain of calibration, reconstruction and γ-hadron separation methods (Ohm et al. 2009). Independent crosschecks for the classical morphological and spectral analyses presented in Section 2.2 and Appendix A were also carried out using the standard HAP software.

Three-dimensional likelihood analysis
Three-dimensional likelihood analysis, routinely used for high energy (HE; 0.2 E γ 100 GeV) γ-ray data processing (Mattox et al. 1996;Abdo et al. 2009), has been recently introduced in the VHE γ-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 γ-ray sources in the (source model), plus the residual hadronic background of γ-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.

Background model and analysis setup
H.E.S.S. data-taking consists of consecutive observations (also called runs), usually of 28 minutes duration. The background model was produced from a large set of empty-field -that is devoid of known γ-ray sources -observations, following an approach similar to the one described in Mohrmann et al. (2019). We first produced a general model in the form of a lookup table, describing the residual hadronic background as a function of few observational parameters -namely, the zenith angle and optical efficiency. We then assigned each observation a specific background model, called its FoV background model, based on a multi-variable interpolation of the general model. In addition, the FoV background model was renormalized run-by-run, to account for possible differences in the level of night sky background (NSB) and atmospheric absorption with respect to the observations that were used to build the general background model.
As a result of the region's observation history (see the introduction of Section 2), standard run selection criteria led to a heterogeneous set of observations in terms of array response, zenith angles and source offsets from the pointing direction. We therefore separated observations obtained before and after the 2017 camera upgrade, for which different IRFs have to be used, and grouped together observations with similar zenith and offset values. More details on the four groups of observations that were defined are reported in Table 1. We stacked observations within each group, thus obtaining four independent datasets. This choice represents a good compromise between the timesaving, due to a decrease in the number of degrees of freedom, and information loss, due to the IRFs averaging, that are connected with the data stacking procedure.
For each of the four observation groups, the list of reconstructed events was reduced to a binned data cube, with spatial dimensions corresponding to an analysis region of interest (RoI) of 4 o × 4 o centered at the HGPS position of HESS J1702-420. This choice represents a compromise between a sufficiently large RoI, to get enough off-source regions for the background estimation, and a sufficiently small RoI, to minimize the number of unrelated sources needing to be modeled. A spatial pixel size of 0.02 o × 0.02 o was adopted, to ensure sufficient per-pixel statistics while still providing good spatial resolution. The third axis of the cube, encoding the reconstructed energy of incident photons, was divided into 20 equally-spaced -in logarithmic scale -bins between 0.5 and 150 TeV. In order to reject poorly reconstructed data, for each observation all events with offset ≥ 2 o from the pointing direction were excluded, together with those whose reconstructed energy was below the safe threshold where E bkg represents the energy at which the maximum rate of hadronic background events occurs and E A eff is the energy at which the effective area at the center of the FoV drops to 10% of its maximal value (e.g., Mohrmann et al. (2019)).
With this setup, we performed a joint-likelihood analysis of the four independent datasets, each one having its own FoV 2017 -2019 * Obtained by averaging over a 0.5 o -radius circle centered at the center of the HGPS position of HESS J1702-420. We observe that runs belonging to the second and third groups were included to better contrain the background level in the source region. The fourth group contains all (and only) runs taken after the 2017 camera upgrade (Ashton et al. 2020). Table 1: Details on the four groups of H.E.S.S. observations, with similar pointing zenith angle and source offsets from the pointing direction, that were used for the 3D likelihood analysis.
background model but all sharing the same source model. This means that we summed the four dataset-specific log-likelihood values, and maximized the total resulting likelihood with respect to the model parameters. For the fit, the Cash statistic (Cash 1979) for Poisson-distributed data with perfectly known background model was used. The 3D analysis was performed in the energy range 2 − 150 TeV. All events with reconstructed energy below 2 TeV were excluded from the likelihood computation, to avoid threshold effects arising from the high energy optimization (above 1 TeV) of the analysis configuration and to ensure that the power law assumption made for the background model spectrum was valid -which is true only well above the background peak. A 0.25 o band around the borders of the RoI was excluded from the analysis, in order to limit possible contamination due to nonmodeled sources outside the RoI. Additionally, a 0.3 o -radius circular region centered at l = 343.35 o and b = −0.93 o , containing a ≈ 3σ significance hotspot, was excluded from the likelihood calculation instead of being modeled. This choice was motivated by the high offset between the hotspot and HESS J1702-420, and the necessity of limiting the number of nuisance parameters of the source model.

Source model derivation and results
The optimal source model for the RoI was determined using a statistical approach based on the improvement of a first-guess model with the iterative addition of new components. As a starting point, we defined a source model including all known VHE sources within the 4 o × 4 o RoI, with the exception of HESS J1702-420. Each iteration then consisted either in the addition of a new source component -described by a symmetric Gaussian morphology and power law spectrum -, or the test of a different assumption on the spatial or spectral shape of an already existing component. Specifically, we looked for the presence of high energy spectral cut-offs or elongated shapes for all components.
Step-by-step, the improvement of the source model was assessed looking at two indicators. Firstly, we used the likelihoodratio test, which allows to estimate the relative significance of nested hypotheses taking into account the number of additional degrees of freedom added at each step. For example, the presence of an additional model component µ with 1 (5) free parameter(s) was considered significant at 5σ confidence level only if TS ≥ 25 (TS ≥ 37.1), where TS is the test statistic Here, L max 0 (L max µ ) represents the maximum likelihood of the model under the null (alternative) hypothesis -that is the absence (presence) of µ. More details can be found in Appendix B. Secondly, we assessed -by visual inspection -the flattening of spatial and spectral residuals toward zero, as a result of the addition of new model components.
The final results of the procedure are shown in Figure 1. The upper left (right) panel of the figure shows the measured (model-predicted) counts map, obtained after stacking the four individual datasets and integrating over the energy axis above 2 TeV. Diagonal line hatches represent portions of the RoI that were excluded from the likelihood computation (Section 2.1.1). The measured counts map is well matched by the prediction, since the spatial distribution of the significance of model residuals (lower left panel) does not contain significant structures. The histogram of significance values (lower right panel) closely matches a standard normal distribution, as expected if residuals are only due to statistical Poisson fluctuations. Additionally, the spatial distribution of model residuals in three independent energy bands is shown in Figure H The upper right panel of Figure 1 also shows the 1σ contours of all components found in the final source model. There are two overlapping objects, called HESS J1702-420A and HESS J1702-420B, that together describe the emission from HESS J1702-420. Being the most relevant model components for the scope of this paper, their details are discussed in the dedicated Section 2.1.3. Two other model components represent the nearby sources HESS J1708-410 and HESS J1708-443. Due to their large angular distance from the center of the RoI, the details of their modeling do not have a strong impact on HESS J1702-420. The fitted model for HESS J1708-410 was found to be consistent with that reported in the HGPS catalog, while the model for HESS J1708-443, being only partially contained in the RoI, was directly fixed to the catalog one. Finally, we found a large-scale component, indicated by the dashed circle in Figure 1 (upper right panel), whose presence was not confirmed by the crosscheck analysis. More details can be found in Appendix C. During the analysis, all parameters describing HESS J1702-420A and HESS J1702-420B, together with the nuisance parameters of all other components and the background model, were left free to vary. Details on the final parameters for the most relevant model components are provided in Tables 2, 3 and 4.

HESS J1702-420A and HESS J1702-420B
The most relevant result for the identification of Galactic Pevatrons is the discovery -with a TS-based confidence level corresponding to 5.4 σ -of a new source component, HESS J1702-420A, hidden under the bulk emission formerly associated with HESS J1702-420. This object has a spectral index of Γ = 1.53 ± 0.19 stat ± 0.20 sys and a γ-ray spectrum that, extending with no sign of curvature up to at least 64 TeV (possibly 100 TeV), makes it a compelling candidate site for the presence of extremely high energy cosmic rays. With a flux above 2 TeV

Galactic Longitude
Galactic Latitude  of (2.08 ± 0.49 stat ± 0.62 sys ) × 10 −13 cm −2 s −1 and a 1σ radius of (0.06 ± 0.02 stat ± 0.03 sys ) o , HESS J1702-420A is outshone below ≈ 40 TeV by the companion HESS J1702-420B. The test of a point-source hypothesis for HESS J1702-420A resulted in a non-convergence of the fit. HESS J1702-420B has a steep spectral index of Γ = 2.62 ± 0.10 stat ± 0.20 sys , elongated shape and a flux above 2 TeV of (1.57 ± 0.12 stat ± 0.47 sys ) 10 −12 cm −2 s −1 that accounts for most of the low-energy HESS J1702-420 emission. By comparing results obtained with the main and crosscheck analysis configurations, we verified that all discrepancies were consistent with the expected level of H.E.S.S. systematic uncertainties (H.E.S.S. Collaboration et al. 2018).
For neither of the two sources did an exponential cut-off function statistically improve the fit with respect to a simple power law (cut-off significance 1σ). The γ-ray spectra of both components are shown in Figure 2, together with spectral points computed under a power law assumption and reoptimizing all the nuisance parameters of the model -see Table H.1 for details. We adapted the binning of the spectral energy distributions to obtain approximately equal counts in each bin. HESS J1702-420B is the brightest component up until roughly 40 TeV, where HESS J1702-420A eventually starts dominating with its Γ ≈ 1.5 power law spectrum up to 100 TeV. The second to last spectral point of HESS J1702-420A (HESS J1702-420B), covering the reconstructed energy range 64 − 113 TeV (36 − 113 TeV), is significant at 4.0σ (3.2σ) confidence level. and HESS J1702-420B (blue solid line), as a function of the incident photon energy E γ . The butterfly envelopes indicate the 1σ statistical uncertainty on the spectral shape. They have been obtained from a 3D fit of the H.E.S.S. data with gammapy (more details in the main text). The spectral points, shown for reference purpose only, have been obtained by rescaling the amplitude of the reference spectral model within each energy bin, reoptimizing at the same time all free nuisance parameters of the model. In the energy bins with less than 3σ excess significance, the 3σ confidence level upper limits are shown.
actually represent two separate sources -superimposed on the same line of sight -or rather different emission zones of a single complex object. Moreover, any morphology assumption based on exact geometric shapes -in this case, two overlapping Gaussian components -represents an idealization, that might differ from the real underlying astrophysical model. In particular, a model assumption based on the energy-dependent morphology of a single source might also be well suited to describe the emission of HESS J1702-420. To address this point, we performed dedicated studies that ultimately provided a confirmation of the 3D analysis results, in that they brought no evidence of energy-dependent variations of the 3D model or spectral softening as a function of the distance from HESS J1702-420A (see Appendix D and A). Therefore, a model describing HESS J1702-420 with a single energydependent component is disfavored, even if it cannot be definitively ruled out.

Flux maps and source morphology
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 γ-like hadronic background, we measured γ-ray flux integrated above 2, 5, 15 and 40 TeV inside a 1.6 o × 1.6 o re-gion 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 γ-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 ± 0.04 stat ± 0.02 sys ) o .

Multi wavelength observations
Even in the absence of a multi wavelength detection of HESS J1702-420, low-energy observations can help to constrain the TeV emission scenarios. Section 3.1 summarizes a dedicated analysis of archival Fermi-LAT data in the HESS J1702-420 region, while Section 3.2 reports our considerations on the surrounding ISM and Section 3.3 discusses archival Suzaku measurements in the context of the new H.E.S.S. results.

Fermi-LAT data analysis and results
Launched in 2008, the Fermi-LAT is a pair-conversion instrument sensitive to the HE γ-ray domain (Atwood et al. 2009). We analyzed ≈ 12 yr of events with energies in the 10 − 900 GeV interval, within a 10 o × 10 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 o from the RoI center. In the model, we also included recent diffuse γ-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 (Γ = 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 H.2 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σ 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σ) 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 [deg] HESS J1702-420A 344.15 ± 0.02 stat ± 0.01 sys −0.15 ± 0.02 stat ± 0.01 sys 0.06 ± 0.02 stat ± 0.03 sys -HESS J1702-420B 344.29 ± 0.03 stat ± 0.01 sys −0.15 ± 0.02 stat ± 0.01 sys 0.32 ± 0.02 stat ± 0.03 sys 0.20 ± 0.02 stat ± 0.03 sys 67.0 ± 5.4 stat ± 9.7 sys * Measured counterclockwise starting from the l = 0, b > 0 axis.

The interstellar medium
Observations 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 H.3), does not exhibit any obvious correlation with the VHE γ-ray maps (see Figure H.4). In particular, no hydrogen cloud clearly correlates with HESS J1702-420A or HESS J1702-420B.

Comparison with X-ray observations of HESS J1702-420
In the X-ray domain, deep Suzaku observations of the HESS J1702-420 region revealed the presence of two faint pointlike 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 keV band) of (1.9±0.7)×10 −14 erg s −1 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.

Discussion
To model the γ-ray emission of HESS J1702-420A and HESS J1702-420B, we replaced the power law spectral functions that were used in the 3D analysis (Section 2.1) with simple physically-motivated non thermal radiative models from naima (Zabalza 2015). We derived the present-age spectral shape of the parent cosmic ray population, exploring both hadronic and leptonic one-zone emission scenarios. Owing to the NaimaSpectralModel class implemented in gammapy, we could forward-fold the naima radiative models directly on the H.E.S.S. 3D data. This represents a significant improvement with respect to a simple fit to precomputed flux points, which is inevitably biased by the spectral assumption made for the flux point computation.
Because of the unclear level of association between HESS J1702-420A and HESS J1702-420B, we modeled them independently. For the hadronic emission models, based on the analytic parametrization of p-p interaction and π 0 decay developed in Kafexhiu et al. (2014), we assumed a fixed target density n H = 100 cm −3 . During the fit, a fiducial distance from Earth of d = 3.5 kpc was assumed. We note that the gas density, as well as the source distance from Earth, do not influence the spectral Article number, page 7 of 24 A&A proofs: manuscript no. paper_giunti_HESSJ1702  (Li & Ma 1983). The cyan markers indicate the position (surrounded by uncertainty ellipses) of the Fermi-LAT sources 4FGL J1702.9-4131 and 2FHL J1703.4-4145. The former is associated with the PSR J1702-4128 (in yellow), the latter with the SNR G344.7-0.1 (in gray). The circle around the SNR represents its angular extension (Giacani et al. 2011). The white pentagon and upward-pointing triangle represent unidentified X-ray Suzaku sources. The orange markers show the positions of nearby X-ray binaries. Finally, the center and 1σ extension of HESS J1702-420A (HESS J1702-420B) are indicated in green (blue). In the bottom-left corner of each panel the 68% containment radius of the H.E.S.S. PSF is shown, which -for the chosen analysis configuration -does not have a strong dependency on the energy. shape and are both degenerate with the source intrinsic luminosity, that may be rescaled a posteriori assuming different values of n H and d (see for example Eq. 5). In the leptonic scenario, based on the analytic approximation presented in Khangulyan et al. (2014), the VHE γ-ray emission was attributed to inverse-Compton up-scattering by electrons of the cosmic microwave background (CMB) and infrared (IR) low-energy photon fields. The uniform CMB field was described as a black-body radiation with energy density of CMB = 0.261 eV cm −3 and temperature of T CMB = 2.73 K. The starlight emission in the near IR ( NIR = 1 eV cm −3 and T NIR = 3000 K) and dust re-emission in the far IR ( FIR = 0.5 eV cm −3 and T FIR = 30 K) were obtained using the 3D interstellar radiation field (ISRF) model from Porter et al. (2018), at the coordinates of HESS J1702-420 and the assumed 3.5 kpc distance. We verified that the level of fluctuations of the ISRF along the line of sight did not significantly impact the modeling conclusions. The results are discussed in Sections 4.1 and 4.2.

HESS J1702-420A
HESS J1702-420A has one of the hardest γ-ray spectra ever detected in a VHE γ-ray source. This means that the spectral indices of the underlying particle distributions, responsible for the γ-ray flux via hadronic or leptonic processes, have to be extremely hard themselves. A pure power law distribution of protons (electrons) with slope Γ p = 1.58 ± 0.14 stat (Γ e = 1.61 ± 0.15 stat ) is well suited to produce the γ-ray emission of HESS J1702-420A, via hadronic (leptonic) radiative processes. The two spectra, with their 1σ butterfly envelopes, are shown in Figure 5 (left panel), where the H.E.S.S. spectral points -see Table H.1 -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 Γ p = 2 and with σ = 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 Γ p = 1.7 and Γ p = 2.3. We found that for a prior centered at Γ 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 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 σ = 0.5. Based on Sironi &Spitkovsky (2011) andWerner et al. (2015), we chose a prior centered at Γ e = 1.5 to probe shock-driven magnetic reconnection in conditions of moderate wind magnetization, Γ e = 2.5 to account for Fermi-like acceleration at the termination shock in conditions of low upstream magnetization, and finally Γ e = 2.0 as an intermediate scenario. Our results showed that assuming Γ 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 γ-ray emission of HESS J1702-420A, were computed integrating the particle spectra above 1 TeV: Given the best-fit proton and electron distributions found for HESS J1702-420A, with power law indices Γ p/e ≈ 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, Γ e ≈ 1.6 , and the electron energy required to power the γ-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 ≈ 100 TeV. Indeed, given the ∝ 1/E e dependence of the synchrotron loss timescale in the Thomson regime (see Eq. G.4), 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 ≈ 0.3 µG (see Figure H.5). 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 γ-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 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 γ-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 TeV γ-ray emission from HESS J1702-420A, together with its proton cut-off energy lower limit at 0.55 − 1.16 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 γ-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 Γ p ≈ 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 γ-ray emitter. In that case, the hard measured proton spectrum could result from the energy-dependent particle escape from a nearby proton Pe-Vatron (Gabici et al. 2009). Alternatively, the γ-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 γ-ray emission of HESS J1702-420A, which would explain the observed nonlinearity between the ISM and TeV maps (see Section 3.2).

HESS J1702-420B
The baseline proton and electron spectra, used to model the γ-ray emission of HESS J1702-420B, are broken power laws of the form whereẼ 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 γ-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, α 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.
The values of proton and electron energetics, necessary to power the γ-ray emission of HESS J1702-420B, were computed integrating the broken power law particle spectra above 1 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 TeV γ-rays; where L [1,10] TeV was obtained considering both HESS J1702-420A and HESS J1702-420B, and assuming the same pulsar's distance from Earth d = 5.2 kpc (Kramer et al. 2003). The result of Eq. 10 is  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 TeVbright 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 γ-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 ≈ 7 TeV is due to energy-dependent cosmic ray escape from the source. In this case, as argued in Section 4.1, the hard γ-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 H correlation at the location of HESS J1702-420A (see Section 3.2).

Distance from Earth and environmental parameters
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 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 γ-ray emission of HESS J1702-420B in each case (see Table H.2 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 GeV (from Eq. 8) do not exceed 10 50 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 ≥ 1 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 µG, and tested different values for the normalization factor of the diffusion coefficient χ -defined in Gabici et al. (2007). It is clear that, if the source lies in the diluted ISM where n H 1 cm −3 , it has to be relatively closed 2 kpc -, unless its proton energy budget exceeds 10 50 erg. If the normalization of the diffusion coefficient is low (χ 0.001), only the three nearest molecular clouds would be apt to harbor the source, whose distance would again be d 2 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 ≥ 1 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 (χ 0.01), then the source has to be relatively close -less than ≈ 3 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.

Conclusions
We present new H.E.S.S. observations of the unidentified source HESS J1702-420, processed using improved techniques, that bring new evidence for the presence of γ-rays up to 100 TeV. The low-level analysis configuration, used to reduce the raw telescope data to lists of γ-like events, was adapted to maximize the telescope's sensitivity at the highest energies. We performed a 3D likelihood analysis -a relatively new high-level technique in the VHE γ-ray domain -with gammapy, to determine the simplest and best suited spatial and spectral models to describe the source and its surroundings. This allowed us to separate for the first time two components -both detected at > 5σ confidence level -inside HESS J1702-420 based on their different morphologies and γ-ray spectra, both of which extend with no sign of curvature up to several tens of TeV (possibly 100 TeV). We report the 4.0σ confidence level detection of γ-ray emission from the hardest component, called HESS J1702-420A, in the energy band 64 − 113 TeV, which is an unprecedented achievement for the H.E.S.S. experiment and brings evidence for the source emission up to 100 TeV. With a spectral index of Γ = 1.53 ± 0.19 stat ± 0.20 sys , this object is a compelling candidate site for the presence of PeV cosmic rays.
We adjusted physically-motivated non thermal radiative models to the H.E.S.S. data, testing simple one-zone hadronic and leptonic models, and determined that the available observations do not allow us to rule out either of the two scenarios. The 95% confidence level energy cut-off of the baseline proton (electron) distribution of HESS J1702-420A was found in the range 0.55 − 1.16 PeV (64 − 152 TeV), depending on the assumption made on the particle spectral index. Remarkably, in a hadronic emission scenario the particle spectral cut-off is at E p > 0.5 PeV, for a range of tested priors. For such a scenario, this implies that the source harbors PeV protons, thus becoming one of the most solid PeVatron candidates detected in H.E.S.S. data. Nevertheless, a leptonic emission scenario for HESS J1702-420A could not be definitively ruled out. We additionally measured the particle energetics that are necessary to power the observed γ-ray emission. We finally discussed possible constraints on the source distance, ambient magnetic field and surrounding gas density.
In the future, the improved angular resolution of the Cherenkov Telescope Array (CTA) and higher energy coverage of the Southern Wide-field Gamma-ray Observatory (SWGO) will possibly close the debate on the nature of HESS J1702-420. In particular, deep measurements in the 100-200 TeV γ-ray band will constrain the spectral shape near the cut-off region, thus probing the hadronic or leptonic origin of the emission and determining whether either of the two detected components operates as a real cosmic ray PeVatron. Observations in the X-ray band, on the other hand, will be important to search for a multi wavelength counterpart of the TeV source, and clarify the relationship between HESS J1702-420A and the unidentified Suzaku src B. γ-ray signal significance above 2 TeV, with contours corresponding to 2σ, 3σ, 5σ, 9σ and 12σ significance levels Li & Ma (1983). The map has been obtained with the Adaptive Ring Background estimation method, and centered at the position of HESS J1702-420A.
Overlaid on the map are the concentric regions -one circle and three annuli -that were used to extract the source spectrum.
Lower panel: Results of the spatially-resolved spectral analysis, showing the spectral index and flux as a function of the distance from HESS J1702-420A.

Appendix A: Spatially-resolved spectral analysis of H.E.S.S. data
With the benefit of an unprecedented level of statistics in this region, we performed a spatially-resolved spectral analysis for HESS J1702-420. A 0.15 o -radius circle and three 0.2 o -radius annuli were used to measure the VHE γ-ray spectrum of the source. We did not adopt narrower extraction regions, in or-der to limit the level of PSF-induced correlation between them. Figure A .1 (upper panel) shows the four nonoverlapping regions, overlaid on a map of the γ-ray flux significance above 2 TeV. All regions are concentric around Galactic coordinates l = 344.15 o and b = −0.15 o , corresponding to the position of HESS J1702-420A. The level of cosmic ray background within each region was computed with the reflected region background estimation technique (Berge et al. 2007), while a forward-folding approach (Piron et al. 2001) was adopted to determine the maximum-likelihood estimates of the spectral slope and flux, in each region, under a power-law assumption.
The detailed results of the spectral analysis are reported in Table A.1, while the spectral variations as a function of the distance from HESS J1702-420A are shown in Figure A.1 (bottom panel). The error bars in the figure represent the statistical errors on the fitted parameters. The level of systematic uncertainties, reported in Table A.1, have been estimated following H.E.S.S. Collaboration et al. (2018). The figure shows that, in this datasets, there is no evidence for significant spectral variations around HESS J1702-420A. This measurement tends to support a two-component approach, with respect to a model based on a single source with energy-dependent morphology. Indeed, in the latter case significant spatially-resolved spectral variations would be expected, as seen for other well known H.E.S.S. sources (e.g., H.E.S.S. Collaboration et al. (2019)).
where L max 0 represents the maximum model likelihood under the null power law -or equivalently E c → ∞ -hypothesis; iv we finally computed the 90%,95% and 99% confidence level lower limits on the particle cut-off energy by finding the values where the TS profile increased from the minimum (which in our case was at infinity) by an amount TS (90%) = 2.706, TS (95%) = 3.841 and TS (99%) = 6.635 respectively.
This procedure, based on Rolke et al. (2005), is partly implemented in the Fit.stat_profile() routine of gammapy.

Appendix G: Cosmic ray diffusion model and energy loss calculation
The source physical size R source depends directly on the distance from Earth d and the measured angular size of the source θ source , as R source = d × tan(θ source ). Here, we assumed θ source = 1.28 o , which corresponds to the major 2σ diameter of HESS J1702-420B. For both hadronic and leptonic cosmic rays, we adopted the energy-dependent diffusion coefficient defined in Gabici et al. (2007), testing different values for the normalization χ. The relation between the diffusion coefficient and the diffusion timescale is given by: Energy losses for protons due to p-p collisions were estimated -neglecting ionization losses that are irrelevant for relativistic protons -as τ pp ≈ 6 × 10 5 n H 100 cm −3 −1 yr , (G.2) as in Gabici et al. (2007). Finally, the electron energy loss timescale was computed using is the synchrotron loss timescale and τ IC ≈ 3 × 10 7 E 10 GeV is the inverse-Compton loss timescale for a given photon field (Ginzburg & Syrovatskii 1964).

Appendix H: Additional material
In this section, additional figures and tables are provided.
Article number, page 17 of 24       : Column density maps of molecular hydrogen in the direction of HESS J1702-420, obtained by integrating the brightness temperature profile of 12 CO(J = 1 → 0) data from the Mopra radio survey within the velocity intervals indicated above each panel (corresponding to the peaks in Figure H.3). The brightness temperature values were converted to H 2 column density assuming the conversion factor X CO = 1.5 × 10 20 cm −2 (K km s −1 ) −1 (Strong et al. 2004). The green (orange) contours indicate the 5 and 12σ (3 and 5σ) significance levels of the TeV γ-ray flux above 2 TeV (40 TeV). The dashed ellipse and solid circle represent the 1σ morphologies of HESS J1702-420B and HESS J1702-420A, respectively. Finally, the red square -centered at the best-fit position of HESS J1702-420A -indicates the extraction region used to produce the profile reported in Figure H