Multiwavelength constraints on the unidentified Galactic TeV sources HESS J1427$-$608, HESS J1458$-$608, and new VHE $\gamma$-ray source candidates

The H.E.S.S. Galactic Plane Survey (HGPS) revealed 78 TeV sources among which 47 are not clearly associated with a known object. We present a multiwavelength approach to constrain the origin of the emission from unidentified HGPS sources. We present a generic pipeline that explores a large database of multiwavelength archival data toward any region in the Galactic plane. Along with a visual inspection of the retrieved multiwavelength observations to search for faint and uncataloged counterparts, we derive a radio spectral index that helps disentangle thermal from nonthermal emission and a mean magnetic field through X-ray and TeV data in case of a leptonic scenario. We also search for a spectral connection between the GeV and the TeV regimes with the Fermi-LAT cataloged sources that may be associated with the unidentified HGPS source. We complete the association procedure with catalogs of known objects and with the source catalogs from instruments whose data are retrieved. The method is applied on two unidentified sources, namely HESS J1427$-$608 and HESS J1458$-$608, for which the multiwavelength constraints favor the pulsar wind nebula (PWN) scenario. We model their broadband nonthermal spectra in a leptonic scenario with a magnetic field $B \lesssim 10$ $\mu$G, which is consistent with that obtained from ancient PWNe. We place both sources within the context of the TeV PWN population to estimate the spin-down power and the characteristic age of the putative pulsar. We also shed light on two possibly significant $\gamma$-ray excesses in the HGPS: the first is located in the south of the unidentified source HESS J1632$-$478 and the second is spatially coincident with the synchrotron-emitting supernova remnant G28.6$-$0.1. The multiwavelength counterparts found toward both $\gamma$-ray excesses make these promising candidates for being new very-high energy $\gamma$-ray sources.


