© H. Abdalla et al. 2021

1. Introduction

The acceleration sites of cosmic rays¹ 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 (The KASCADE-Grande Collaboration 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 2016, 2018a), 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 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 2018b), 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 2018b) 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 2012, 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² 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° ×  0.15° 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 of 15σ based on 9.5 hr of observations. It simplified however the source morphology to a 0.2° 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° 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 Large Area Telescope (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 2018b). 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³⁵ erg s⁻¹, 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 2018c). 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 (Sect. 2), with focus on the three-dimensional (3D) likelihood anaysis (Sect. 2.1) and a morphological study made with a more classical background estimation technique (Sect. 2.2). Then, Sect. 3 reports on the analysis of archival multi wavelength observations of the region, while Sect. 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, Sect. 5 summarizes the main conclusions of the paper.

2. H.E.S.S. data analysis and results

H.E.S.S. is an array of five imaging atmospheric Cherenkov telescopes (IACTs) located in the Khomas Highland of Namibia, 1800 m above sea level. The design of the original array (H.E.S.S. I), operating since 2003, involved four 12 m diameter 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° from the HGPS position of HESS J1702-420 (l  =  344.30°, b  =   − 0.18°) and averaging over a 0.5° 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³ (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, 2020; Nigro et al. 2019). 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 Sect. 2.2 and Appendix A were also carried out using the standard HAP software.

2.1. 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 (Sect. 2.1.1) and results (Sect. 2.1.2), while Sect. 2.1.3 focuses on the most relevant components of the model, that is those describing HESS J1702-420.

2.1.1. 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 Sect. 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 time-saving, 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° ×  4° 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° ×  0.02° 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° from the pointing direction were excluded, together with those whose reconstructed energy was below the safe threshold

$$\begin{array}{c}\hfill E_{\mathrm{threshold}}=\max \{E_{\mathrm{bkg}},E_{A_{\mathrm{eff}}}\},\end{array}$$(1)

where E_bkg represents the energy at which the maximum rate of hadronic background events occurs and E_Aeff 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 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° band around the borders of the RoI was excluded from the analysis, in order to limit possible contamination due to non-modeled sources outside the RoI. Additionally, a 0.3°-radius circular region centered at l  =  343.35° and b  =   − 0.93°, 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.

2.1.2. 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° ×  4° 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 likelihood-ratio 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

$$\begin{array}{c}\hfill \mathrm{TS}=-2ln\left(\frac{{\mathcal{L}}_0^{\max }}{{\mathcal{L}}_{\mathrm{\mu }}^{\max }}\right).\end{array}$$(2)

Here, ${\mathcal{L}}_0^{\max }$ (${\mathcal{L}}_{\mathrm{\mu }}^{\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 Fig. 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 (Sect. 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 Fig. H.1.

Fig. 1. Upper left panel: image of the RoI obtained by integrating the binned cube of measured counts over the energy axis (E > 2 TeV), and correlating it with a 0.1°-radius top-hat kernel. The hatched regions were excluded from the likelihood computation. The bright area around l ≳ 345° results from deep observations of RX J1713-3946 and HESS J1708-410. Upper right panel: energy-integrated (E > 2 TeV) map of model-predicted counts, with names and 1σ shapes of all model components overlaid. The large-scale discarded component is indicated by the dashed circle – see the main text for more details. Lower left panel: spatial distribution of model residuals, showing the statistical significance – in units of Gaussian standard deviations – of counts – model fluctuations. The image was obtained assuming Cash statistic for Poisson-distributed signals with perfectly known background model (Cash 1979). Lower right panel: histogram containing the number of occurrences of each significance value (assuming Cash statistic), from the lower left panel. The adjustment of a Gaussian function to the histogram is shown, together with a reference standard normal distribution.

The upper right panel of Fig. 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 Sect. 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 Fig. 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–4.

2.1.3. 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 of (2.08 ± 0.49_stat ± 0.62_sys) × 10⁻¹³ cm⁻² s⁻¹ and a 1σ radius of (0.06 ± 0.02_stat ± 0.03_sys)°, 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⁻¹² cm⁻² s⁻¹ 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 2018b).

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 Fig. 2, together with spectral points computed under a power law assumption and re-optimizing 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.

Fig. 2. Power law spectra of HESS J1702-420A (red solid line) 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, re-optimizing 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.

We explicitly point out that, based on this dataset, it is impossible to tell whether HESS J1702-420A and HESS J1702-420B 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 Appendices D and A). Therefore, a model describing HESS J1702-420 with a single energy-dependent component is disfavored, even if it cannot be definitively ruled out.

2.2. 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 (Sect. 2.1). After subtracting the γ-like hadronic background, we measured γ-ray flux integrated above 2, 5, 15 and 40 TeV inside a 1.6° ×  1.6° region encompassing HESS J1702-420. The result is shown in Fig. 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 (Sect. 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)°.

Fig. 3. γ-ray flux maps of the HESS J1702-420 region, computed with the Ring Background Method, above 2 (top left), 5 (top right), 15 (bottom left) and 40 (bottom right) TeV. All maps are correlated with a 0.1°-radius top-hat kernel, and the color code is in unit of γ-ray flux (cm−2 s−1) per smoothing area. The white contours indicate the 3σ and 5σ H.E.S.S. significance levels (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.

3. 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 Sect. 3.2 reports our considerations on the surrounding ISM and Sect. 3.3 discusses archival Suzaku measurements in the context of the new H.E.S.S. results.

3.1. 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° ×  10° 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° 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 Fig. 4 shows that, within the source region, there is no evidence for a significant excess, but some low-significance fluctuations are present. For comparison, Fig. H.2 shows a TS map computed before removing the contribution from point sources.

Fig. 4. Residual TS map after source modeling in the RoI (see the main text for details). The yellow and white contours represent the TS = 4 (2σ) and TS = 9 (3σ) significance levels, respectively. The red ellipse (orange circle) correspond to the 1σ shape of HESS J1702-420B (HESS J1702-420A). The positions of nearby 4FGL and 2FHL sources are shown as green and cyan circles, respectively.

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 (Sect. 2.1.3). Its 1σ contour is indicated by the red ellipse in Fig. 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 with the HESS J1702-420B template shape, at the level of

$$\begin{array}{c}\hfill {\left(E^2\frac{\mathrm{d}N}{\mathrm{d}E}\right)}_{E=E_{\mathrm{ref}}}\le 7.6\times {10}^{-9}\mathrm{GeV}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1},\end{array}$$(3)

where E_ref ≈ 95 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 (Sect. 4.2).

3.2. The interstellar medium

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 λ = 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° ×  0.15° ellipse centered at l  =  344.30° and b  =   − 0.18°. 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 (Sect. 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 Fig. H.3), does not exhibit any obvious correlation with the VHE γ-ray maps (see Fig. H.4). In particular, no hydrogen cloud clearly correlates with HESS J1702-420A or HESS J1702-420B.

3.3. 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 point-like objects (src A and src B, indicated in Fig. 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 Fig. 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⁻¹⁴ erg s⁻¹ cm⁻², 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 Fig. 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.

4. 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 (Sect. 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 π⁰ decay developed in Kafexhiu et al. (2014), we assumed a fixed target density n_H = 100 cm⁻³. 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 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⁻³ and temperature of T_CMB = 2.73 K. The starlight emission in the near IR (ϵ_NIR = 1 eV cm⁻³ and T_NIR = 3000 K) and dust re-emission in the far IR (ϵ_FIR = 0.5 eV cm⁻³ 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 Sects. 4.1 and 4.2.

4.1. 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 Fig. 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.

Fig. 5. Models of γ-ray emission based on hadronic (red) and leptonic (blue) one-zone scenarios, for HESS J1702-420A – left panel – and HESS J1702-420B – right panel. The best-fit spectra, under the assumption of simple power law distribution of the underlying particle populations, are shown as solid lines, while the shaded areas and dotted lines represent the 1σ statistical error envelope and extrapolations outside the fit range, respectively. The H.E.S.S. and Fermi-LAT flux points are also shown, for reference purpose. Fit results were obtained with a 3D fit of H.E.S.S. data.

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) and Werner 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:

$$\begin{array}{c}\hfill W_{\mathrm{p}/\mathrm{e}}(E_{\mathrm{p}/\mathrm{e}}>1\mathrm{TeV})={\displaystyle {\int }_{1\mathrm{TeV}}^{\infty }{E}_{\mathrm{p}/\mathrm{e}}\frac{\mathrm{d}N}{\mathrm{d}{E}_{\mathrm{p}/\mathrm{e}}}\mathrm{d}{E}_{\mathrm{p}/\mathrm{e}}.}\end{array}$$(4)

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:

$$\begin{array}{c}\hfill W_{\mathrm{p}}(E_{\mathrm{p}}>1\mathrm{TeV})\gtrsim 1.8\times {10}^{47}{\left(\frac{d}{3.5\mathrm{kpc}}\right)}^2{\left(\frac{n_{\mathrm{H}}}{100{\mathrm{cm}}^{-3}}\right)}^{-1}\mathrm{erg}\end{array}$$(5)

$$\begin{array}{c}\hfill W_{\mathrm{e}}(E_{\mathrm{e}}>1\mathrm{TeV})\gtrsim 8.1\times {10}^{45}{\left(\frac{d}{3.5\mathrm{kpc}}\right)}^2\mathrm{erg}.\end{array}$$(6)

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 2018c). 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 Sect. 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 Fig. H.5). In other words, if the Suzaku measurement is reliable (see Sect. 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 PeVatron (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 Sect. 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 Sect. 3.2).

4.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

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}E}\propto \{\begin{array}{cc}{(E/E_0)}^{-{\alpha }_1},\hfill & \mathrm{if}E<\stackrel{\sim }{E}\hfill \\ {(\stackrel{\sim }{E}/E_0)}^{{\alpha }_2-{\alpha }_1}{(E/E_0)}^{-{\alpha }_2},\hfill & \mathrm{if}E>\stackrel{\sim }{E}\hfill \end{array}\end{array}$$(7)

where $\stackrel{\sim }{E}$ and E₀ 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 (Sect. 3.1). The first power law index, α₁, 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 α₁ = 1.6 (1.4) and α₂ = 2.66  ± 0.11_stat (3.39  ± 0.11_stat), and with break energy of $\stackrel{\sim }{E}=(6.77\pm 3.{64}_{\mathrm{stat}})\mathrm{TeV}$ ( (4.19  ± 1.25_stat) 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 γ-ray emission of HESS J1702-420B, were computed integrating the broken power law particle spectra above 1 GeV. They are:

$$\begin{array}{c}\hfill W_{\mathrm{p}}(E_{\mathrm{p}}>1\mathrm{GeV})\approx 2.8\times {10}^{48}{\left(\frac{d}{3.5\mathrm{kpc}}\right)}^2{\left(\frac{n_{\mathrm{H}}}{100{\mathrm{cm}}^{-3}}\right)}^{-1}\mathrm{erg}\end{array}$$(8)

$$\begin{array}{c}\hfill W_{\mathrm{e}}(E_{\mathrm{e}}>1\mathrm{GeV})\approx 4.5\times {10}^{47}{\left(\frac{d}{3.5\mathrm{kpc}}\right)}^2\mathrm{erg}.\end{array}$$(9)

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 Sect. 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;

$$\begin{array}{c}\hfill ϵ=\frac{L_{[1,10]\mathrm{TeV}}}{\dot{E}}\approx 19\%,\end{array}$$(10)

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 well above the efficiency of all other PWNe identified by H.E.S.S. in the same energy range (H.E.S.S. Collaboration 2018c). 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 2012, 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 γ-ray emission is detected by H.E.S.S. near the pulsar position – see Fig. 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 Sect. 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 Sect. 3.2).

4.3. 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 Fig. 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⁵⁰ 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⁻³, it has to be relatively close – d ≲ 2 kpc –, unless its proton energy budget exceeds 10⁵⁰ 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.

Fig. 6. Possible constraints on the gas density, magnetic field and distance of HESS J1702-420, based on multi-wavelegth observations, under the assumption of simple one-zone hadronic (left panel) or leptonic (right panel) scenarios. More details are given in the main text.

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 Fig. 6 can be found in Appendix G.

5. 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.

Acknowledgments

The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of H.E.S.S. is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the German Research Foundation (DFG), the Helmholtz Association, the Alexander von Humboldt Foundation, the French Ministry of Higher Education, Research and Innovation, the Centre National de la Recherche Scientifique (CNRS/IN2P3 and CNRS/INSU), the Commissariat à l’énergie atomique et aux énergies alternatives (CEA), the UK Science and Technology Facilities Council (STFC), the Knut and Alice Wallenberg Foundation, the National Science Centre, Poland grant no. 2016/22/M/ST9/00382, the South African Department of Science and Technology and National Research Foundation, the University of Namibia, the National Commission on Research, Science and Technology of Namibia (NCRST), the Austrian Federal Ministry of Education, Science and Research and the Austrian Science Fund (FWF), the Australian Research Council (ARC), the Japan Society for the Promotion of Science and by the University of Amsterdam. We appreciate the excellent work of the technical support staff in Berlin, Zeuthen, Heidelberg, Palaiseau, Paris, Saclay, Tubingen and in Namibia in the construction and operation of the equipment. This work benefited from services provided by the H.E.S.S. Virtual Organisation, supported by the national resource providers of the EGI Federation.

References

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°-radius circle and three 0.2°-radius annuli were used to measure the VHE γ-ray spectrum of the source. We did not adopt narrower extraction regions, in order 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° and b  =   − 0.15°, 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.

Fig. A.1. Upper panel: Map of the H.E.S.S. γ-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.

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 Fig. 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 (2018b). 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 2019).