Introduction
Supernova remnants (SNRs) and pulsar wind nebulae (PWNe) are considered the best candidates to accelerate the bulk of Galactic cosmic rays (CRs) at least up to the knee of the CR spectrum (∼ 3 × 10 15 eV). Although no evidence for particle acceleration up to PeV energies has been found so far, efficient particle acceleration in SNRs and PWNe has been observed and can be explained by the widely adopted diffusive shock acceleration mechanism (DSA; Bell 1978). Accelerated electrons and positrons radiate synchrotron emission and also γ rays through Bremsstrahlung and inverse Compton (IC) scattering on photon fields while accelerated protons and heavier nuclei colliding with ambient matter emit γ rays through neutral pion decay. Thus, γray studies probe the population of both electrons and protons and allow us to disentangle different origins of the emission. However, γ-ray instruments usually have lower angular resolution than most radio and X-ray instruments, and their observations often suffer from source confusion. Thus, probing the accelerated leptons through the expected synchrotron emission mainly in the radio and X-ray bands helps constrain the nature of γ-ray sources.
The identification of γ-ray emitting sources in the highenergy (HE; 0.1 < E < 100 GeV) and very-high energy (VHE; 0.1 TeV < E < 100 TeV) domains often relies on spatial correlation with identified objects at other wavelengths. While catalogs of known γ-ray emitters in the Galaxy such as SNRs, PWNe, and pulsars can provide hints of the origin of the emission, multiwavelength data exploration is often necessary to pinpoint the nature of the HE/VHE sources (Aliu et al. 2014;MAGIC Collaboration et al. 2014. For instance, the TeV shell HESS J1534−571 was found to be spatially coincident with the SNR candidate G323.7−1.0 (discovered through a close inspection of archival radio images) and was thus identified as the VHE counterpart of this SNR (H. E. S. S. Collaboration et al. 2018a). One of the brightest TeV sources, HESS J1825−137, was also identified as a PWN owing to its association with the X-ray PWN G18.0−0.7 and its γ-ray energy-dependent morphology shrinking to the energetic pulsar PSR J1826−1334 at the highest energies (Aharonian et al. 2005). Finally, the source Article number, page 1 of 19 arXiv:2101.07775v1 [astro-ph.HE] 19 Jan 2021 HESS J1356−645 was classified as an evolved PWN as a result of the exploitation of multiwavelength data that revealed a nonthermal, center-filled radio and X-ray emission surrounding the energetic pulsar PSR J1357−6429 and spatially coincident with the H.E.S.S. source (H.E.S.S. Collaboration et al. 2011).
Several databases have been developed to gather multiwavelength observations that can be easily retrieved. The Centre d'Analyse des Données Étendues 1 and the Space Science Data Center 2 provide interactive databases, such as the Skyview 3 tool, which enables the production of images at multiple wavelengths for any position in the sky. The interactive tool gamma sky 4 permits us to superimpose catalogs such as those of the Fermi-LAT and known SNRs.
However, using multiwavelength data to identify γ-ray sources is not always straightforward. The best example is the nearby and evolved PWN Vela X, which exhibits a changing morphology at different wavelengths. The brightest part of the VHE emission is spatially coincident with the X-ray PWN, but its largest extent is consistent with that of the surrounding extended radio nebula of 2 • × 3 • (Abramowski et al. 2012). Emission in the north of the pulsar is also seen with the Fermi-LAT (Grondin et al. 2013), illustrating how the interpretation of multiwavelength data can be difficult and puzzling. Moreover, multiwavelength associations relying on spatial correlations is challenging if the γ-ray emission is produced by an ancient PWN or a molecular cloud (MC) illuminated by CRs that have escaped from a nearby SNR. In the former case, the pulsar can be largely offset from the TeV emission produced by old electrons (de Jager et al. 2009) and, in the latter case, the MC can also be distanced from the SNR from which CRs have escaped (Gabici et al. 2009). Thus, there may be no spatial correlation between the driven-phenomenon of the emission (here pulsar or SNR) and the TeV emission itself, making the association with multiwavelength data difficult. Finally, nonthermal radio and Xray emission from ancient PWNe and MCs illuminated by CRs is expected to be faint, challenging the identification process.
With 2700 hours of observations, the H.E.S.S. experiment, through a Galactic Plane Survey (HGPS; H. E. S. S. Collaboration et al. 2018b) covering the inner part of the Galaxy (l = 250 • to 65 • and |b|< 3.5 • ), has reported 78 TeV sources, among which 12 PWNe, 16 SNRs (or composite SNRs), 3 binaries, and 47 sources that are not firmly identified. With a point spread function (PSF) of ∼ 0.08 • (68% containment radius) and a pointsource sensitivity of 1.5% of the Crab Nebula VHE flux, this Galactic scan was performed from 0.2 to 100 TeV with an unprecedented spatial coverage. Besides positional evidence with cataloged objects, the firm identification of the VHE sources was based on correlated multiwavelength variability, matching multiwavelength morphology and energy-dependent γ-ray morphology. The association procedure made use of the 3FGL and 2FHL Fermi-LAT catalogs (Acero et al. 2015;Ackermann et al. 2016), the "SNRcat" (Ferrand & Safi-Harb 2012), the ATNF pulsar catalog (Manchester et al. 2005, version 1.54), and 20 external analyses involving data at other wavelengths. When a firm identification was not possible, the H.E.S.S. sources were associated with all objects whose cataloged position is at an angular distance smaller than the H.E.S.S. source spectral extraction radius (noted R spec ). For a total of 47 unidentified sources, the HGPS association procedure reported 11 sources that are not associated with any cataloged sources and 36 sources for which the origin of the emission is unclear, mainly owing to several possible scenarios and source confusion.
Given the large number of TeV sources without firm identification, we developed a generic pipeline that retrieves archival multiwavelength data toward any region of the Galactic plane to constrain the origin of the γ-ray emission. Although extragalactic sources can be located near the Galactic plane (as exemplified by the blazar HESS J1943+213, H. E. S. S. Collaboration et al. 2018b, and references therein), we are primarily interested in Galactic sources that are either extended at TeV energies or associated with an extended component seen at other wavelengths. Thus, in this study we also discarded high-mass binary systems, some of which are known to be HE and VHE emitters such as PSR B1259−63 (e.g., Chernyakova & Malyshev 2020). We first describe our multiwavelength approach in Section 2. In Section 3, we apply this method on two unidentified sources, namely HESS J1427−608 and HESS J1458−608, and place them within the context of the PWN population seen at TeV energies. Finally, Section 4 sheds light on possibly significant (but not cataloged) TeV excesses in two confused regions: in the south of HESS J1632−478 and toward the SNR G28.6−0.1, which both have multiwavelength counterparts, strengthening the interest in carrying out dedicated γ-ray analyses to confirm these as new VHE sources.

Multiwavelength approach
We present in this section a generic pipeline that retrieves archival radio, X-ray, infrared, and GeV data to search for uncataloged, presumably faint counterparts of the TeV unidentified sources. The association procedure is made with the instrument source catalogs whose data are retrieved and with catalogs of known objects that are more numerous than those used for the HGPS association procedure. After automatically downloading all the archival data and performing the association procedure, we took prior attention to the visual inspection of all the archival multiwavelength images to search for faint excesses within the TeV source extent and to define the regions of interest for flux extraction. This is the only step that needs manual intervention. After this, the pipeline automatically derives a radio spectral index that allows us to disentangle thermal from nonthermal emission and a mean magnetic field using X-ray and TeV data, assuming that the γ-ray emission is produced by IC scattering off the cosmic microwave background (CMB) photon field. Finally, the pipeline automatically plots the spectra of the Fermi-LAT cataloged sources found within the TeV source extent together with that of the H.E.S.S. source to search for a spectral connection between the HE and VHE regimes.

Archival data retrieval
Given a position in the sky and a search radius, we developed a code that automatically extracts radio continuum data from 13 single-dish telescopes or interferometers that survey the northern and the southern sky. This includes the second epoch of the Molongo Galactic Plane Survey (MGPS-2, 843 MHz; Murphy et al. 2007), the VLA Galactic Plane Survey (VGPS, 1.4 GHz; Stil et al. 2006), the Southern Galactic Plane Survey (SGPS, 1.4 GHz, a combination of the Australia Telescope Compact Array and the Parkes Radio Telescope; McClure-Griffiths et al. 2005), the TIFR GMRT Sky Survey (TGSS, 150 MHz;Intema et al. 2017), the NRAO-VLA Sky Survey (NVSS, 1.4 GHz;Condon et al. 1998), the Canadian Galactic Plane Survey (CGPS, 1.4 GHz;Taylor et al. 2003), the Green Bank Telescope (8.35 at 14.35 GHz and the 87GB survey at 4.85 GHz; Langston et al. 2000;Gregory & Condon 1991), and the HI/OH/Recombination line survey of the inner Milky Way (THOR, 1, 1.3, 1.4, 1.6, and 1.8 GHz; Beuther et al. 2016). We also used data from the Parkes 64 m radio telescope (2.4 GHz; Duncan et al. 1995) and those from the Parkes-MIT-NRAO survey 1 (PMN, 4.85 GHz;Griffith & Wright 1993). Archival data from the HI Parkes All-Sky Survey (HIPASS) and the HI Zone of Avoidance (HIZOA) survey were reprocessed into a new continuum map that we retrieved (CHIPASS, 1.4 GHz; Calabretta et al. 2014). We also retrieved data from the Multi-Array Galactic Plane Imaging Survey (MAGPIS, 0.325, 1.5 and 5 GHz;Helfand et al. 2006). The technical details of the 13 radio instruments and data retrieval are given in Appendix A (and references therein). In addition, data from the Planck satellite are also used (30, 44, 70, 100, and 143 GHz, Planck Collaboration et al. 2016).
For the association procedure, we used catalogs of objects reporting pulsars, PWNe, SNRs, H ii regions, etc. This includes the ATNF pulsar catalog (Manchester et al. 2005, version 1.58), the "SNRcat" (Ferrand & Safi-Harb 2012), the Galactic SNR catalog (Green 2017), and the catalog of 76 Galactic SNR candidates revealed with THOR data (Anderson et al. 2017). We used Fermi-LAT pulsar, point-source, and extended-source catalogs (2PC, 3FGL, 2FHL, 3FHL, FGES, and 4FGL catalogs 2 ) as well as the HAWC catalog (2HWC; Abeysekara et al. 2017). We note that energetic pulsars are considered indirect associations since we do not expect them to produce TeV emission but rather to generate a PWN that could be seen at TeV energies. We also used the catalogs of H ii regions obtained with WISE data (Anderson et al. 2014), of MCs in the Milky Way (Rice et al. 2016), of Galactic O stars (Maíz Apellániz et al. 2013, and information reported in the TeV source catalog 3 . The details of the catalogs are given in Appendix A. To reduce potential missing associations, we automatically requested the Simbad 4 database, which provides information on the listed astronomical objects. 1 This survey was made using the Parkes 64 m radio telescope with the NRAO multibeam receiver at a frequency of 4.85 GHz.

Nonthermal radio emission
When extended radio emission is found toward the unidentified TeV source, we calculated a radio spectral index α (defined as S ν ∝ ν α , where S ν and ν are the flux density and frequency, respectively) that helps identify thermal from nonthermal (i.e., synchrotron) emission. While α 0 in thermal sources such as H ii regions, its average value amounts to ∼ −0.4/−0.5 in SNRs (Green 2017) and ∼ −0.3/0 in PWNe (de Jager et al. 2009). Although radio emission from PWNe and H ii regions can exhibit a similar spectral index, they differ in shape with a shell-like and a Gaussian-like morphology for H ii regions and PWNe, respectively. In each archival image, we masked the cataloged radio sources to estimate the noise and we calculated the flux in a given region of interest (called R ON ) by summing up the flux in the pixels and correcting from the instrument beam. The flux extraction region was chosen after a thorough visual inspection of the data. The background-corrected fluxes, calculated at different frequencies, were fit with a power law to derive the radio spectral index. Upper limits on flux were derived at the 3σ confidence level. We followed the same method as for the PWN HESS J1356−645 (H.E.S.S. Collaboration et al. 2011), which is now implemented in a generic way using all the retrieved radio observations. Detailed explanations on the radio spectral index derivation, with an illustration on the PWN HESS J1356−645, are given in Appendix B.

Mean magnetic field
Nonthermal X-ray emission probes the presence of the highestenergy electrons that also radiate TeV γ rays through IC scattering on photon fields. Although ROSAT/PSPC detects soft X-ray photons (up to 2.4 keV), which can be heavily absorbed by the interstellar gas, we took advantage of its large field of view (∼ 6 • ×6 • ) to derive constraints on the X-ray flux within any spectral extraction region defined in the HGPS. We masked the cataloged sources (reported in the 1RXS and 2RXS catalogs; Voges et al. 1999;Boller et al. 2016) to estimate the background with the ring and reflected background methods (Berge et al. 2007). The significance is calculated following Li & Ma (1983), taking a correlation radius equal to the spectral extraction radius (R spec ) used in the HGPS. For a given X-ray spectral index Γ X and a column density N H , the source count rate (or its 5σ upper limit 5 ) in the 0.9−2.4 keV band is converted into a flux value (or upper limit) using the HEASARC tool WebPimms 6 . Assuming a one-zone model and that the TeV emission is produced by IC scattering on the CMB in the Thomson regime (Γ X = Γ TeV = Γ), the X-ray and TeV flux measurements (F sync , F IC ) can be used to constrain the mean magnetic field, which is expressed as where Γ = (p + 1)/2 is the photon spectral index and p is the particle spectral index. More details on the mean magnetic field estimate are given in Appendix C, with an illustration on the SNR RX J1713.7−3946. Although nonthermal X-ray filaments detected in several SNRs imply a magnetic field amplification up to hundreds of µG (Parizot et al. 2006), the expected mean magnetic field value in case of linear DSA is roughly B ∼ 15−20 µG for SNRs (Renaud 2009) and B ∼ 3−8 µG for evolved PWNe (Torres et al. 2014).

Gamma-ray spectra
The GeV spectra show different features depending on the origin of the emission. γ-ray emission from pulsars is usually described by a power law with an exponential cutoff occurring below tens of GeV or by a logarithmic parabola. The second Fermi-LAT pulsar catalog (2PC; Abdo et al. 2013) contains 117 pulsars for 2796 radio-detected pulsars. This difference can be explained by the pulsar spin-down powerĖ, which needs to be high enough to provide efficient particle acceleration up to γ-ray emitting energies. The emission from PWNe is usually characterized by a power law with a hard spectral shape (Γ < 2, where Γ is the photon spectral index) originating from IC scattering on photon fields. This hard spectral shape limits the detection of PWNe with the Fermi-LAT and makes PWNe the most numerous VHEemitting objects in our Galaxy. While a hard spectral shape is indicative of a leptonic origin of the emission, a roll-off in the photon spectrum below 100−200 MeV is the signature of accelerated protons colliding with ambient matter. In this case, the photon spectrum at GeV energies is usually found to be soft (Γ 2) in SNRs interacting with MCs such as IC 443, W44, or W51C (Jogler & Funk 2016, and references therein). Given the crucial information brought by GeV spectra, we retrieved the Fermi-LAT spectrum of all cataloged sources located within the R spec of the HGPS unidentified source to search for a spectral connection between GeV and TeV energies.  Aharonian et al. (2008) derived a size of σ = 0.063 • ± 0.01 • while the HGPS reported σ = 0.048 • ± 0.009 • , a value just below the significant extension threshold. The upper limit at the 95% confidence level on the Gaussian extent is σ = 0.067 • . Taking a spectral extraction radius R spec = 0.15 • , the spectrum is represented by a power law with a spectral index Γ TeV = 2.85 ± 0.22 and an integrated energy flux F 1−10 TeV = (1.43 ± 0.34) × 10 −12 erg cm −2 s −1 . The source is associated with an extended nonthermal X-ray emission (Suzaku J1427−6051, σ = 0.9 ± 0.1 ) with an energy flux F X = 8.9 +3.6 −2.0 × 10 −13 erg cm −2 s −1 and a spectral index Γ X = 3.1 +0.6 −0.5 (Fujinaga et al. 2013). Such a soft spectrum indicates that the photon cutoff energy is likely below the Suzaku energy band. The derived column density N H = 1.1 +0.29 −0.25 ×10 23 cm −2 points toward a possible large distance that could explain why it appears as point-like with H.E.S.S. No associated X-ray point source with similar column density was reported using XMM-Newton data (Fujinaga et al. 2013). Above 3 GeV, a Fermi-LAT point source was detected at the center of HESS J1427−608 with a hard spectral index Γ GeV = 1.87 ± 0.17 and an integrated flux of Φ GeV = (3.06 ± 0.73) × 10 −10 photon cm −2 s −1 (Guo et al. 2017). Figure 1 (left) shows the HGPS significance map with the associated X-ray and GeV sources overlaid. The extended X-ray emission and the hard spectral shape of the Fermi-LAT source, both spatially coincident with the TeV emission, make the leptonic scenario plausible as the origin of the TeV emission.

Multiwavelength constraints
We extracted archival radio continuum data toward the region of HESS J1427−608. Figure 1 (right) shows the MGPS-2 map with the extended nonthermal emission detected by Suzaku overlaid. Two MGPS-2 sources (MGPS J142755−605038 and MGPS J142714−604755) are located inside the R spec of HESS J1427−608 and a H ii region (reported in the PMN catalog) lies in the northwest (outside) of the H.E.S.S. source extent. The compact source MGPS J142755−605038 (F = 34.5 ± 6.4 mJy) was considered as the possible radio counterpart by Vorster et al. (2013). While the authors could reproduce the X-ray and TeV data with a leptonic model implying B = 4 µG, the predicted radio flux was significantly larger than that of MGPS J142755−605038. On the other hand, they could explain the radio and TeV data with B = 0.4 µG, but the model fails to account for the flux measured with Suzaku. We defined the ON region such as it encompasses the brightest part of the emission (located between the two MGPS-2 sources) and MGPS J142755−605038 (lying at the center of the Suzaku source), excluding the compact source MGPS J142714−604755 of unknown origin. The radio fluxes were calculated using MGPS-2, SGPS and PMN data and are reported in Table 1. We found a MGPS-2 source with flux roughly six times higher than that of MGPS J142755−605038, which was considered as the radio counterpart in the modeling of Vorster et al. (2013). We also note that the flux upper limit derived with the PMN data is underestimated as a result of the masking of MGPS J142714−604755, which leads to a large number of pixels in the ON region that are not taken into account in the flux calculation. The source masking is then the main limitation of this method when using moderate angular resolution instruments (the PSF of the PMN is 4.9 , compared to 45 with the MGPS-2). The best-fit power law gives a spectral index of α = −0.38 ± 0.36, consistent with nonthermal emission. We also derived the spectral index in the region encompassing only the brightest part of the emission (located between the two MGPS-2 sources) and we found a softer spectrum α = −0.99 ± 0.34 that may be due to contamination from the compact source MGPS J142714−604755. The radio morphology in the ON region is more similar to that of a center-filled PWN than a shell-type SNR, and the radio spectral index is compatible with that expected from PWN emission. The softness of the X-ray spectrum reported by Fujinaga et al. (2013) indicates a relatively low maximum energy reached by particles, which points toward an ancient PWN scenario. In evolved PWNe, radio and X-ray morphologies can be uncorrelated since they originate from particle populations with different injection and radiative loss timescales. Freshly injected particles produce X-ray emission while the oldest particles have already cooled and diffused away from the pulsar, producing a relic radio PWN. We therefore naturally expect the radio extent of an evolved PWN to be larger than the X-ray extent, which is some- what at odds with that seen in HESS J1427−608. However, it should be noted that radio interferometers cannot reveal large components as extended emission can be partially suppressed during the data reduction; for example, the MOST telescope does not detect structures on angular scales larger than 20-30'. If Suzaku J1427−6051 and the nonthermal radio emission are both associated with HESS J1427−608, the latter is then probably affected by this instrumental limitation. The displacement of the centroid of the radio emission from that of the X-ray emission is found in other evolved PWNe such as Vela X and can be explained by different particle populations. This could explain why the brightest part of the radio emission is displaced from the Suzaku J1427−6051 centroid, even though it is coincident with the northern tip of the extended X-ray emission. It could thus be connected to the compact source MGPS J142714−604755 if the latter is assumed to be the pulsar at the origin of this multiwavelength wind nebula.
The ROSAT/PSPC significance map shows no emission toward HESS J1427−608. The X-ray flux upper limit is therefore not constraining (F X < 1.29 × 10 −8 erg cm −2 s −1 ). This is consistent with the high column density value derived in Fujinaga et al. (2013), which leads to a large absorption of soft X-ray photons along the line of sight. Thus, we derived a loose constraint on the mean magnetic field, that is, B < 263.4 µG at the 5σ confidence level with N H < 1.39 × 10 23 cm −2 (Fujinaga et al. 2013) and Γ TeV < 3.7. Assuming that Suzaku J1427−6051 is the X-ray counterpart of HESS J1427−608, accelerated electrons can be considered in the Thompson regime (Γ X ≈ Γ TeV , with Γ X = 3.1 +0.6 −0.5 and Γ TeV = 2.85 ± 0.22) and B amounts to 10.9 +4.5 −3.2 µG (from Equation 1), which is compatible with that obtained for PWNe and SNRs. The latest XMM-Newton source catalog contains 4XMM J142756.7−605214 located near the center of HESS J1427−608 (black cross in Figure 1, left), which was not detected in the study of Fujinaga et al. (2013). A dedicated X-ray analysis of this region would allow us to determine if 4XMM J142756.7−605214 could be associated with Suzaku J1427−6051.
Three Fermi-LAT sources are located inside the R spec of HESS J1427−608 (3FHL J1427.9−6054, 4FGL J1427.8−6051 and the source detected by Guo et al. 2017, as illustrated in Figure 1, left). Figure 2 (left) shows the corresponding spectral energy distributions (SEDs) reported in the Fermi-LAT catalogs. The spectrum of the point source detected above 3 GeV (Guo et al. 2017) smoothly connects to the H.E.S.S. measurements while the 4FGL catalog reports a source with a pulsar-like spectrum represented by a logarithmic parabola. If the dominant emission above 3 GeV originates from a PWN powered by a putative pulsar whose emission is reported in the 4FGL catalog, we note that the PWN spectrum is very soft in the TeV range and has a relatively low energy cutoff (in the ∼ 500 GeV range), which also points toward an evolved PWN scenario.

Discussion
Since HESS J1427−608 is unresolved at TeV energies, a blazar or binary origin could have been considered. However, given the existence of an extended X-ray emission spatially coincident with the H.E.S.S. source, these two scenarios are clearly disfavored. With an extended and center-filled nonthermal radio emission detected within HESS J1427−608 and a 4FGL Fermi-LAT cataloged source exhibiting a spectrum reminiscent of a pulsar, our multiwavelength data exploitation strengthens the scenario of an evolved PWN.
To constrain the physical parameters in more detail, we modeled the broadband nonthermal spectrum of HESS J1427−608 using the Naima package (Zabalza 2015) assuming a one-zone leptonic model. The electron spectrum is in the form dN where p is the particle spectral index. The minimum and maximum energy of the particles are set to E min = 1 GeV and E max = 1 PeV. Since the distance d is unknown, we arbitrarily set d = 8 kpc (as in Fujinaga et al. 2013). For the IC scattering, we considered the CMB and the infrared and optical photon fields whose temperature and energy density were estimated using the data from Porter et al.  Table 1. Estimated radio fluxes toward HESS J1427−608 and HESS J1458−608 with the associated statistical errors. The flux upper limits are given at the 3σ confidence level. For HESS J1458−608, the values correspond to those derived within the regions R ON1 and R ON2 shown in Figure  3 (right).  (2008) included within the GALPROP 1 project. We obtained T IR = 25 K, U IR = 0.82 eV cm −3 , T opt = 2005 K, and U opt = 1.07 eV cm −3 . We used the radio fluxes reported in Table 1 (excluding the upper limit derived with PMN data, see Section 3.1.2), the X-ray measurements from Suzaku (Fujinaga et al. 2013), the GeV data points from Guo et al. (2017), and the TeV SED reported in the HGPS. To simultaneously fit the radio and X-ray data, a hard particle spectral index such as p = 1.9 is required, and the best fit gives E cut = 12.3 +1.2 −1.2 TeV and B = 10.1 +1.0 −0.7 µG, which is consistent with the magnetic field value inferred in Section 3.1.2. As shown in Figure 2 (right), the low-energy Fermi-LAT data points cannot be reproduced with such a hard particle spectral index. Using p = 2.5 and E min = 25 GeV, the broadband nonthermal emission can be explained with B = 9.9 +0.8 −0.8 µG and E cut = 19.1 +1.7 −1.7 TeV (Figure 2, right). Although the low-energy cutoff in the electron spectrum is questionable, a value of E min = 13 GeV was also required to simultaneously explain the radio and X-ray measurements of the TeV PWN HESS J1356−645 (H.E.S.S. Collaboration et al. 2011). Such values of E min are also in the range of those considered in Kennel & Coroniti (1984) and Ackermann et al. (2011). We note however that significant radio emission could have been suppressed during the data reduction. Finally, it is worth noting that if HESS J1427−608 is an ancient PWN, a more detailed model is likely needed to explain its broadband nonthermal emission, as this is the case for the Vela X PWN whose multiwavelength modeling requires two electron populations (Hinton et al. 2011).
Deeper radio observations and dedicated TeV analyses would be helpful to obtain better insight into the source morphology. Pulsation searches on 4FGL J1427.8−6051 are also necessary to possibly assess the nature of this Fermi-LAT source as a pulsar and dedicated X-ray observations and analyses are required to search for this putative compact source. If a pulsar is detected with a spin-down luminosity high enough to power a detectable nebula despite its likely large distance, HESS J1427−608 could thus be firmly identified as a TeV PWN.

Source presentation
In the HGPS, HESS J1458−608 is described by a twodimensional symmetric Gaussian with σ = 0.37 • ± 0.03 • and has a significance of √ TS = 11.5 (H. E. S. S. Collaboration et al. 2018b). The TeV spectrum is hard with Γ TeV = 1.81 ± 0.14 and F 1−10 TeV = (5.28 ± 1.30) × 10 −12 erg cm −2 s −1 obtained with R spec = 0.5 • . At the center of HESS J1458−608 lies a Fermi-LAT γ-ray pulsar PSR J1459−6053 with a spin-down power oḟ E = 9.1 × 10 35 erg s −1 (Ray et al. 2011) and a characteristic age of τ c = 64 kyr (Marelli et al. 2011). Data taken with Swift/XRT (6 ks) revealed a point source, offset from the γ-ray pulsar by ∼ 9.8" and this source was suggested as the X-ray counterpart of PSR J1459−6053 (Ray et al. 2011). XMM-Newton observations of 50 ks confirmed the nonthermal X-ray emission from the pulsar with a derived column density of N H = 6.4 +8.9 −5.1 × 10 21 cm −2 (Pancrazi et al. 2012). Since the pulsar is not detected in the radio band, the distance of the system is unknown. No associated X-ray and γ-ray emissions originating from a PWN have been detected so far. Figure 3 (left) shows the H.E.S.S. significance map, indicating that the TeV morphology may be more complex than a symmetric Gaussian and that two Fermi-LAT cataloged sources lie in the western part of HESS J1458−608.

Multiwavelength constraints
We extracted archival radio continuum data that show emission in the vicinity of the pulsar and in the western part of HESS J1458−608, the latter being spatially coincident with the Fermi-LAT sources 3FGL J1456.7−6046 and 4FGL J1456.7−6050c (Figure 3, right). We defined two regions that encompass the diffuse emissions within the H.E.S.S. source: one in the vicinity of the pulsar (called R ON1 ) and one toward the western part of the H.E.S.S. source (called R ON2 ). The estimated radio fluxes using MGPS-2, Parkes, and PMN are given in Table 1 whose best-fit power law gives α 1 = −0.44 ± 0.18 (for R ON1 ) and α 2 = −0.28 ± 0.20 (for R ON2 ). Both spectral indexes indicate nonthermal radio emission and it is not clear whether these two regions originate from the same system. The spectral indexes of R ON1 and R ON2 are compatible with those expected from both PWNe and shell-type SNRs. The spectral index of R ON1 is even closer to that of a SNR than a PWN, but the absence of a shell-like morphology and the location of the radio emission close to the energetic PSR J1459−6053 makes a PWN scenario more likely.
The pulsar is seen in Chandra, Swift/XRT, XMM-Newton, and Suzaku data but no extended emission appears in these images. The ROSAT/PSPC data show no significant emission and give a tight constraint on the X-ray flux F X < 2.13 × 10 −11 erg cm −2 s −1 (with N H < 1.5 × 10 22 cm −2 , Pancrazi et al. 2012, and Γ TeV < 1.95), giving B < 10.1 µG at the 5σ confidence level.
The pulsar PSR J1459−6053 is reported in the 2PC, 3FGL, 3FHL, and 4FGL catalogs. The Fermi-LAT catalogs also contain two sources located in the western part of HESS J1458−608 (3FGL J1456.7−6046 and 4FGL J1456.7−6050c, the latter being labeled as confused). Figure 4 (left) shows that the best-fit spectrum of 4FGL J1456.7−6050c in the whole energy range is hard (Γ GeV = 1.30 ± 0.28) and connected to the spectrum of HESS J1458−608. The SED at low energy of 4FGL J1456.7−6050c is likely contaminated by the nearby γ-ray pulsar, as might be the case for the SED of 3FGL J1456.7−6046. However, we keep in mind that the detec-tion of 4FGL J1456.7−6050c may be spurious since the source is labeled as confused.

Discussion
HESS J1458−608 is a puzzling source, with a complex and elongated TeV morphology, and shelters the energetic X-ray and γ-ray pulsar PSR J1459−6053. With nonthermal radio emission located in the western part of the H.E.S.S. source, which is spatially coincident with a GeV source exhibiting a hard spectrum (4FGL J1456.7−6050c), part of the emission from HESS J1458−608 could originate from a PWN. We also reported nonthermal radio emission in the vicinity of PSR J1459−6053. These extended nonthermal radio emissions, possibly connected to each other, are spatially coincident with the western part of the elongated shape of HESS J1458−608, as seen in Figure 3.
As for HESS J1427−608, we modeled the broadband nonthermal spectrum of HESS J1458−608 assuming a one-zone leptonic scenario and a distance of d = 8 kpc. The electron spectrum is defined as in Section 3.1.3 and the temperature and energy density of the infrared and optical photon fields are T IR = 25 K, U IR = 0.88 eV cm −3 , T opt = 2005 K and U opt = 1.94 eV cm −3 (estimated using GALPROP code). We used the radio fluxes (from both regions) listed in Table 1, the upper limit on the X-ray flux derived in Section 3.2.2 and the TeV SED from the HGPS. Owing to the lack of X-ray data points and the absence of a VHE cutoff, B and E cut cannot be accurately derived. As shown in Figure 4 (right), the data can be reproduced with E min = 1 GeV, E cut = 60 TeV, B = 10 µG, p = 2 and W e = 7 × 10 47 ×(d/8 kpc) 2 erg. For illustration, the best-fit spectrum of 4FGL J1456.7−6050c is represented in Figure 4 (right), keeping in mind that it was derived for a point-source morphology and that the source detection needs to be confirmed.
Deeper X-ray observations and dedicated TeV analyses are required to provide better insight into the source morphology and to better constrain the broadband nonthermal spectrum. A study at HE would also be useful to confirm the detection of a source (and to unveil its morphology) exhibiting a hard  spectrum connected to that of HESS J1458−608.
3.3. HESS J1427−608 and HESS J1458−608 within the context of the population of TeV PWNe We found extended nonthermal radio emission toward HESS J1427−608 and HESS J1458−608. For both sources, the faint radio emission, although not clearly correlated with the TeV morphology, points toward an ancient PWN. The mean magnetic field (inferred from Equation 1 and in more detail through broadband modeling) is consistent with such a scenario and we found a pulsar-like and a PWN-like spectrum in the vicinity of HESS J1427−608 and HESS J1458−608. We note however that HESS J1458−608 could be a composition of multiple sources, since the hard spectrum seen by the Fermi-LAT (although this detection needs to be confirmed) originates from the western edge of the TeV source with a nonthermal radio counterpart. Therefore, it is not clear whether this part is connected to the main emission region of HESS J1458−608, where the pulsar PSR J1459−6053 is located. Assuming a PWN scenario for HESS J1427−608 and HESS J1458−608, we want to estimate the pulsar spin-down power and its characteristic age to place both sources within the context of the population of TeV PWNe. The H.E.S.S. collaboration undertook a PWN population study (using a total of 37 sources including firmly identified PWNe and PWN candidates, H. E. S. S. Collaboration et al. 2018a) and derived a relation between the TeV surface brightness S and the pulsar spin-down powerĖ as follows: where S andĖ are in units of erg pc −2 s −1 and 10 36 erg s −1 , respectively. The surface brightness can be written as where L 1−10 TeV , F 1−10 TeV , R PWN and σ are, respectively, the luminosity (erg s −1 ) and the integral energy flux (erg pc −2 s −1 ) between 1 and 10 TeV, the radius of the PWN (pc) and the Gaussian σ extent (radians). Since no pulsar was detected in the vicinity of HESS J1427−608 and the distance of PSR J1459−6053 is unknown, HESS J1427−608 and HESS J1458−608 were not considered in the H.E.S.S. PWN population study. For HESS J1427−608, we have F 1−10 TeV = (1.43 ± 0.34) × 10 −12 erg cm −2 s −1 and σ < 0.067 • , which givesĖ > 1.38 × 10 36 erg s −1 . Assuming that Suzaku J1427−6051 is the X-ray counterpart of HESS J1427−608, the TeV extent should be larger than the X-ray extent (σ X = 0.9') as it originates from older electrons that have diffused farther away from the pulsar than those emitting X-ray synchrotron. With σ > 0.9', we obtainedĖ < 6.91 × 10 38 erg s −1 . For HESS J1458−608, with F 1−10 TeV = (5.28 ± 1.30) × 10 −12 erg cm −2 s −1 and σ = 0.37 • ± 0.03 • , we obtainedĖ = 2.1 +3.4 −1.6 × 10 35 erg s −1 , which is only 2σ away from the measurement ofĖ = 9.1 × 10 35 erg s −1 for PSR J1459−6053 (Ray et al. 2011). If PSR J1459−6053 and HESS J1458−608 are associated, this difference could indicate that either the integral energy flux of the TeV source was underestimated or its Gaussian σ extent was overestimated. Given the complex TeV morphology of HESS J1458−608 (Figure 3, left), possible source confusion could lead to a σ extent larger than that associated with the PWN produced by PSR J1459−6053.
The characteristic age of a pulsar τ c is expressed as where n is the braking index, t age the age of the system, and τ 0 the initial characteristic age of the pulsar (Gaensler & Slane 2006). Assuming n = 3 (appropriated for magnetic dipoles), we have τ c ≈ t age for evolved systems (t age >> τ 0 ). Since these relic PWNe no longer inject particles into the system, the maximum energy reached by particles is limited by radiative losses. The break energy E b (above which electrons significantly suffer from synchrotron losses) is found equating t age = τ sync , where τ sync is the synchrotron loss time where E TeV and B 100µG are the particle energy and mean magnetic field in units of TeV and 100 µG, respectively (Parizot et al. 2006). In evolved systems limited by radiative losses, we can therefore assume E b = E cut . Our estimate on τ c relies on the assumption that HESS J1427−608 and HESS J1458−608 are evolved PWNe and that we assume τ c = t age . For HESS J1427−608, we used B = 10.0 +1.1 −0.9 µG and E cut = 15.7 +5.1 10 0 10 1 10 2 10 3 Characteristic age c (kyr) TeV (combining the constraints obtained with E min = 1 and 25 GeV in Section 3.1.3) and we deduced an age τ c = t age = 8.0 [4.9, 13.6] kyr (Equation 5). For HESS J1458−608, we derived a lower limit of E cut > 30 TeV since no break or cutoff is visible in Figure 4 (right). In order to bound the cutoff energy and the mean magnetic field, we assumed E cut < 100 TeV and B > 3 µG, which are appropriate for evolved PWNe. With B < 10.1 µG (derived in Section 3.2.2), we found 1.2 < τ c = t age < 46.3 kyr. Figure 5 depicts the spin-down power and the characteristic age of pulsars with either a firmly identified TeV PWN or a PWN candidate (reproduced from H. E. S. S. Collaboration et al. 2018a, with the addition of HESS J1427−608 and HESS J1458−608). The estimated values for HESS J1427−608 and HESS J1458−608 are contained within those of the sample of identified PWNe and PWN candidates. However, HESS J1427−608 and HESS J1458−608 appear to be similar to a relatively young PWN (with a highĖ) and an evolved PWN (with anĖ estimate lying among the lowest values of the TeV PWN population), respectively, whereas their γ-ray spectrum (shown in Figure 2 and 4) may suggest the opposite. It is nevertheless possible that the extent of HESS J1458−608 was overestimated as a result of source confusion and that the hard TeV spectral shape originates from multiple source contribution. Owing to the inherent uncertainties in such approach, we cannot draw precise conclusions except that HESS J1427−608 and HESS J1458−608 seem to broadly follow the general trend of the population of TeV PWNe.
Finally, we note that the characteristic age of HESS J1427−608 is comparable to that of the evolved PWN HESS J1356−645 (τ c = 7.3 kyr). Assuming that the intrinsic size of HESS J1427−608 is similar to that of HESS J1356−645 (R ≈ 8.4d 2.4 pc, where d 2.4 is the distance in units of 2.4 kpc), the distance to HESS J1427−608 would be d ≈ 10 kpc using an angular TeV size of σ = 0.048 • (or d > 7 kpc for σ < 0.067 • ). HESS J1458−608 would be more comparable to the PWN HESS J1825−137 (τ c = 21 kyr, R ≈ 50d 4 pc). Assuming that their intrinsic sizes are similar, that would place HESS J1458−608 at a distance d ≈ 8 kpc for an angular TeV size of σ = 0.37 • . The distances to HESS J1427−608 and HESS J1458−608 could therefore be compatible with the relatively large column densities derived for their associated X-ray sources (Suzaku J1427−6051 and PSR J1459−6053, see Sections 3. 1.1 and 3.2.1). Regardless of the dominant particle transport mechanism (diffusion or advection, as discussed for

New VHE γ-ray source candidates
Although very robust, the HGPS detection pipeline, relying on an iterative template fitting in which the γ-ray sources are treated as two-dimensional symmetric Gaussians, faced the difficulty of revealing multiple components in regions with large source confusion. The visual inspection of the HGPS images led us to focus on two possible γ-ray excesses lying in such complex regions, which eluded detection by the algorithm. The first is located south of the unidentified source HESS J1632−478 and the other lies at the position of the synchrotron-emitting SNR G28.6−0.1, at the edge of HESS J1843−033. Recently, Remy et al. (2020) developed an alternative algorithm based on image processing and pattern recognition techniques to detect γ-ray sources without strong prior morphological assumptions: structural information is extracted using an edge detection operator and the detected objects are found as local maxima after applying the Hough circle transform. This algorithm, applied on the HGPS maps, has provided a list of objects that warrant further investigation. Each of the two above-mentioned excesses, which we independently identified after examining the HGPS maps, has a counterpart in the catalog of Remy et al. (2020) as a detected object uncataloged in the HGPS (labeled HC_147 and HC_382, respectively). Although this strengthens the interest of these two excesses, only a dedicated H.E.S.S. data analysis could confirm these as proper VHE sources. In the following, we consider them as such and investigate their origin through the existing multiwavelength data.

The south of HESS J1632−478
The unidentified source HESS J1632−478 is described in the HGPS by a Gaussian component with σ = 0.18 • ± 0.02 • . The HGPS significance map is shown in Figure 6 (left) and shows a γ-ray excess with a peak significance at about 8σ in the south of HESS J1632−478 (outside of the HGPS R spec ). The pulsar PSR J1633−4805 (Ė = 8.5 × 10 33 erg s −1 and d = 6.2 kpc, Manchester et al. 2005) and the uncertain SNR G336.7−0.3 (Green 2017) are located near this γ-ray excess.
We explored archival radio continuum data toward the south of HESS J1632−478. The MGPS-2 map is given in Figure 6 (right) and shows a bright emission with some distorted shelllike structures. Several H ii regions reported in the WISE catalog (Anderson et al. 2014) and one MC (Rice et al. 2016) are spatially coincident with the γ-ray excess. Using data from the MGPS-2, SGPS and Parkes, we derived a radio spectral index of α = 0.64 ± 0.04 in the ON region encompassing the radio emission (white circle in Figure 6, right). The spectral index indicates Article number, page 9 of 19 A&A proofs: manuscript no. TeVsources that the radio emission is thermal in nature and could originate from the H ii regions. Since these H ii regions do not constitute the entire radio emission, it is possible that there is a subdominant nonthermal emission.
We also found a X-ray point-source detected by ROSAT/PSPC (2RXS J163352.2−480643) located in the western part of the ON region, which is also seen by ASCA/GIS although not cataloged. Four XMM-Newton extended sources are also reported in this region (see The Fermi-LAT extended source FGES J1633.0−4746 (described by a disk with a radius r = 0.61 • ; Ackermann et al. 2017) encompasses the TeV γ-ray excess but its large extent points toward possible source confusion. The Fermi-LAT residual count map above 10 GeV shows emission that is spatially coincident with the TeV γ-ray excess, but a dedicated analysis would be required to potentially reveal an extended component.
To conclude, we found radio emission that is spatially coincident with the significant γ-ray excess in the south of HESS J1632−478, exhibiting a complex morphology reminiscent of a distorted shell-like structure with a thermal spectrum. Several H ii regions located within the γ-ray excess indicate that the TeV emission could arise from a star-forming region, which can accelerate particles up to γ-ray emitting energies (with colliding stellar winds). Such regions were recently proposed as sources of Galactic CRs . Extended Xray sources are found in the western part of this region with un-known origin. Deeper X-ray observations would help understand the nonthermal processes occurring in this region. Detailed HE and VHE γ-ray data analyses would be of great interest to investigate the morphology and spectrum of this γ-ray excess and potentially reveal a new VHE source. The unidentified source HESS J1843−033 is defined as a large extended Gaussian component (σ = 0.24 • ± 0.06 • ) and results from the merging of two Gaussian components, which were previously detected by the HGPS pipeline, as illustrated in the HGPS significance map (Figure 7, left). Several SNR candidates revealed in the THOR radio survey (Anderson et al. 2017) overlap with HESS J1843−033 (Figure 7, right). The radio SNR G28.6−0.1 and the energetic pulsar PSR J1844−0346 (Ė = 4.2 × 10 36 erg s −1 and τ c = 11.6 kyr with unknown distance; Clark et al. 2017) are located within one of the previously detected Gaussian components. Figure 7 highlights the complexity of the region and shows a peak significance at nearly 6σ at the position of the synchrotron-emitting SNR G28.6−0.1.
The eastern and southern parts of the shell of G28.6−0.1 were first detected in radio, with a flux of ≈ 2 Jy at 1.4 GHz (Helfand et al. 1989). Bamba et al. (2001) revealed nonthermal emission from an extended ASCA source spatially coincident with G28.6−0.1, with a power-law spectral index Γ X = 2.1 +0.3 −0.4 , a column density of N H = 2.6 +0.8 −0.6 × 10 22 cm −2 and a distance estimate between 6 kpc and 8 kpc. Chandra observations of G28.6−0.1 also revealed thermal and nonthermal X-ray emission components with a faint shell-like morphology (Ueno et al. 2003). The association between G28.6−0.1 and the nearby X-ray pulsar PSR J1844−0346 was considered owing to the column density of N H = 5 +4 −3 × 10 22 cm −2 , similar to that obtained for the SNR (Zyuzin et al. 2018). The transverse velocity of the pulsar can be written as v t = 4.74µD km s −1 , where µ is the proper velocity of the pulsar (mas/yr) and D the distance (kpc) (Lyne & Lorimer 1995). Since µ has not been measured, we can approximate µ ≈ θ τ c , where θ is the angular distance between the current position of the pulsar and its place of birth and τ c its characteristic age. Assuming a distance of 4.3 kpc (estimated from the empirical relation obtained for γ-ray pulsars; Saz Parkinson et al. 2010) and 6 kpc (the closest estimate for G28.6−0.1), the transverse velocity of the pulsar is v t ≈ 1020 km s −1 and ≈ 1422 km s −1 (Zyuzin et al. 2018), while typical pulsar velocities range between 400−500 km s −1 (Lyne & Lorimer 1995). The association between PSR J1844−0346 and G28.6−0.1 is therefore unlikely.
Using data from the TGSS, NVSS, and PMN surveys, we investigated the radio spectrum of G28.6−0.1 and we derived a spectral index of α = −0.44 ± 0.01, fully consistent with that usually measured for shell-type SNRs. The brightest parts of the shell are reported as sources in the TGSS catalog (at 150 MHz). By summing up their respective fluxes we obtained ∼ 8 Jy, which is broadly consistent with our estimate of 9.01 ± 0.17 Jy obtained for the entire SNR. At X-ray energies, the inset in Figure 7 (right) depicts the Chandra image, which shows a faint shell-like morphology, previously revealed by Ueno et al. (2003). No Fermi-LAT cataloged sources are reported close to the SNR. We also explored the archival multiwavelength data toward the pulsar PSR J1844−0346 and we did not find any radio or X-ray counterpart that could indicate a possible PWN. Radio data from MAGPIS and infrared data from Spitzer show a bright extended emission southeast of PSR J1844−0346. The Simbad database indicates the presence of the star-forming region N49 containing H ii regions, whose distance is estimated to be 5.1 kpc or 9.7 kpc; the former is the most likely distance (see Anderson & Bania 2009;Dirienzo et al. 2012). With an angular distance between N49 and PSR J1844−0346 of θ = 181.5", the proper velocity of the pulsar is µ = 15.65 t −1 11.6 mas yr −1 , assuming t ≈ τ c = 11.6 kyr. Taking the near distance to N49 (d = 5.1 kpc), the transverse velocity would be v t ∼ 378 d 5.1 t −1 11.6 km s −1 , making the association of PSR J1844−0346 with N49 more likely than with G28.6−0.1.
We shed light on a significant TeV excess that is spatially coincident with the synchrotron-emitting shell-type SNR G28.6−0.1. We investigated its radio spectrum and we derived a spectral index of α = −0.44 ± 0.01. Multiwavelength data do not show any emission from a potential PWN originating from the nearby and energetic pulsar PSR J1844−0346. While the association between PSR J1844−0346 and G28.6−0.1 was previously considered, we found that the pulsar more likely originates from the star-forming region N49. The radio and X-ray synchrotronemitting SNR G28.6−0.1 is a promising source to be studied at HE and VHE since it should normally be able to accelerate particles up to γ-ray emitting energies in a leptonic scenario. More generally, the unidentified source HESS J1843−033 is an interesting source to investigate at VHE, given the numerous overlapping SNR candidates. Among these, we note that there is another excess (at ∼4−5σ in the HGPS maps) lying on the eastern part of the shell of G28.78−0.44. With unprecedented angular resolution and sensitivity, the next generation of Cherenkov telescopes CTA (Cherenkov Telescope Array Consortium et al. 2017) will give better insight into this region and might be able to confirm the VHE γ-ray emission from G28.6−0.1 (and potentially from G28.78−0.44 as well), making it a new TeV-emitting shell-type SNR.

Conclusions
Given the large number of HGPS sources that are not firmly identified, we developed a generic code aiming to constrain the origin of their TeV emission through the exploitation of multiwavelength data. The algorithm is based on an automatically retrieval of archival data from radio continuum, X-ray, infrared, and GeV instruments toward any region of the sky to search for faint counterparts through a careful visual inspection in every available image. The pipeline was completed by an association procedure using the catalogs of known objects (SNRs, PWNe, H ii regions, etc.) and those related to the instruments whose data are exploited. The main constraints on the origin of the TeV emission are obtained through the derivation of a radio spectral index that helps us disentangle thermal from nonthermal emission and through the estimate of the mean magnetic field under the assumption of a leptonic scenario. Finally, the Fermi-LAT cataloged source spectra are used to search for a smooth connection between GeV and TeV energies. This algorithm is well suited for isolated sources but faces some limitations for large and complex regions, although it gives an overall insight into the possible multiple components as the origin of the TeV emission.
We applied this pipeline on two unidentified sources reported in the HGPS: HESS J1427−608 and HESS J1458−608. We found extended nonthermal radio emission toward both sources and a mean magnetic field value consistent with that expected for evolved PWNe. Fermi-LAT data revealed a pulsar-like spectrum and a PWN-like spectrum in the vicinity of HESS J1427−608 and HESS J1458−608, respectively. We modeled the broadband nonthermal emission of both sources in a leptonic scenario with B 10 µG, which is consistent with that obtained from ancient PWNe. We estimated the spin-down power and the characteristic age of the putative pulsars and we found that these are broadly in line with those of the population of known TeV PWNe. Deeper X-ray observations are necessary to potentially reveal a compact source in the vicinity of HESS J1427−608. A detailed search for pulsation on 4FGL J1427.8−6051, which exhibits a pulsarlike spectrum, could help reveal its nature. Follow-up X-ray observations could also reveal an extended source associated with HESS J1458−608, and dedicated HE and VHE γ-ray data analyses are required to reveal in details the morphology of this unidentified VHE source.
We also shed light on a possibly significant, yet uncataloged, γ-ray excess in the HGPS data, located in the south of the unidentified source HESS J1632−478, for which we found an extended radio counterpart with a complicated shell-like structure. The radio emission was found to be thermal and could be associated with a star-forming region, which could accelerate particles up to γ-ray emitting energies. We also highlighted another uncataloged γ-ray excess, that is spatially coincident with the synchrotron-emitting shell-type SNR G28.6−0.1, lying near the unidentified source HESS J1833−043. Both sources have multiwavelength counterparts and could be revealed as new VHE sources in the future.
This code was developed with the purpose of applying a uniform method on the unidentified TeV sources to constrain the origin of their emission. Given the constraints obtained with archival data, this work also aims to advocate deeper multiwavelength observations, searches for pulsation, and dedicated analyses to be able to firmly identify TeV sources. Remy et al. (2020) recently investigated blind search methods for VHE γ-ray source detection based on image processing and pattern recognition techniques and provided a list of promising HGPS TeV source candidates. Future work will include the application of this pipeline on these source candidates. If interesting multiwavelength counterparts are to be found, detailed VHE data reanalyses would then be warranted to firmly confirm these candidates as new VHE sources. Our next improvements of the code will include the use of H i and CO data to assess the environment of the sources and the analysis of X-ray observations to derive better constraints on the synchrotron spectrum. The convolution of the radio and X-ray maps with the H.E.S.S. PSF would be also needed to quantify the spatial correlations at different wavelengths.
The next generation of Cherenkov telescopes CTA (Cherenkov Telescope Array Consortium et al. 2017) will certainly detect a significant number of new sources, leading to a higher degree of source confusion and hence, a larger number of unidentified sources than what is currently observed. Dubus et al. (2013) estimated that 20−70 shell-type SNRs and 300−600 PWNe will be detected during the CTA Galactic Plane Survey. With significantly improved sensitivity and angular resolution compared to current Cherenkov instruments, CTA will provide new insights into the Galactic γ-ray sky and the synergy with multiwavelength observations will be necessary to constrain the origin of the observed TeV emission. The large field-of-view Xray instrument eROSITA 1 , recently launched, will soon give an unprecedented high-energy map up to 10 keV (Merloni et al. 2012), allowing us to search for faint X-ray counterparts along the whole Galactic plane. Data from the Square Kilometre Array (Weltman et al. 2018) are also very promising to help identify γ-ray sources. Thus, future VHE observations and multiwavelength data exploitation supported by new generation instruments will probably reveal the nature of a significant number of unidentified TeV sources in order to assess their importance within the issue of the origin of the Galactic CRs.  Table A.1. Properties of the radio surveys. The type "S" and "I" denotes a "single-dish telescope" or "interferometers", respectively. The pointsource sensitivity is given, as well as the FWHM. If the FWHM is asymmetric, the major axis is given. For the TGSS and the CGPS, the FWHM depends on the declination (indicated by the asterisk). The links to retrieve these data are labeled with numbers and given below. Numbers followed by an asterisk indicate that a source catalog, associated with the instrument, is available and used.
Archival data are retrieved at the following links:  Table A.2 summarizes the description of the X-ray instruments whose archival data are retrieved through the pipeline. We automatically extract all the observations whose center is comprised within an angular distance from the source of interest. For ROSAT/PSPC and XMM-Newton, we create background-subtracted and exposure-corrected maps, while the count maps of Swift/XRT, ASCA and Suzaku are only corrected from the exposure because the corresponding background maps are not available. Except for large field-of-view instruments (such as ROSAT/PSPC, Integral/IBIS-ISGRI, and Swift/BAT), we create a mosaic of these images in case of multiple observations. For Chandra data, we directly retrieve stacked and subtracted-background images (combining multiple observations), which are also corrected from the exposure and the vignetting. For Integral/IBIS-ISGRI and Swift/BAT, we extract a subregion of the flux, error flux, and significance maps.  6, 4.5, 5.8, 8, 21, 24, 870, and 1100 µm (available on the MAGPIS website https://third.ucllnl.org/gps/). Using ten years of Fermi-LAT data, we also created a binned count map with energies between 10 and 500 GeV, setting a pixel size of 0.05 • and 10 energy bins (version 1.0.10 of the Fermitools). The map encompasses the HGPS with latitude |b| ≤ 5 • . The SOURCE event class was selected, which is a compromise between background rejection and statistics and we imposed a maximum zenith angle on the photon arrival direction of 90 • to reduce the contamination of the Earth limb. We excluded time intervals during which the satellite passed through the South Atlantic Anomaly and when the rocking angle was more than 52 • . We modeled the contribution of the Galactic and isotropic diffuse emissions using the files gll_iem_v07.fits and iso_P8R3_SOURCE_V2.txt 1 (with the normalizations fixed to 1 and the spectral index to 0) using the instrument response functions P8R3_SOURCE_V2, and we created a residual count map between 10 GeV and 500 GeV.
J. Devin et al.: Multiwavelength constraints on unidentified TeV sources flux due to a too small number of remaining unmasked pixels. The source masking is then the main limitation of this method 1 . Figure B.3 depicts the distribution of the unmasked pixels in the OFF region (used for the background estimation) of the MGPS-2 data (left), and the SED with the best-fit power law on the calculated fluxes. We found a radio spectral index of α = −0.03 ± 0.06 that is compatible with the value obtained in H.E.S.S. Collaboration et al. (2011). The difference can be explained by a slightly different background estimation method and a small difference in the size and position adopted for the ON region.

Appendix C: Mean magnetic field estimation
We use ROSAT/PSPC data in the 0.9-2.4 keV band to estimate the flux (or the flux upper limit) and to constrain the value of the mean magnetic field, given the measured flux at TeV energies. We take advantage of this large field-of-view all-sky instrument, which covers the entire HGPS region. We use the 1RXS catalog reporting the brightest sources detected with ROSAT/PSPC (Voges et al. 1999) and the 2RXS catalog (Boller et al. 2016) to mask the sources located outside the ON region. We use the ring background and the reflected background methods, and we calculate the significance following the prescription in Li & Ma (1983), using a correlation radius equal to the flux extraction radius defined in the HGPS (R spec ). For the ring background method, the OFF region is defined between R in = R spec + 1 pixel and R out = (10R 2 spec + R 2 in ) • so that the OFF to ON area ratio amounts to ∼ 10 for large sources. For the reflected background method, we use four circular regions (with a radius R spec ) located at an angular distance of 3 × R spec from the source center. If the X-ray emission is not significant, an upper limit on the flux is calculated at a 5σ confidence level. The X-ray flux is simulated using the HEASARC tools WebPimms 1 and defined as where F 0 is the non-absorbed flux, σ(E) the photoelectric absorption cross section (in units of cm 2 ) and N H the column density (in units of cm −2 ). For each simulated flux, a combination of (Γ, N H ) corresponds to a count rate between 0.9 and 2.4 keV that we simulate in the image with Monte-Carlo sample. The simulation stops when the significance seen in the ROSAT/PSPC data is reached (or when the 5σ-threshold is reached in case of a non-detection) and gives the X-ray flux (or upper limit on flux) in the parameter space (Γ, N H ). Then, we estimate the mean magnetic field assuming a one-zone model and that electrons radiating synchrotron and IC emission (scattering on the cosmic microwave background) are in the Thompson regime (Γ sync = Γ IC = Γ, where Γ is the photon spectral index). If we describe the electron spectrum with a power law (Kγ −p , γ and p being the Lorentz factor and the electron spectral index), the ratio of the synchrotron to IC fluxes is given by where U B ∼ 2.5B 2 −5 eV cm −3 and U CMB = 0.26 eV cm −3 are the magnetic field and cosmic microwave background energy densities respectively (with B −5 = B 10 µG ). The synchrotron and IC radiations occur at the characteristic energies of γ sync 1.4 × 10 8 E 1/2 sync,keV B −1/2 −5 and γ IC 3.6 × 10 7 E 1/2 IC,TeV (C.3) Using the above equations, the mean magnetic field can be expressed as where F is the flux between E 1 and E 2 , G(Γ) (0.1 × 15 Γ−2 ) 1/Γ and Γ = (p + 1)/2. If Γ = 2, the mean magnetic field is simply expressed as F sync F IC 10B 2 −5 (C.5) We can thus estimate the mean magnetic field using the X-ray flux derived with ROSAT/PSPC data and the TeV flux and spectral index reported in the HGPS. Below the method is applied on the SNR RX J1713.7−3946 for illustration.
The spectral extraction radius of the SNR RX J1713.7−3946 is R spec = 0.6 • and the spectral parameters reported in the HGPS are F 1 TeV = (2.13 ± 0.94) × 10 −11 cm −2 s −1 TeV −1 and Γ TeV = 2.32 ± 0.02. We masked in the ROSAT/PSPC images all the 1RXS and 2RXS cataloged sources whose centroid are located outside R spec . Figure C.1 (left) shows the significance map using the reflected background method, which is similar to that obtained using the ring background method. As expected, the significance is well above the 5σ threshold. The mean value of the magnetic field is given in Figure C.1 (right) with respect to the TeV photon spectral index and the column density. Following Sano et al. (2015), who estimated the variation of the column density over the entire SNR, we took N H = (0.3 − 1.0) × 10 22 cm −2 . Within these TeV spectral index and column density ranges, we determined a mean magnetic field of B = [10.48 − 16.37] µG, which is consistent with B = 14.26 ± 0.16 µG, obtained when modeling the broadband nonthermal emission of the SNR with a leptonic model (H. E. S. S. Collaboration et al. 2018b). (Right) Mean magnetic field as a function of the TeV photon spectral index and the column density. The full and dashed horizontal lines are the TeV spectral index and its associated statistical errors, as reported in the HGPS. The vertical dashed lines correspond to the column density range derived from Sano et al. (2015).