Appendix B: Hypothesis testing for nested parametric models

According to Wilks’ theorem (Wilks 1938), the TS defined in Eq. (2) is distributed as a ${\chi }_{\nu }^2$, where ν is the number of additional degrees of freedom of the alternative hypothesis with respect to the null hypothesis. The theorem is valid under the assumptions – always satisfied in our analysis – of high statistics and nested models. Thanks to this theorem, the statistical significance of the alternative hypothesis can be directly estimated from the TS value, by determining the corresponding right-tail p-value of a ${\chi }_{\nu }^2$ distribution. To convert the significance into units of Gaussian standard deviations (σ), it is then sufficient to compute z-score of a Gaussian distribution corresponding to the given p-value.

Appendix C: Non-confirmed large-scale emission component

In the main analysis, a large-scale – ≈0.5° in radius – model component was detected around l  =  345.23° and b  =   − 0.01°. At that position, the borders of several runs partially overlap, resulting in a boosted exposure level but also strong systematics due to edge effects. In the crosscheck analysis, the exposure level at the position of the large-scale component is lower, resulting in a decreased sensitivity. Accordingly, we verified that its inclusion or exclusion in the source model of the crosscheck analysis did not have any relevant impact on the predicted number of counts at its position. We therefore could not confirm the detection of this new large-scale emission component nearby HESS J1702-420. We point out that the HGPS already reported the presence of a large-scale component with similar position and size, called HGPSG 041, that was similarly discarded due to a non-detection in the crosscheck analysis. In the future, new dedicated observations of the region with more uniform exposure will ultimately probe its presence and nature.

Appendix D: 3D analysis of H.E.S.S. data in independent energy bands

According to the 3D model developed in Sect. 2.1.2, HESS J1702-420 is best described by the superposition of two independent components. To further validate this fact, we repeated the 3D analysis within three independent – that is nonoverlapping – energy bands defined by the edges 2.0, 3.7, 15.3 and 150 TeV. These were chosen to ensure a roughly constant level of γ-ray flux. During the fit, the spectral indices of all model components were fixed to the values obtained in the whole energy range (Sect. 2.1.2). This was meant to prevent poor spectral modeling, due to the limited lever arm and insufficient number of photons reconstructed within each individual energy band. The spectral normalizations were instead left free to vary, together with all spatial and background parameters. For each energy band, we used the likelihood ratio test (see Eq. (2)) to compare the statistical significance of three nested hypotheses:

• H₀

  Null hypothesis, with no model component describing HESS J1702-420;

• H₁

  HESS J1702-420 is described by one Gaussian component, with the spectral index of HESS J1702-420B. For this component, we left the spatial eccentricity and rotation angle free to vary, for a total of six⁴ free model parameters;

• H₂

  HESS J1702-420 is described by two Gaussian components, with the spectral indices of HESS J1702-420B and HESS J1702-420A. For the latter, we considered a strictly symmetric Gaussian morphology, for a total of four⁵ free model parameters.

The relative significance of each hypothesis, for all the energy bands, is reported in Table D.1. It turns out that HESS J1702-420B is significant in all energy bands, while HESS J1702-420A is significant only in the 15.3−150 TeV band.

The spatial and spectral shapes of the two components within each energy band are shown in Fig. D.1, where the reference result from the whole 2−150 TeV fit range are reported in black. The upper panel shows the 1σ contours of HESS J1702-420A and HESS J1702-420B. HESS J1702-420A is not drawn in the two lowest energy bands (i.e., for E ≤ 15.3 TeV), because it is not significant (see Table D.1). In each energy band, the best-fit morphologies of both components are consistent with the reference results obtained in Sect. 2.1. The lower panel of the figure compares the reference spectra of HESS J1702-420A and HESS J1702-420B (in black) with the spectra obtained in different energy bands. The energy-resolved results are well connected and in agreement with the reference power laws.

Fig. D.1. Upper panel: the figure, centered at the approximate position of HESS J1702-420, shows contours corresponding to 150, 200, 250 and 300 counts detected by H.E.S.S. above 2 TeV per smoothing area. Overlaid on the map are the 1σ extension contours of the components HESS J1702-420A and HESS J1702-420B, as obtained from the 3D fit in separate energy bands. Lower panel: spectral results of the energy-resolved 3D analysis, for HESS J1702-420A and HESS J1702-420B. Vertical lines separate the energy bands that were independently used to perform the source modeling. In both panels, the reference results obtained over the full energy range (see Sect. 2.1.3) are indicated in black.

To summarize, in each independent energy band HESS J1702-420 is well described by a simple model, based on either one or two components with Gaussian morphologies and power law spectra. HESS J1702-420B is significant in all energy bands, with stable morphology and spectrum. HESS J1702-420A instead is significant only in the highest energy band, due to its exceptionally hard spectrum. This is precisely what would be expected in case the emission is due to two separate components. Additionally, we summed the log-likelihood values obtained in each independent energy band, and estimated that globally a two-component model is better than a one-component model with a confidence level of 5.3σ. This value is consistent with the 5.4σ significance for the presence of HESS J1702-420A obtained from the 3D analysis in the full 2−150 TeV energy range (see Table 3). All these facts support the validity of the simple two-component approach that resulted naturally from the iterative procedure described in Sect. 2.1.2. However, the possibility that a single source component with energy-dependent morphology could provide a better fit of the data – at the expense of a large number of free parameters describing the variation of source size, eccentricity and center position as a function of the energy – cannot be ruled out at this stage.

Appendix E: Fermi-LAT analysis details

The data (photon event file and spacecraft file) were retrieved from the LAT data server, through a query defined by the parameters in Table E.1. We adopted the event selection cuts described in Table E.2. For the analysis, we defined a square 10° ×  10° RoI, fully inscribed within the events selection circle. Events were binned spatially using 0.05° ×  0.05° spatial pixels, and spectrally using 8 bins per energy decade.

During the maximum likelihood fit, the spectral index and normalization of all sources within 3° from the RoI center and having a TS value higher than 25 were left free to vary. Additionally, the spectral normalization of all sources with TS > 30 within the whole 10° ×  10° was also adjusted. The Galactic diffuse emission model (gll_iem_v07.fits) was left free to vary, while the extra-Galactic diffuse model was considered fixed to the default one (iso_P8R3_SOURCE_V2_v1.txt).

Appendix F: Method for the derivation of lower limits on the particle cut-off energy

The lower limits on the cut-off energies in the spectra of the proton or electron parent populations were obtained as follows:

• (i)

  we defined an array of trials particle cut-off energies $\{E_1^{\text{c}},\dots ,E_N^{\text{c}}\}$;

• (ii)

  for each fixed cut-off energy $E_i^{\text{c}}$, we adjusted a power-law with exponential cut-off to the 3D H.E.S.S. data. The spectral normalization and index, together with all free nuisance parameters of the model, were optimized at each step. For each trial cut-off energy, we stored the likelihood value of the fit, ${\mathcal{L}}_{(E_i^{\text{c}})}^{\text{\hspace{0.05em}max}}$;

• (iii)

  then we computed the profile of

  $$\begin{array}{c}\hfill \mathrm{TS}(E_i^{\mathrm{c}})=-2ln\left(\frac{{\mathcal{L}}_0^{\max }}{{\mathcal{L}}_{(E_i^{\mathrm{c}})}^{\max }}\right),\end{array}$$(F.1)

  where ${\mathcal{L}}_0^{\max }$ 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°, 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:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{diff}}\approx \frac{{[R_{\mathrm{source}}(d)]}^2}{6D(E,B)}.\end{array}$$(G.1)

Energy losses for protons due to p–p collisions were estimated – neglecting ionization losses that are irrelevant for relativistic protons – as

$$\begin{array}{c}\hfill {\tau }_{\mathit{pp}}\approx 6\times {10}^5{\left(\frac{n_{\mathrm{H}}}{100{\mathrm{cm}}^{-3}}\right)}^{-1}\mathrm{yr},\end{array}$$(G.2)

as in Gabici et al. (2007). Finally, the electron energy loss timescale was computed using

$$\begin{array}{c}\hfill {\tau }_{\mathrm{loss}}={(\frac{1}{{\tau }_{\mathrm{syn}}}+{\displaystyle \sum _i}\frac{1}{{\tau }_{\mathrm{IC}}^i})}^{-1},\end{array}$$(G.3)

where

$$\begin{array}{c}\hfill {\tau }_{\mathrm{syn}}\approx 1.3\times {10}^5{\left(\frac{E}{1\mathrm{TeV}}\right)}^{-1}{\left(\frac{B}{10\mathrm{\mu }\mathrm{G}}\right)}^{-2}\mathrm{yr}\end{array}$$(G.4)

is the synchrotron loss timescale and

$$\begin{array}{c}\hfill {\tau }_{\mathrm{IC}}\approx 3\times {10}^7{\left(\frac{E}{10\mathrm{GeV}}\right)}^{-1}{\left(\frac{U_{\mathrm{rad}}}{1\mathrm{eV}{\mathrm{cm}}^{-3}}\right)}^{-1}\mathrm{yr}\end{array}$$(G.5)

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.

Fig. H.1. Spatial distributions of the significance of model residuals (left column) and histograms of significance values (right column), computed in the energy bands 2.0−5.0 TeV (first row), 5.0−15.0 TeV (second row) and 15.0−150 TeV (third row).

Fig. H.2. Map of test statistic (TS) for the presence of an additional source in the RoI, with respect to a source model containing only the galactic and isotropic γ-ray diffuse sources. The color bar, contours and markers are identical to Fig. 4.

Fig. H.3. Brightness temperature peaks obtained by integrating the data from the Mopra radio survey (Braiding et al. 2018) within a 0.3°-side square window centered at the best-fit position of HESS J1702-420A (shown in red in Fig. H.4).

Fig. H.4. Column density maps of molecular hydrogen in the direction of HESS J1702-420, obtained by integrating the brightness temperature profile of 12CO(J = 1 → 0) data from the Mopra radio survey within the velocity intervals indicated above each panel (corresponding to the peaks in Fig. H.3). The brightness temperature values were converted to H2 column density assuming the conversion factor XCO = 1.5 × 1020 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 Fig. H.3.

Fig. H.5. multi wavelength modeling of HESS J1702-420A, under the assumption of a one-zone leptonic scenario powered by a power law distribution of electrons. The synchrotron emission was computed assuming a magnetic field value of 0.3 and 3 μG. The black error bars represent the VHE flux points of HESS J1702-420A, while the red one indicates the flux of the unidentified Suzaku src B (Fujinaga et al. 2011).

All Tables

All Figures

Fig. 1. Upper left panel: image of the RoI obtained by integrating the binned cube of measured counts over the energy axis (E > 2 TeV), and correlating it with a 0.1°-radius top-hat kernel. The hatched regions were excluded from the likelihood computation. The bright area around l ≳ 345° results from deep observations of RX J1713-3946 and HESS J1708-410. Upper right panel: energy-integrated (E > 2 TeV) map of model-predicted counts, with names and 1σ shapes of all model components overlaid. The large-scale discarded component is indicated by the dashed circle – see the main text for more details. Lower left panel: spatial distribution of model residuals, showing the statistical significance – in units of Gaussian standard deviations – of counts – model fluctuations. The image was obtained assuming Cash statistic for Poisson-distributed signals with perfectly known background model (Cash 1979). Lower right panel: histogram containing the number of occurrences of each significance value (assuming Cash statistic), from the lower left panel. The adjustment of a Gaussian function to the histogram is shown, together with a reference standard normal distribution.

In the text

Fig. 2. Power law spectra of HESS J1702-420A (red solid line) 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, re-optimizing 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.

In the text

Fig. 3. γ-ray flux maps of the HESS J1702-420 region, computed with the Ring Background Method, above 2 (top left), 5 (top right), 15 (bottom left) and 40 (bottom right) TeV. All maps are correlated with a 0.1°-radius top-hat kernel, and the color code is in unit of γ-ray flux (cm−2 s−1) per smoothing area. The white contours indicate the 3σ and 5σ H.E.S.S. significance levels (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.

In the text

	Fig. 4. Residual TS map after source modeling in the RoI (see the main text for details). The yellow and white contours represent the TS = 4 (2σ) and TS = 9 (3σ) significance levels, respectively. The red ellipse (orange circle) correspond to the 1σ shape of HESS J1702-420B (HESS J1702-420A). The positions of nearby 4FGL and 2FHL sources are shown as green and cyan circles, respectively.
In the text

Fig. 5. Models of γ-ray emission based on hadronic (red) and leptonic (blue) one-zone scenarios, for HESS J1702-420A – left panel – and HESS J1702-420B – right panel. The best-fit spectra, under the assumption of simple power law distribution of the underlying particle populations, are shown as solid lines, while the shaded areas and dotted lines represent the 1σ statistical error envelope and extrapolations outside the fit range, respectively. The H.E.S.S. and Fermi-LAT flux points are also shown, for reference purpose. Fit results were obtained with a 3D fit of H.E.S.S. data.

In the text

	Fig. 6. Possible constraints on the gas density, magnetic field and distance of HESS J1702-420, based on multi-wavelegth observations, under the assumption of simple one-zone hadronic (left panel) or leptonic (right panel) scenarios. More details are given in the main text.
In the text

Fig. A.1. Upper panel: Map of the H.E.S.S. γ-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.

In the text

Fig. D.1. Upper panel: the figure, centered at the approximate position of HESS J1702-420, shows contours corresponding to 150, 200, 250 and 300 counts detected by H.E.S.S. above 2 TeV per smoothing area. Overlaid on the map are the 1σ extension contours of the components HESS J1702-420A and HESS J1702-420B, as obtained from the 3D fit in separate energy bands. Lower panel: spectral results of the energy-resolved 3D analysis, for HESS J1702-420A and HESS J1702-420B. Vertical lines separate the energy bands that were independently used to perform the source modeling. In both panels, the reference results obtained over the full energy range (see Sect. 2.1.3) are indicated in black.

In the text

	Fig. H.1. Spatial distributions of the significance of model residuals (left column) and histograms of significance values (right column), computed in the energy bands 2.0−5.0 TeV (first row), 5.0−15.0 TeV (second row) and 15.0−150 TeV (third row).
In the text

	Fig. H.2. Map of test statistic (TS) for the presence of an additional source in the RoI, with respect to a source model containing only the galactic and isotropic γ-ray diffuse sources. The color bar, contours and markers are identical to Fig. 4.
In the text

	Fig. H.3. Brightness temperature peaks obtained by integrating the data from the Mopra radio survey (Braiding et al. 2018) within a 0.3°-side square window centered at the best-fit position of HESS J1702-420A (shown in red in Fig. H.4).
In the text

Fig. H.4. Column density maps of molecular hydrogen in the direction of HESS J1702-420, obtained by integrating the brightness temperature profile of 12CO(J = 1 → 0) data from the Mopra radio survey within the velocity intervals indicated above each panel (corresponding to the peaks in Fig. H.3). The brightness temperature values were converted to H2 column density assuming the conversion factor XCO = 1.5 × 1020 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 Fig. H.3.

In the text

Fig. H.5. multi wavelength modeling of HESS J1702-420A, under the assumption of a one-zone leptonic scenario powered by a power law distribution of electrons. The synchrotron emission was computed assuming a magnetic field value of 0.3 and 3 μG. The black error bars represent the VHE flux points of HESS J1702-420A, while the red one indicates the flux of the unidentified Suzaku src B (Fujinaga et al. 2011).

In the text