# Constraining the Origin of the Puzzling Source HESS J1640–465 and the PeVatron Candidate HESS J1641–463 Using Fermi-Large Area Telescope Observations

A. Mares, M. Lemoine-Goumard, F. Acero, C.J. Clark, J. Devin, S. Gabici,

J.D. Gelfand, D.A. Green, M.-H Grondin

► To cite this version:

A. Mares, M. Lemoine-Goumard, F. Acero, C.J. Clark, J. Devin, et al. Constraining the Origin of the Puzzling Source HESS J1640–465 and the PeVatron Candidate HESS J1641–463 Using Fermi-Large Area Telescope Observations. The Astrophysical Journal, 2021, 912 (2), pp.158. 10.3847/1538-4357/abef62. hal-03203620

HAL Id: hal-03203620

https://hal.science/hal-03203620v1

Submitted on 6 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Constraining the Origin of the Puzzling Source HESS J1640–465 and the PeVatron Candidate HESS J1641–463 Using Fermi-Large Area Telescope Observations

A. Mares<fnref n="1" />, M. Lemoine-Goumard<fnref n="1" />, F. Acero<fnref n="2" />, C. J. Clark<fnref n="3" />, J. Devin<fnref n="1" />, S. Gabici<fnref n="2" />, J. D. Gelfand<fnref n="4" />, D. A. Green<fnref n="5" />, and M.-H. Grondin<fnref n="1" />

Received 2020 December 16; revised 2021 March 15; accepted 2021 March 15; published 2021 May 17

## Abstract

There are only a few very-high-energy sources in our Galaxy that might accelerate particles up to the knee of the cosmic-ray spectrum. To understand the mechanisms of particle acceleration in these PeVatron candidates, Fermi-Large Area Telescope (LAT) and High-Energy Stereoscopic System (H.E.S.S.) observations are essential to characterize their γ-ray emission. HESS J1640–465 and the PeVatron candidate HESS J1641–463 are two neighboring (0.25°) γ-ray sources, spatially coincident with the radio supernova remnants (SNRs) G338.3–0.0 and G338.5+0.1. Detected both by H.E.S.S. and the Fermi-LAT, we present here a morphological and spectral analysis of these two sources using 8 yr of Fermi-LAT data between 200 MeV and 1 TeV with multiwavelength observations to assess their nature. The morphology of HESS J1640–465 is described by a 2D Gaussian ($σ = 0.053° ± 0.011°_{stat} ± 0.03°_{syst}$) and its spectrum is modeled by a power law with a spectral index $Γ = 1.8 ± 0.1_{stat} ± 0.2_{syst}$. HESS J1641–463 is detected as a point-like source and its GeV emission is described by a logarithmic-parabola spectrum with $α = 2.7 ± 0.1_{stat} ± 0.2_{syst}$ and significant curvature of $β = 0.11 ± 0.03_{stat} ± 0.05_{syst}$. Radio and X-ray flux upper limits were derived. We investigated scenarios to explain their emission, namely, the emission from accelerated particles within the SNRs spatially coincident with each source, molecular clouds illuminated by cosmic rays from the close-by SNRs, and a pulsar/pulsar wind nebula origin. Our new Fermi-LAT results and the radio and flux X-ray upper limits pose severe constraints on some of these models.

*Unified Astronomy Thesaurus concepts: Pulsars* (1306); *Gamma-ray astronomy* (628); *Supernova remnants* (1667); *Cosmic rays* (329)

## 1. Introduction

HESS J1640–465 is a very-high-energy (VHE) γ-ray source located in the Galaxy discovered by the High-Energy Stereoscopic System (H.E.S.S.) during its Galactic Plane Survey (Aharonian et al. 2006). First detected as a point source, further observations revealed a Gaussian extension of $σ = 4.32 ± 0.18$, with a spectrum described by an exponential cutoff power law ($Γ = 2.11 ± 0.09$; Abramowski et al. 2014a) with a cutoff energy of $E_c = 6.0_{1.2}^{2.0}$ TeV. HESS J1640–465 coincides with a dense molecular cloud (G338.4+0.1) overlapping the northwestern part of the supernova remnant (SNR) shell of G338.3–0.0 (Abramowski et al. 2014a). The age of this radio shell-type SNR ($σ = 7'$; Whiteoak & Green 1996) remains uncertain due to the large amount of surrounding interstellar material that could impact the evolution of the SNR, and is estimated to be between 1 and 8 kyr. X-ray observations of the region with XMM-Newton (Funk et al. 2007) revealed an extended nonthermal source ($σ = 27'' ± 3''$) close to the center of the remnant and coincident with the TeV emission. Further observations with Chandra (Lemiere et al. 2009) found a compact X-ray object that was suggested to be the associated pulsar. No X-ray emission was detected from the SNR G338.3–0.0 because of the low temperature of the SNR ejecta and/or the large column density in the line of sight. Later, X-ray pulsations were detected by Gotthelf et al. (2014) at the position of this compact object using the Nuclear Spectroscopic Telescope Array. Later on, Archibald et al. (2016) calculated the braking index of this pulsar using a 2.3 yr long phase-connected timing solution. With a characteristic

age of $τ_c = 3350$ yr and a spin-down power of $\dot{E} = 4.4 × 10^{36}$ erg s⁻<fnref n="1" />, the pulsar PSR J1640-4631 has similar properties to pulsars powering TeV pulsar wind nebulae (PWNe; Klepser et al. 2013). No radio counterpart was found, and only upper limits on the radio fluxes in this region were derived. No γ-ray or radio pulsation has been detected so far.

A spatially coincident GeV γ-ray source (1FGL J1640.8-465) was then revealed by the Fermi-Large Area Telescope (LAT) with a soft spectral index of $Γ = 2.3 ± 0.1$ (Slane et al. 2010). A distance of 8–13 kpc was derived based on 21 cm HI absorption by Lemiere et al. (2009), which implied that HESS J1640–465 was the most luminous VHE γ-ray source known in the Galaxy at that time. The later discovery of the neighboring VHE source HESS J1641–463 (Abramowski et al. 2014b) changed the GeV soft spectral index from 2.3 to a harder spectral index of $Γ = 1.99 ± 0.04$ (Lemoine-Goumard et al. 2014) from a dedicated analysis of both sources in the 0.1–300 GeV energy band. Using high-resolution multifrequency radio data of G338.3–0.0, Castelletti et al. (2011) derived upper limits on the radio flux from a potential extended radio nebula as the multifrequency radio observations did not detect radio pulsations toward the point source. More recently, Xin et al. (2018) obtained a hard index of $Γ = 1.42 ± 0.19$ during a dedicated analysis of HESS J1640–465 (from 10–500 GeV).

In addition to HESS J1640–465, there is a nearby (0.25°) weak source HESS J1641–463, which could only be detected by H.E.S.S. using an energy threshold of 4 TeV (Abramowski et al. 2014b). Characterized as a point-like source with a

<footnote n="1">¹ Université Bordeaux, CNRS, CENBG, UMR 5797, F-33170 Gradignan, France; mares@cenbg.in2p3.fr, lemoine@cenbg.in2p3.fr</footnote>

<footnote n="2">² AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="3">³ Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, The University of Manchester, M13 9PL, UK</footnote>

<footnote n="4">⁴ New York University Abu Dhabi, Abu Dhabi, United Arab Emirates</footnote>

<footnote n="5">⁵ Cavendish Laboratory, Cambridge CB3 0HE, UK</footnote>
spectral index of $\Gamma = 2.07 \pm 0.11$ and significant detection of photons up to 30 TeV, HESS J1641-463 is one of the best H.E. S.S. PeVatron candidates. Its location is coincident with the radio SNR G338.5+0.1 whose age is estimated between 1.1 and 3.5 kyr or between 5 and 17 kyr (Kothes & Dougherty 2007). The 843 MHz radio data from the Molongo Observatory Synthesis Telescope revealed a bright dense HII region (G338.4 +0.0) that connects in projection the two SNRs G338.3-0.0 and G338.5+0.1 (Whiteoak & Green 1996). No other known counterpart is found to be compatible with HESS J1641-463. Lemoine-Goumand et al. (2014) reported the discovery of a GeV source using 5 yr of the Fermi-LAT data coinciding with this TeV source. This analysis shows a spectral index of $\Gamma = 2.47 \pm 0.05$ for HESS J1641-463 in the 0.1-300 GeV energy range, which differs significantly from the spectral shape in the TeV energy range ($\Gamma \sim 2.1$).

At first, the detected TeV emission of HESS J1640-465 was interpreted as accelerated particles in the pulsar wind nebula (PWN)-emitting γ-rays through inverse Compton (IC) scattering. This scenario was reinforced by the detection of the extended X-ray source and the associated X-ray pulsar PSR J1640-4631. Xin et al. (2018) showed that the leptonic hypothesis of particle acceleration in a PWN was the best model to reproduce the γ-ray emission of HESS J1640-465. They modeled the particles with a power law in the high-energy part of their spectrum. The cutoff in the TeV energy band seen by H.E.S.S. was then easily explained by the Klein-Nishima regime.

Then, taking into account the radio data, the TeV emission of both neighboring sources could be explained by a self-consistent scenario in which protons might be accelerated up to hundreds of TeVs in either G338.5+0.1 or G338.3-0.0. In this case, the observed broadband spectrum of HESS J1640-465 would come from proton-proton interactions inside the dense HII region G338.4+0.1 located in the northern part of the SNR shell coincident with HESS J1640-465, as first proposed by Abramowski et al. (2014b) and then studied in more detail by Tang et al. (2015). Protons of higher energy could then escape from the accelerating zone in the SNR shock wave, travel toward the molecular cloud coinciding with HESS J1641-463 and interact with it. This would explain the relatively high brightness of the γ-ray emission from HESS J1641-463 in comparison with HESS J1640-465 above 4 TeV. This scenario is possible if both sources are at the same distance. This hypothesis is supported by the radio "bridge," detected by Whiteoak & Green (1996), connecting HESS J1640-465 and HESS J1641-463.

Another possibility is that HESS J1641-463 could be a binary system, with a pulsar-like spectrum at GeV energies that does not connect to the TeV spectrum coming from the wind shocks. However, no γ-ray variability has been detected using the Fermi-LAT data (Gotthelf et al. 2014), and no star is seen at any wavelength with a position close to HESS J1641-463.

This work presents an analysis of HESS J1640-465 and HESS J1641-463 using 8 yr of Fermi-LAT data. Section 2.1 describes the Fermi-LAT data used for the analysis, while Sections 2.2 and 2.3 present the results obtained from a morphological and spectral analyses of HESS J1640-465 and HESS J1641-463. Finally, in Section 3, we discuss the main implications of these results concerning the origin of the detected γ-ray signal. Our extensive analyses provide new constraints on the origin of the γ-ray emission as well as the

efficiency of these two teraelectron sources to accelerate protons and contribute to the Galactic cosmic-ray (CR) flux around the knee.

## 2. Fermi-LAT Analysis of HESS J1640-465 and HESS J1641-463

### 2.1. Fermi-LAT Observation and Data Analysis

The Fermi-LAT is an electron/positron pair conversion telescope launched in 2008, observing the γ-ray sky between 20 MeV to more than 500 GeV (Atwood et al. 2009). The following analyses take full advantage of the improvement in the quality of the data provided by the Pass 8 data upgrade (Atwood et al. 2013; Bruel et al. 2018). Using 8 yr of data, we performed morphological and spectral analyses above 1 GeV and 200 MeV, respectively, and up to 1 TeV for both sources within a $15^\circ \times 15^\circ$ region centered on HESS J1640-465. For these studies, we used all γ-rays, except those coming from a zenith angle above 90° to the detector axis to avoid Earth's limb contamination (Abdo et al. 2009). Time intervals when Fermi overflew the South Atlantic Anomaly were also removed. Version 11-07-00 of ScienceTools was used, together with the Fermi\'y (Wood et al. 2017) package version 0.18.0. We used the version P8R3_SOURCE_V2 of the instrument response functions (IRFs), released by the Fermi-LAT Collaboration, with the SOURCE event-class. This event-class is the most favorable for analysis of extended and point-like sources on medium to long timescales, while providing lower charged particle residual background contamination. Both analyses were made using *gtlike*. Included in Fermi ScienceTools,<fnref n="6" /> *gtlike* performs unbinned and binned maximum likelihood analyses (Mattox et al. 1996) by fitting a source model to the Fermi-LAT data. The diffuse emission in the Galaxy is described by the standard LAT diffuse model *gl>_iem_v07.fits* calculated using templates from GALPROP<fnref n="7" /> (Moskalenko et al. 2017). Extragalactic diffuse γ-rays, unresolved extragalactic sources, and residual CR emission are modeled by the isotropic spectral template in the analysis (*iso_P8R3_SOURCE_V2_v1.txt*).<fnref n="8" /> During this analysis, the normalization of the diffuse Galactic and isotropic emission were fitted. Our maximum likelihood analysis combines the four point-spread function (PSF) event types<fnref n="9" /> in a joint likelihood function. Using the catalog of the Fermi-LAT Collaboration (hereafter referred to as the 4FGL; Abdollahi et al. 2020) based on more than 8 yr of Fermi-LAT data, we took into account sources up to $20^\circ$ from HESS J1640-465. The spectral parameters of sources within $3.5^\circ$ of the center are left free, while the others are fixed at the value of the 4FGL catalog. Energy dispersion was taken into account for both morphological and spectral analyses.

<footnote n="6">⁶ The documentation of the Fermi ScienceTools is available at https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/.</footnote>

<footnote n="7">⁷ GALPROP is a numerical code for calculating the propagation of relativistic charged particles and the diffuse emissions produced during their propagation. The GALPROP code incorporates as much realistic astrophysical input as possible together with latest theoretical developments. The code calculates the propagation of CR nuclei, antiprotons, electrons, and positrons, and computes diffuse γ-rays and synchrotron emission in the same framework.</footnote>

<footnote n="8">⁸ Available from the Fermi Science Support Center: https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html.</footnote>

<footnote n="9">⁹ Each PSF event types indicate the quality of the reconstructed direction, the data is divided into quartiles, from the lowest quality quartile (PSF0) to the best quality quartile (PSF3).</footnote>
## 2.2. Morphological Analysis

As the PSF of the LAT improves with energy, we chose to perform the morphological analysis of HESS J1640-465 and HESS J1641-463 above 1 GeV to benefit from good angular resolution, without loosing too many photons. Two sources are detected in the 4FGL catalog (4FGL J1640.6-4632 and 4FGL J1641.0-4619), whose position coincides with HESS J1640-465 and HESS J1641-463. Figure 1 shows images of the sources in two energy bands (top: $E > 1$ GeV; bottom: $E > 10$ GeV) showing that HESS J1641-463 only appears below 10 GeV. Using *Fermi*, we have computed the associated residual test statistic (TS) map of the region, with a binning of 0.04°. The TS is defined as twice the log-likelihood difference between the model with a putative additional source and the reference model. It follows a chi-square distribution with *N* (number of parameters associated with the additional source) degrees of freedom assuming the reference model and the *N* parameters are nondegenerate. This TS map was obtained by testing, at the center of each pixel, the presence of a point-like source with a power-law spectrum with a fixed index of 2. The map contains the TS value for each bin and gives the significance of the detection in excess of the background at a given position. We then looked for potential new sources in the region by fitting the prefactors and the indices of the hotspots with a TS above 16. If the TS of a source was above 25 after the fit, we added them in the final source model of the region. Six sources were found by this method above 1 GeV. Two additional sources were added this way above 200 MeV (see Table 2). One of the added sources, PS J1703.2-4145, is close to a source (4FGL J1704.1-4140) detected in the second release of the 4FGL catalog (Ajello et al. 2020). It should be noted here that, unlike our analysis, the 4FGL introduced weights that give more importance to high-energy photons compared to the low-energy ones. This allows mitigating the effect of systematic errors due to our imperfect knowledge of the Galactic diffuse emission but decreases the TS value of sources with a soft spectrum. Since the sources reported in Table 2 have relatively soft spectra and a TS below 98, they could have been missed due to the use of weights in the 4FGL and 4FGL-DR2 catalogs. It might also be that these sources are artificial structures from the diffuse emission model in this region. A more careful analysis would be needed to confirm these soft and faint sources.

The spatial analysis of both sources was performed iteratively with *Fermi*. Using the updated source model, we relocalized HESS J1640-465 and HESS J1641-463 using a point-like source assumption. Then, we fixed the estimated position of HESS J1641-463 and looked for a potential extension of HESS J1640-465 (and computed its position). Finally, we used the new position and extension of HESS J1640-465 to relocalize HESS J1641-463 and search for its extension.

To find the best morphological model, the different spatial shapes were tested: (A) two point-like sources relocalized by *Fermi*, (B) a 2D Gaussian for HESS J1640-465 plus a point source for HESS J1641-463, relocalized by *Fermi*, and (C) the best positions and extension obtained by H.E.S.S. for both sources. For all these tests, a power law and a log-parabola spectral shape were assumed for HESS J1640-465 and for HESS J1641-463, respectively.

Table 1 gives the positions and extensions together with the associated number of degrees of freedom and the TS$_{ext}^{10}$ value.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Fermi-LAT residual TS maps of a 1.5° × 1.5° region around HESS J1640-465 above 1 GeV (top) and above 10 GeV (bottom). Sources from the 4FGL catalog are indicated with gray crosses (×). Our best localization of the point-source HESS J1641-463 is marked with the black plus (+), while the best localization and extension of the Gaussian HESS J1640-465 are represented by the blue circle. The red diamond indicates the position of the pulsar PSR J1640-4631. Radio continuum at 843 MHz is shown as cyan contours (Whiteoak & Green 1996). All sources are included in the background model except HESS J1640-465 and HESS J1641-463. The green circle indicates the position of the closest 4FGL source 4FGL J1639.8-4642 discussed in Section 2.3.</figcaption></figure>

This is the first time that the extension of HESS J1640-465 is significantly measured above 1 GeV using the Fermi-LAT data with a TS$_{ext}$ = 35, which corresponds to a significance of ~6σ. The best-fit position found using *Fermi* for HESS J1640-465 as a 2D Gaussian is $\alpha(J2000) = 250.17° \pm 0.01°_{\text{stat}} \pm 0.01°_{\text{syst}}$, $\delta(J2000) = 46.55° \pm 0.01°_{\text{stat}} \pm 0.01°_{\text{syst}}$ with an extent of $\sigma = 0.053° \pm 0.011°_{\text{stat}} \pm 0.008°_{\text{syst}}$. The extension measured by H.E.S.S. of 0.072° ± 0.003° (Abramowski et al. 2014a) is consistent with our best estimate, considering the uncertainties on both measurements.

<footnote n="10">10 TS$_{ext}$ is defined as twice the log-likelihood difference between an extended model of the source and the point-like source model.</footnote>
<table><caption>**Table 1**

Best Morphological Parameters (Positions and Extension) for HESS J1640–465 and HESS J1641–463</caption><thead><tr><th></th><th></th><th>Fermipy Position (A)</th><th>Fermipy Best Fits (B)</th><th>H.E.S.S. Position/Extension (C)</th></tr></thead><tbody><tr><td rowspan="3">HESS J1640–465</td><td>R.A. (°)</td><td>250.18 ± 0.01<sub>stat</sub></td><td>250.17 ± 0.01<sub>stat</sub> ± 0.01<sub>syst</sub></td><td>250.180</td></tr><tr><td>Decl. (°)</td><td>-46.55 ± 0.01<sub>stat</sub></td><td>-46.55 ± 0.01<sub>stat</sub> ± 0.01<sub>syst</sub></td><td>-46.530</td></tr><tr><td>σ (°)</td><td>...</td><td>0.053 ± 0.011<sub>stat</sub> ± 0.008<sub>syst</sub></td><td>0.072</td></tr><tr><td rowspan="2">HESS J1641–463</td><td>R.A. (°)</td><td>250.26 ± 0.01<sub>stat</sub></td><td>250.25 ± 0.01<sub>stat</sub> ± 0.01<sub>syst</sub></td><td>250.26</td></tr><tr><td>Decl. (°)</td><td>-46.32 ± 0.01<sub>stat</sub></td><td>-46.29 ± 0.01<sub>stat</sub> ± 0.03<sub>syst</sub></td><td>-46.30</td></tr><tr><td>Degree of freedom (DoF)</td><td></td><td>9</td><td>10</td><td>5</td></tr><tr><td>TS<sub>ext</sub> (point versus model)</td><td></td><td>...</td><td>35</td><td>34</td></tr></tbody></table>

No significant extension was found for HESS J1641–463 using *Fermipy* (with a TS<sub>ext</sub> = 6 for a 2D Gaussian, equivalent to 2.45σ). The best-fit position for HESS J1641–463 as a point-like source is α(J2000) = 250.25° ± 0.01°<sub>stat</sub>, δ(J2000) = 46.29° ± 0.01°<sub>stat</sub>. The good agreement between the estimated positions of HESS J1640–465 and HESS J1641–463 in the GeV band by *Fermipy* (B) and in the TeV range covered by H.E.S.S. (C) tends to show that H.E.S.S. and the Fermi-LAT detect the same sources.

We used the best positions and extension obtained by H.E.S. S. for HESS J1640–465 and HESS J1641–463 (model C) to derive their spectra, since the best model obtained by *Fermipy* is not significantly better than the H.E.S.S. extension and position values extracted from Abramowski et al. (2014a, 2014b).

## 2.3. Spectral Analysis

The following spectral analysis was performed using *gtlike* in a summed likelihood binned analysis with a combination of the four PSF event types between 200 MeV up to 1 TeV. The source model of the region contains the eight additional sources presented in Table 2, as well as the best morphological parameters from H.E.S.S. for HESS J1641–463 and HESS J1640–465. In order to obtain the best model of the observed emission, we made a comparison between multiple spectral shapes for both sources (see Table 3).

Assuming the best spectral shape, the γ-ray sources observed by the LAT above 200 MeV are detected with a TS of 232 for HESS J1640–465 and 813 for HESS J1641–463. The spectrum of HESS J1641–463 is well described by a log-parabola ($dN/dE \propto E^{-(\alpha+\beta \log(E/E_b))}$) with $\alpha = 2.7 \pm 0.1_{\text{stat}} \pm 0.2_{\text{syst}}$ and a significant curvature ($4\sigma$ improvement with respect to a power-law model) of $\beta = 0.11 \pm 0.03_{\text{stat}} \pm 0.05_{\text{syst}}$. An integrated flux of ($4.54 \pm 0.31_{\text{stat}} \pm 1.53_{\text{syst}}$) $\times 10^{-11}$ erg cm$^{-2}$ s$^{-1}$ is obtained between 200 MeV and 1 TeV.

The spectrum of HESS J1640–465 is well described by a power law ($dN/dE \propto E^{-\Gamma}$) with a hard spectral index of $\Gamma = 1.8 \pm 0.1_{\text{stat}} \pm 0.2_{\text{syst}}$ and integrated flux between 200 MeV and 1 TeV of ($40.0 \pm 0.5_{\text{stat}} \pm 0.3_{\text{syst}}$) $\times 10^{-11}$ erg cm$^{-2}$ s$^{-1}$. The spectrum does not show any significant curvature (see Table 3).

The results presented here are given with statistical and systematic uncertainties. The systematic uncertainties are due to our imperfect model of the Galactic diffuse emission (Err<sub>adm</sub>), to the spatial shape of HESS J1640–465 (Err<sub>shape</sub>) and to our knowledge of the Fermi-LAT IRFs (Err<sub>irf</sub>). The systematic uncertainties are added in a quadrature as Err<sub>syst</sub> = $\sqrt{(Err_{\text{adm}})^2 + (Err_{\text{irf}})^2 + (Err_{\text{shape}})^2}$.

The first uncertainty was estimated using the results of fits taking into account eight alternative interstellar diffuse emission models and comparing the results to the nominal fit, as in Acero et al. (2016). The uncertainties on the morphology

<table><caption>**Table 2**

Names, Spectral Indices, and TS of the Additional Sources Found by *Fermipy* in the Region of 15° × 15° Centered on HESS J1640–465 above 200 MeV (†) and above 1 GeV</caption><thead><tr><th>Name</th><th>Index</th><th>TS</th></tr></thead><tbody><tr><td>PS J1617.5–5104<sup>†</sup></td><td>Γ = 2.02 ± 0.02</td><td>37</td></tr><tr><td>PS J1708.8–4007<sup>†</sup></td><td>Γ = 3.02 ± 0.02</td><td>55</td></tr><tr><td>PS J1633.7–4755</td><td>Γ = 2.71 ± 0.17</td><td>98</td></tr><tr><td>PS J1638.4–4715</td><td>Γ = 2.62 ± 0.02</td><td>33</td></tr><tr><td>PS J1640.2–4803</td><td>Γ = 3.26 ± 0.34</td><td>34</td></tr><tr><td>PS J1647.4–4541</td><td>Γ = 2.22 ± 0.17</td><td>32</td></tr><tr><td>PS J1652.3–4433</td><td>Γ = 2.84 ± 0.02</td><td>54</td></tr><tr><td>PS J1703.2–4145</td><td>Γ = 3.45 ± 0.40</td><td>59</td></tr></tbody></table>

of HESS J1640–465 was estimated by fitting the γ-ray emission with the best spatial models given by *Fermipy*. The third uncertainty was estimated by using modified IRFs whose effective areas bracket the nominal ones.<sup>11</sup> The same uncertainties were taken into account for HESS J1641–463, but without the uncertainties on the morphology Err<sub>shape</sub>.

Figures 2 and 3 present the spectra obtained for HESS J1640–465 and HESS J1641–463 by dividing the 0.2–1000 GeV interval into 11 logarithmically spaced energy bins. A maximum likelihood spectral analysis was performed in each interval, considering both sources as a power-law shape of index $\Gamma = 2$. Only the prefactor of the sources within 3.5° from HESS J1640–465, and the normalization of the diffuse Galactic and isotropic emission were left free. When the statistical significance was below $1\sigma$, a 95% confidence level upper limit was computed. As we did not detect any signal with the Fermi-LAT data in the 0.1–1 TeV band for HESS J1641–463, an upper limit was obtained for the combination of the three last bins. Both statistical (red) and statistical plus systematic (black) uncertainties are shown in Figures 2 and 3.

There are several differences between this analysis and the previous analysis reported in Lemoine-Goumard et al. (2014): the use of 8 yr of Pass 8 data, the new 4FGL source list, updated Galactic and isotropic diffuse emission models, and a curved spectral shape for HESS J1641–463. Furthermore, the lower energy range is different between the two analysis, which as a strong impact on the first point of the spectral energy distribution (SED) (from UL to flux). Two hypotheses were tested in order to recover the results of the previous analysis of 2014. The first test was to remove the curvature from HESS J1641–463. Doing so, the index of HESS J1640–465 became softer (with an index of $\Gamma = 1.85$) but not enough to retrieve the index of ~2. The biggest difference between the two analyses came from the neighboring source 4FGL J1639.8-4642 (see

<footnote n="11">11 https://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html</footnote>
<table><caption>**Table 3**
Comparison of Spectral Models (200 MeV–1 TeV)</caption><thead><tr><th>HESS J1640–465</th><th>HESS J1641–463</th><th>Index<sup>a</sup></th><th>β</th><th>DoF</th><th>ΔTS<sub>modelvs PLPL</sub></th></tr></thead><tbody><tr><td>Power law</td><td>Power law</td><td>1.7 ± 0.1/2.7 ± 0.06</td><td>–/–</td><td>4</td><td>...</td></tr><tr><td>Power law</td><td>Log-parabola</td><td>1.8 ± 0.1/2.7 ± 0.1</td><td>–/0.11 ± 0.03</td><td>5</td><td>16</td></tr><tr><td>Log-parabola</td><td>Power law</td><td>1.7 ± 0.1/2.7 ± 0.08</td><td>0.0003 ± 0.005/–</td><td>5</td><td>2</td></tr><tr><td>Log-parabola</td><td>Log-parabola</td><td>1.7 ± 0.2/2.6 ± 0.1</td><td>0.0004 ± 0.005/0.17 ± 0.05</td><td>6</td><td>18</td></tr></tbody></table>

**Note.**

<figure><img src="image_3.png" /><figcaption>**Figure 2.** γ-ray SED of HESS J1640–465 between 200 MeV and 1 TeV with the Fermi-LAT data, using the best Gaussian spatial model from H.E.S.S. The purple thin line and the open diamonds are from the latest publication by the Fermi-LAT Collaboration (Lemoine-Goumard et al. 2014). The gray line and open stars are from Xin et al. (2018). The red line and dashed area represent our best-fit and 68% confidence band of the fitted Fermi-LAT spectrum. The red data points (circle) are the flux estimated in each energy bins together with their statistic (red) and statistic plus systematic (black) uncertainties. A 95% confidence level upper limit (UL) is computed when the significance in an energy bin is lower than 1. The blue diamonds and line represent the spectrum from the H.E.S.S. analysis extracted from Abramowski et al. (2014b). The thin blue line and the open crosses are from the HESS Galactic Plane Survey (Abdalla et al. 2018a).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** γ-ray SED of HESS J1641–463 between 200 MeV and 1 TeV with the Fermi-LAT data, using the best H.E.S.S. localization, using the same conventions as in Figure 2.</figcaption></figure>

Figure 1). This source is 0.2° away from HESS J1640–465 and was not detected in the 3FGL source list used in the previous analysis. Removing this source from the source model of the region degraded the log-likelihood value of the fit by

ΔTS = 102. Not taking this source into account results in a soft spectral index for HESS J1640–465 of ≈2.1. The log-parabola shape of this nearby source, with a significant curvature of β = 0.15 ± 0.058, probably contaminated the emission of HESS J1640–465 at low energy, yielding the softer<fnref n="12" /> spectral index of ~2 reported in the previous analysis. This was also the case in the analysis of Slane et al. (2010) with additional contamination coming from HESS J1641–463, separated by 0.25° as well, giving a spectral index of ~2.3 for HESS J1640–465.

### 3. Multiwavelength Context

#### 3.1. X-Ray Counterpart

The region has been observed both with XMM-Newton and the Chandra telescopes for follow-up observation of HESS J1640–465. In the XMM-Newton observation (ObsID 0302560201; exposure of 22 ks), HESS J1641–463 lies in the border of the field of view and its upper limit of 3' on its extent is not fully covered by the camera. In the Chandra observation (ObsID 12508; exposure of 19 ks), the source and its upper limit on its extent is fully contained in the ACIS-S3 chip and the analysis that follows is performed with the Chandra data only.

As shown in Figure 4, no significant sources are detected within the source extent, and we computed a conservative upper limit assuming an extended source based on the HESS position, α(J2000) = 16h41m2.1s, δ(J2000) = −46°18′13″, and using the H.E.S.S. upper limit of 3' as the extent. To compare with other multiwavelength measurements, the X-ray upper limit requires an assumed value of the X-ray absorption (N<sub>H</sub>) along the line of sight. As this N<sub>H</sub> value and the distance to the source are unknown, we estimated the upper limit above 2 keV to limit this systematic effect as above this energy the X-ray absorption effect becomes negligible. For the upper bound, the energy range was limited to 6 keV as above this energy the instrumental background becomes dominant. The astrophysical and instrumental background were derived using an OFF region (shown in Figure 4) taken from the same chip.<fnref n="13" /> In this restricted energy range, the OFF region is well modeled with two power laws. The ON region is modeled with an absorbed power law with a spectral index fixed to 2 and the N<sub>H</sub> value is taken to the Galactic HI value at this position (2 × 10<sup>22</sup> cm<sup>−2</sup>). As the unabsorbed flux is estimated above 2 KeV, the choice in the N<sub>H</sub> value has a negligible impact on the flux upper limit. With the models described above, the upper limit is derived by computing a likelihood profile on the

<footnote><sup>a</sup> The left value is for HESS J1640–465 and the right value is for HESS J1641–463. The Index column is for both the index (power law), and for α (log-parabola).</footnote>

<footnote n="12">¹² Compared to this analysis.</footnote>

<footnote n="13">¹³ As the S3 chip is back illuminated, the instrumental background is slightly different and cannot be estimated from the other front-illuminated chips in this observation.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Counts map in the 0.5–6 keV energy band of the Chandra observation smoothed with a 6" kernel. The scale bar shows photons per pixel. The circle is centered at the best-fit position for HESS J1641–463, with a radius that shows the upper limit on its extent as derived in the H.E.S.S. Galactic Plane Survey. The definition of the OFF region used to estimate the astrophysical and instrumental background is also shown.</figcaption></figure>

flux jointly fitting the flux and the background from the OFF region. The resulting 95% confidence unabsorbed flux upper limit is $1.58 \times 10^{-13} \text{ erg cm}^{-2} \text{ s}^{-1}$.

## 3.2. Radio Counterpart

The Southern Galactic Plane Survey (SGPS; Haverkorn et al. 2006) covers G338.5+0.1 at 1.4 GHz, with a resolution of 100". The SGPS image does not show any obvious local maximum at the position of HESS J1641–463. In order to constrain the level of any radio emission associated with HESS J1641–463, Gaussian models with varying amplitudes, each with an FWHM of 3', were added to the observed SGPS image. From these an upper limit of ~2 Jy for any radio emission was estimated.

## 4. Discussion

In this paper, we have investigated the GeV emission from HESS J1640–465 and HESS J1641–463. Multiple scenarios can explain the γ-ray emission from these sources. These scenarios include the emission from accelerated particles within the SNRs coincident with each source, molecular clouds illuminated by cosmic rays (CRs) from close-by SNR (either from the SNR G338.3–0.0 for HESS J1640–465 or the SNR G338.5-0.1 for HESS J1641–463), a PWN (detected in X-ray for HESS J1640–465 and hypothetical for HESS J1641–463), and an active galactic nucleus (AGN) for HESS J1641–463. These scenarios are discussed below.

### 4.1. HESS J1640–465

The detection of the extension of HESS J1640–465 in the GeV energy band using the Fermi-LAT data, together with the

new spectrum, help us constrain the origin of its emission. The hard spectrum reported in this work smoothly connects with the TeV spectral points extracted from Abramowski et al. (2014a) as shown in Figure 2. Previous studies of this source have interpreted its high-energy emission as being produced by high-energy leptons accelerated in a PWN (e.g., Slane et al. 2010; Abramowski et al. 2014b) or high-energy hadrons accelerated in the associated SNRs. In Section 4.1.1, we analyze the properties of this source assuming a PWN origin for the observed γ-rays, while in Section 4.1.3 we assume that the observed γ-rays are primarily produced as the result of interactions between high-energy protons or electrons accelerated in the SNR and the surrounding medium.

#### 4.1.1. Pulsar Wind Nebula Scenario

Since the first detection of γ-ray emission from the Crab Nebula (Fazio et al. 1972), numerous γ-ray sources have been associated with PWNe (e.g., Klepser et al. 2013). In these sources, the rotational energy of a strongly magnetized, rapidly rotating neutron star (i.e., a pulsar) powers a highly relativistic $e^{\pm}$ wind (e.g., Goldreich & Julian 1969). The confinement of this “pulsar wind” by the surrounding medium creates a PWN (e.g., Pacini & Salvati 1973), inside which leptons are accelerated to high energies (e.g., Nishikawa et al. 2005; Spitkovsky 2008). In such systems, γ-rays are believed to be produced through IC scattering from these energetic particles. For the clarity of this section, the details on the equations and on the modeling were placed in the Appendix.

When a pulsar is young it resides in the SNR created by the interaction between the material ejected in the progenitor explosion and the surrounding medium (e.g., Kennel & Coroniti 1984). During this time, the evolution of its PWN is sensitive to the spin-down evolution of the pulsar, the content of the pulsar wind, and the evolution of the surrounding SNR (e.g., Reynolds & Chevalier 1984). As a result, the properties of such a PWN can be used to determine the initial period of the associated neutron star, the magnetization and particle spectrum of its pulsar wind, and the mass and initial kinetic energy of the supernova explosion. Currently, this is best done by reproducing the dynamical and radiative properties of a PWN inside an SNR with a model for its evolution, as described in a recent summary by Gelfand (2017). This has been successfully done for numerous systems (e.g., Chevalier 2005; Bucciantini et al. 2011; Tanaka & Takahara 2011; Torres et al. 2014; Zhu et al. 2018)—including HESS J1640–465 using previous measurements of its γ-ray spectrum (e.g., Gotthelf et al. 2014). As described in detail in the Appendix, we determined if such models can reproduce the updated γ-ray spectrum derived in Section 2.3. The combination of model parameters that best fit the observed properties of this source ($\chi^2 \approx 38$, or reduced $\chi^2 \approx 1.9$ given the 20 degrees of freedom) are given in Table 4, with the predicted values for the observed quantities for this combination also given in Table 6, in the Appendix. As shown in Figure 5, the SED predicted by our model provides an excellent fit to the observed quantities, especially its γ-ray emission.

While our model produces a statistically good fit to the observed properties of this source, determining if it is plausible requires considering the resultant physical properties, as listed in Tables 4 and 5. While our fit suggests the mass ($M_{ej} \approx 10 M_{\odot}$) and initial kinetic energy ($E_{sn} \approx 10^{51}$ erg) of the supernova ejecta are close to “typical” values (e.g.,
<table><caption>**Table 4**
"Best-fit" Model Parameters for HESS J1640-465</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Supernova explosion energy E<sub>sn</sub></td><td>1.4 × 10<sup>51</sup> erg</td></tr><tr><td>Supernova ejecta mass M<sub>ej</sub></td><td>11 M<sub>☉</sub></td></tr><tr><td>ISM density n<sub>ism</sub></td><td>0.009 cm<sup>-3</sup></td></tr><tr><td>Spin-down timescale τ<sub>sd</sub></td><td>4.15 yr</td></tr><tr><td>Magnetization η<sub>B</sub></td><td>0.10</td></tr><tr><td>Minimum energy of injected e<sup>±</sup>, E<sub>min</sub></td><td>270 GeV</td></tr><tr><td>Break energy of injected e<sup>±</sup>, E<sub>break</sub></td><td>16 TeV</td></tr><tr><td>Maximum energy of injected e<sup>±</sup>, E<sub>max</sub></td><td>0.3 PeV</td></tr><tr><td>Low-energy index of injected e<sup>±</sup> p<sub>1</sub></td><td>1.94</td></tr><tr><td>High-energy index of injected e<sup>±</sup> p<sub>2</sub></td><td>2.94</td></tr><tr><td>Number of added IC photon field</td><td>≈ 2</td></tr><tr><td>Temperature of added photon fields</td><td>5.4 K, 27,000 K</td></tr><tr><td>Norm of added photon fields</td><td>0.16, 5.2 × 10<sup>-13</sup></td></tr><tr><td>Distance d</td><td>12.1 kpc</td></tr></tbody></table>

Note. The observed properties of HESS J1640-465 used in this fit are listed in Table 6.

Baade & Zwicky 1934), it requires that the resultant pulsar (PSR J1640-4631) has an extremely short spin-down timescale<fnref n="14" />: $\tau_{sd} \approx 4$ yr instead of $\tau_{sd} \sim 1000$ yr inferred for other pulsars associated with PWNe (see the Appendix) (e.g., Bucciantini et al. 2011; Tanaka & Takahara 2011; Gelfand et al. 2014, 2015; Torres et al. 2014; Zhu et al. 2018). If correct, then the pulsar was born spinning very rapidly, with an initial spin period (e.g., Goldreich & Julian 1969; Pacini & Salvati 1973; Gaensler & Slane 2006):

$$P_0 = P \left( 1 + \frac{t_{\text{age}}}{\tau_{\text{sd}}} \right)^{-\frac{1}{p-1}} \approx 9.5 \text{ ms,} \quad (1)$$

for a current period of $P \approx 206$ ms (Gotthelf et al. 2014). The initial period derivative of this pulsar $\dot{P}_0$ predicted by this solution is

$$\dot{P}(t) = \frac{P_0}{\tau_{\text{sd}}(p-1)} \left(1 + \frac{t}{\tau_{\text{sd}}}\right)^{\frac{2-p}{p-1}}, \quad (2)$$

$$\dot{P}_0 = \frac{P_0}{\tau_{\text{sd}}(p-1)} \approx 3.38 \times 10^{-11} \frac{\text{s}}{\text{s}}, \quad (3)$$

which suggests as initial surface dipole magnetic field strength $B_{s0} \equiv 3.2 \times 10^{19} \sqrt{P_0 \dot{P}_0} \approx 1.8 \times 10^{13}$ G—only ~30% larger than the $B_s \approx 1.4 \times 10^{13}$ G inferred from its current spin-down properties (Gotthelf et al. 2014). However, this model requires this pulsar had an extremely high initial spin-down luminosity $\dot{E}_0 \approx 1.55 \times 10^{42}$ erg s⁻¹, ~100× higher than that estimated in analyses of other PWNe (e.g., Torres et al. 2014). However, the total rotational energy injected into the SNR at early times, $\dot{E}_{0sd} \approx 2.0 \times 10^{50}$ erg, is considerably lower than the $\gtrsim 10^{51}$ erg explosion energy. Lastly, the magnetization of the wind produced by this pulsar $\eta_B \sim 10^{-1}$ is significantly larger than the magnetization $\eta_B \sim 10^{-3}$ inferred from many other pulsars (e.g., Bucciantini et al. 2011; Tanaka & Takahara 2011; Torres et al. 2014; Gelfand et al. 2015; Zhu et al. 2018), but is comparable to the $\eta_B \sim 0.04$ inferred by a similar analysis of

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Observed (color points) and predicted (black line) SED of the PWN in HESS J1640-465. The model parameters used in the prediction are listed in Table 4.</figcaption></figure>

<table><caption>**Table 5**
Physical Properties of HESS J1640-465 Predicted by Our Model with the Parameters Listed in Table 4</caption><tr><td>True age t<sub>age</sub></td><td>3094 yr</td></tr><tr><td>Initial spin period P<sub>0</sub></td><td>9.5 ms</td></tr><tr><td>Initial period derivative P&#775;<sub>0</sub></td><td>3.38 &times; 10<sup>-11</sup> s s<sup>-1</sup></td></tr><tr><td>Initial surface magnetic field B<sub>s0</sub></td><td>1.8 &times; 10<sup>13</sup> G</td></tr><tr><td>Initial spin-down luminosity E&#775;<sub>0</sub></td><td>1.5 &times; 10<sup>42</sup> erg s<sup>-1</sup></td></tr><tr><td>Energy density of added photon fields</td><td>0.64 eV cm<sup>-3</sup>, 1.3 keV cm<sup>-3</sup></td></tr></table>

Kes 75 (Gelfand et al. 2014), whose associated neutron star also has $B_s \gtrsim 10^{13}$ G. It is also worth noting that PSR J1640-4631 is one of the few pulsars with a measured braking index, $3.15 \pm 0.03$, and the first one with $p > 3$, the canonical value for pure magnetic dipole radiation (Archibald et al. 2016), measured from a 2.3 yr long phase-connected timing solution. As such, it is perhaps not surprising that it would have an initial spin period $P_0$ and initial spin-down luminosity $\dot{E}_0$ faster and higher, respectively, than most neutron stars.

In addition, our model requires that this SNR be expanding into an environment with an extremely intense UV photon field, suggesting it is near a collection of massive stars capable of producing the $T \sim 25,000$ K photon field required by our model to reproduce the GeV emission from this source (Table 4). In fact, a cluster of young massive stars Mc81 is ~8' away from HESS J1640-465 (Davies et al. 2012), with a distance ($d = 11 \pm 2$ kpc) and high extinction ($A_V \sim 45$), which are consistent with the distance $d \sim 12$ kpc favored by our modeling of the PWN (Table 4) and the considerable $N_H \sim (1-2) \times 10^{23}$ cm⁻² along the line of sight to this source inferred from its X-ray spectrum (Gotthelf et al. 2014). Since the projected distance between this stellar cluster and HESS J1640-465 is ~30 pc, this stellar cluster is close enough to provide the intense photon field our modeling required to reproduce the observed γ-ray spectrum.

If the progenitor of HESS J1640-465 was originally a member of this cluster but dynamically ejected before it exploded, as suggested by Davies et al. (2012), it likely was its most massive member since it was the first to explode. In this case, its initial mass should have been $\gtrsim 65 M_{\odot}$ inferred for a Wolf-Rayet star in this cluster (Mc81-3; Davies et al. 2012). If

<footnote n="14">¹⁴ $\tau_{sd}$ is the spin-down timescale, equivalent to the “characteristic age” of the pulsar at birth.</footnote>
correct, the $M_{ej} \sim 11 M_{\odot}$ of material ejected in this explosion (Table 4) suggests the progenitor lost most of its mass before exploding. In this case, its powerful stellar winds would create a low-density bubble into which the ejecta initially expands. This is supported by the low interstellar matter (ISM) density ($n_{ism} \sim 0.01 \text{ cm}^{-3}$) preferred by our modeling (Table 4). While seemingly at odds with the dense molecular material detected nearby (e.g., Lau et al. 2017), the short lifetime of such massive stars suggests they explode in the vicinity of their birth molecular cloud—consistent with the suggested interaction between the SNR shell and surrounding molecular material.

### 4.1.2. Pulsar Scenario

The significant extension detected by our analysis above 1 GeV, in agreement with that previously reported by the H.E. S.S. Collaboration, together with our hard spectral index of 1.8 strongly disfavor the pulsar hypothesis for HESS J1640–465. The fact that no γ-ray pulsations have been detected coming from an X-ray pulsar with such a high $\dot{E}$ (4.4 × 10$^{36}$ erg s$^{-1}$) could be explained as in Johnston et al. (2020). These X-ray pulsars have a selection bias making them at great distance, deep in the Galactic plane, and thus below LAT sensitivity. Furthermore, this pulsar is co-located with HESS J1640–465 and HESS J1641–463, making the detection of any signal from this pulsar even more difficult. On top of that, their geometry may be particularly unfavorable. Finally, an X-ray timing solution of 10 or 12 yr is necessary for the faintest pulsars (Smith et al. 2019), which is far from the 2.3 yr timing solution used for PSR J1640-4631.

### 4.1.3. SNR Scenario

Now, we consider the hypothesis that the γ-ray emission detected from HESS J1640–465 is produced by particles accelerated at the shocks of the SNR G338.3–0.0. Since an SNR origin for such emission has been proposed by a number of authors (Abramowski et al. 2014a, 2014b; Tang et al. 2015; Lau et al. 2017; Supan et al. 2016), we limit the discussion to a few essential issues that emerged from the new observational results presented in this paper.

In the GeV domain, the γ-ray spectrum of HESS J1640–465 can be described by a power law with a slope of 1.8, which is significantly harder than what was previously reported (Lemoine-Goumand et al. 2014). High-energy γ-ray spectra harder than $E^{-2}$ are often interpreted as the result of IC scattering emission from electrons accelerated at the SNR shock. This is the case, for example, for the SNR RX J1713.7–3946, (Abdo et al. 2011), whose TeV spectrum very closely resembles that of HESS J1640–465 (Abramowski et al. 2014a). An estimate of the energy required in accelerated electrons can be obtained by analogy with RX J1713.7–3946. Porter et al. (2006) showed that a leptonic (IC) fit to γ-ray data from RX J1713.7–3946 would require an energy of $\approx 5 \times 10^{47}$ erg in the form of accelerated electrons. This is a small fraction ($\approx 2 \times 10^{-3}$) of the typical supernova explosion energy ($\sim 10^{51}$ erg). The TeV lux of HESS J1640–465 is a factor of $\lesssim 10$ smaller than that of RX J1713.7–3946, but its distance is 10 times larger. This implies that a leptonic interpretation of the γ-ray spectrum of HESS J1640–465 would require a total energy in relativistic electrons equal to $\gtrsim 5 \times 10^{48}$ erg. In other words, if the SNR shock converts the canonical $\sim 10\%$ fraction of the explosion energy into CR protons, the leptonic

interpretation would require an electron-to-proton ratio of the order of $K_{ep} \gtrsim 5 \times 10^{-2}$. Such a large value of the electron-to-proton ratio, though not excluded, is significantly larger than what is usually assumed in modeling the broadband emission from SNRs, i.e., $K_{ep} \approx 10^{-5}...10^{-2}$ (see, e.g., Cristofari et al. 2013). A similarly large value of $K_{ep} \approx 0.1$ was also found by Abramowski et al. (2014a) (see their Figure 5, where they fit the multiwavelength spectrum adopting a time-dependent single-zone model).

On the other hand, as pointed out by Zirakashvili & Aharonian (2010) and Inoue & Yamazaki (2012), an alternative, hadronic scenario may also be invoked to explain the hard γ-ray spectra that are observed from a number of SNRs. This is because, if the SNR shock expands in a clumpy medium, propagation effects can shape the γ-ray spectrum. The reason is that, due to the finite age of the SNR, only CRs of sufficiently large energy have enough time to diffusively penetrate into dense clumps, interact with the gas, and produce gamma-rays from neutral pion decay. The exclusion of lower energy CRs from the clumps will result in a very hard γ-ray spectrum (Gabici & Aharonian 2014; Celli et al. 2019) that might easily match the observed spectrum of HESS J1640–465. Note that a more conventional value for the electron-to-proton ratio $K_{ep} \approx 10^{-2}$ was adopted in the hadronic scenario described in Abramowski et al. (2014a).

## 4.2. HESS J1641–463

A possibility is that HESS J1641–463 is an AGN, given its point-like morphology and curved spectrum seen in the GeV energy band. If it is a blazar, the spectral index of 2.7 (for a log-parabola spectrum) indicates that it should be a flat spectrum radio quasar (FSRQ). With a TS of 316 in the 4FGL catalog, a FSRQ will have a probability >0.95 to be seen as variable in the 4FGL and a probability of almost 90% to present a power-law spectrum instead of a log-parabola spectrum (see Figures 8 and 12 in Ajello et al. 2020). In the second release of the 4FGL catalog, the variability index<fnref n="15" /> is $\sim 7.4$. Therefore, no variability is detected in the 4FGL-DR2 using 10 yr of data. Whereas the FSRQ hypothesis cannot be ruled out, it is unlikely. Instead we consider SNR, pulsar, and PWN hypotheses for HESS J1641–463.

### 4.2.1. Hadronic Scenario (SNR or Illuminated Molecular Cloud Scenario)

Interpreting the γ-ray spectrum of HESS J1641–463 is a more difficult task than for HESS J1640–465, because it is unusual in having a curved spectrum at GeV energies, and a hard spectrum at TeV energies. The curved spectrum observed by Fermi could result either from the direct acceleration of CR protons at the shock of an aged SNR (Lau et al. 2017), or from the reacceleration of ambient CRs, if the shock interacts with a cloud (i.e., cloud crushed model, see Uchiyama et al. 2010). In the latter case, the shock becomes radiative and compresses (crushes) simultaneously the gas and the accelerated CRs.

As a result, both the γ-ray and radio emission are enhanced significantly. However, the very large ratio between γ-ray and radio emission (see radio upper limit in Figure 6) casts serious doubts on the viability of such a scenario. Concerning an

<footnote n="15">¹⁵ Parameter constructed from the value of the likelihood in the null hypothesis, i.e., the flux of the source is constant over the time, and the value under the hypothesis where the flux in each bin is optimized.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Multiwavelength SED and γ-ray spectrum of HESS J1641-463 between 200 MeV and 1 TeV with the Fermi-LAT data, using the best H.E.S.S. localization, using the same conventions as in Figure 3. The estimated 95% confidence X-ray flux upper limit between 2 and 6 keV and the radio flux upper limit at 1.4 GHz (see Sections 3.1 and 3.2 for more details) are also shown.</figcaption></figure>

hadronic origin of the emission, the ratio between the radio flux upper limit and the GeV flux is in agreement with what is observed for other SNRs (Figure 12 in Acero et al. 2016). But the X-ray upper limit also constraints severely the scenario of direct acceleration at the SNR shock. Indeed, the large ratio between the X-ray upper limit and the GeV flux (GeV/X-ray ~ 100) is problematic for a hadronic emission. This is because, due to the large gas density in the region, a bright thermal X-ray emission would be expected (Katz & Waxman 2008). Such an emission is not observed here. For these reasons, a non-SNR origin of the GeV emission remains a very plausible hypothesis.

Finally, the hard TeV spectrum shown by HESS J1641-463 has been often interpreted as the result of the illumination of ambient dense gas by CRs that escaped from the SNR G338.3-0.0 (Abramowski et al. 2014b; Tang et al. 2015; de Wilt et al. 2017; Lau et al. 2017). In order to test the validity of this picture, a detailed knowledge of the three-dimensional distribution of the gas in the region is needed (see, e.g., Pedaletti et al. 2015). In this context, de Wilt et al. (2017) performed an accurate study of various gas tracers and concluded that, though possible, the illumination scenario would imply that the peak of the TeV γ-ray emission should correlate with the peak observed in the gas distribution, which is not observed. Another prediction of the illumination scenario is the presence of a faint X-ray emission, due to the synchrotron emission from the secondary electrons produced in the interactions of the runaway CRs with the dense gas (Gabici et al. 2009). However, this is expected to be faint and could be absorbed. Large ratios between the γ-ray and X-ray fluxes are indeed expected, independently on the choice of model parameters (see the discussion in Gabici et al. 2009), and thus this scenario is consistent with the X-ray upper limit reported here (see Figure 6). In this case, the GeV spectrum of HESS J1641-463 would remain unexplained.

### 4.2.2. Pulsar Scenario

The curved GeV spectrum of HESS J1641-463 reported in Section 2.3 is similar to other γ-ray pulsars seen with Fermi-LAT (Abdo et al. 2013). We use the spectral parameters

obtained here to calculate the γ-ray luminosity of the source in the GeV energy band assuming it is a pulsar.

We first attempt to fit the spectrum of HESS J1641-463 using a power law with a super exponential cutoff($dN/dE \propto E^{-\Gamma_1} \times \exp(-(E/E_{cutoff})^{\Gamma_2})$) generally used to characterize pulsar spectra. The best fit yields $\Gamma_1 = 2.1 \pm 0.1$, $\Gamma_2 = 0.34 \pm 0.12$ and $E_{cutoff} = 573.4 \pm 15.1$ MeV. Using the value of the energy flux in the 0.1–100 GeV band ($G_{100}$), we can estimate the luminosity of the pulsar candidate, as done in the second Fermi-LAT pulsar catalog (hereafter referred to as 2PC; Abdo et al. 2013):

$$L_{\gamma,2PC} = 4\pi d^2 f_{\Omega} G_{100}, \quad (4)$$

where d is the distance of the source in kpc (estimated to 11 kpc) and $f_{\Omega}$ is the beam correction factor (Romani & Watters 2010). We assume $f_{\Omega} \approx 1$ as in 2PC, which corresponds to an outer magnetosphere fan-like beam sweeping the entire sky. The resulting luminosity is $L_{\gamma,2PC} = 9 \times 10^{35} (\frac{d}{11 \text{ kpc}})^2 \text{ erg s}^{-1}$.

The “heuristic” γ-ray pulsar luminosity is defined in the 2PC as

$$L_{\gamma}^{h} = \sqrt{10^{33} \dot{E}} \text{ erg s}^{-1}. \quad (5)$$

As the γ-ray luminosity cannot exceed the available power of the pulsar, the spin-down power cannot be lower than the estimated luminosity of ~9 × 10³⁵ erg s⁻¹. Following the maximum value of the spin-down power derived from the heuristic luminosity lower diagonal in Figure 7), we can say that the spin-down power of HESS J1641-463 is between ~10³⁵ and ~10³⁹ erg s⁻¹ for distances between 8 and 13 kpc (see diagonals in Figure 7).

As shown in Figure 7, the γ-ray luminosity is of the same order of magnitude as that of the Crab pulsar, while the spin-down power $\dot{E}$ could be either like the Crab or a few orders of magnitude lower ($L_{\gamma, Crab} = 1.96 \times 10^{35} \text{ erg s}^{-1}$ and $\dot{E}_{\text{Crab}} = 4.5 \times 10^{38} \text{ erg s}^{-1}$).

To test this pulsar hypothesis, we performed a search for γ-ray pulsations using 11 yr of Fermi-LAT data. This blind search used the two-stage approach described by Pletsch & Clark (2014), and employed in Clark et al. (2017), in which an efficient "semicoherent" search is performed to scan the entire parameter space for candidate signals, which are then followed up using a more sensitive fully coherent search. The lag window used for the semicoherent stage had width $\tau = 2^{22} \text{ s} \approx 49$ days. To increase the sensitivity of the search, we weighted the contribution of each photon to the pulsation TS (Kerr 2011). These weights represent the probability of each photon being emitted by the target source as opposed to being from a nearby source or from the diffuse background. The weights were calculated using `gtsrcprob`, using the best-fitting spectral model for the region of interest reported in Section 2.1. We searched for signals with frequency $\nu < 48$ Hz, spin-down rate of $\dot{\nu} > -10^{-10} \text{ Hz s}^{-1}$, and sky positions within a conservative $3\sigma$ range around the best-fitting position reported in Table 1. No significant candidate signals were detected. We note that the GeV spectrum of HESS J1641-463 falls off significantly more slowly at high energies than is typical for a γ-ray pulsar ($\Gamma_2 = 0.306$ versus $\Gamma_2 = 0.667$ assumed for pulsars in 4FGL). If a pulsar is present within this source, it is therefore likely to only contribute a fraction of the total emission, with additional contributions, especially at higher energies, from a putative PWN.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** γ-ray luminosity $L_{\gamma} = 4\pi f_{\Omega} d^2 G_{100}$ in the 0.1–100 GeV energy band vs. spin-down power $\dot{E}$ (taken from Abdo et al. 2013). The upper diagonal line indicates 100% conversion of spin-down power into γ-ray flux: for pulsars above this line, the distance *d* may be smaller, and/or the assumed beam correction $f_{\Omega} \equiv 1$ is wrong. The lower diagonal line indicates the heuristic luminosity $L_{\gamma}^{\text{he}} = \sqrt{10^{33} \dot{E}} \text{ erg s}^{-1}$. The black dashed lines are the dispersion values of the spin-down power of HESS J1641–463 based on the associated calculated luminosity for a distance of 8 (lower) and 13 kpc (upper).</figcaption></figure>

The lack of detected pulsations allows us to constrain this fraction, although this requires making several restrictive assumptions on the pulsar’s spin-down behavior that may not hold. Using the sensitivity estimation method described in Clark et al. (2017), we estimate a 95% confidence upper limit on the pulsed photon flux fraction from HESS J1641–463 of 42%. This estimate assumes that the putative pulsar spins down at a constant rate, and does not exhibit a significant long-term second frequency derivative, strong timing noise, or glitches, as these cause a significant loss of sensitivity for our search. Unfortunately, for a young pulsar of similar age to the SNR, these phenomena are likely to be present in the spin history over the 11 yr data set searched here. The maximum frequency and spin-down rates searched here were chosen with this in mind. The limiting values correspond to a maximum spin-down power of $\dot{E} \sim 2 \times 10^{38} \text{ erg s}^{-1}$. Assuming magnetic dipole braking (which has a braking index $n=3$ where $\dot{v} \propto v^n$), such an energetic pulsar would have an expected spin-frequency second derivative of $\dot{v} = n \dot{v}^2/v = 6 \times 10^{-22} \text{ Hz s}^{-2}$, while the semicoherent search method employed here loses sensitivity for $\dot{v} \gtrsim 10^{-22} \text{ Hz s}^{-2}$, making it unproductive to search for higher values of $v$ or $\dot{v}$.

### 4.2.3. PWN Scenario

The connection between the curved GeV and the hard TeV spectra remains unclear. This suggests that two mechanisms are producing the radiation in each energy band.

In the hypothesis where the GeV emission is produced by a pulsar, the TeV emission of HESS J1641–463 could be explained as its PWN emission via IC scattering.

First, we can use the integrated energy flux between 1 and 10 TeV from the H.E.S.S. Galactic Plane Survey (Abdalla et al. 2018a) to see if the putative pulsar could generate a PWN. The TeV luminosity of HESS J1641–463 obtained using Equation (4) is $L_{\gamma,\text{HGPS}} = 1.05 \times 10^{34} \text{ erg s}^{-1}$ for a distance of 11 kpc. This luminosity can easily be achieved using the pulsar spin-down power interval calculated in the pulsar scenario

described above. This is consistent with its pre-selection as a PWN candidate in Abdalla et al. (2018b).

However, the upper value of the $\dot{E}$ interval obtained in Figure 7, which could be in agreement with the lower boundary of the estimated age of the SNR (1.1–3.5 kyr), would be only valid for a young pulsar. In such a case, we would expect to detect an associated PWN at radio or X-ray wavelengths. Unfortunately, as reported in Sections 3.1 and 3.2, only upper limits on the X-ray and radio fluxes could be derived from this source. This is similar to the case of the relic PWN HESS J1303-631 (Beilicke et al. 2005) or HESS J1825-137 (Grondin et al. 2011). The ratio of the synchrotron and IC fluxes mainly depends on the magnetic field of the source assuming that both emission come from the same region (one-zone model). Therefore, we use the formula in Aharonian et al. (2009) where the ratio between the γ-ray flux between 1 and 10 TeV and the X-ray upper limit between 2 and 6 keV is used to estimate an upper limit on the magnetic field of the source of ~4.3 μG. This magnetic field value is similar to relic PWNe and would mean that the associated pulsar would be older than 3.5 kyr.

## 5. Conclusion

With more than 8 yr of Fermi-LAT data, our morphological and spectral analyses allow us to bring further constraints on the origin of the γ-ray emission of the neighboring sources HESS J1640–465 and HESS J1641–463. The morphological analysis of both sources above 1 GeV reveal a significant extension for HESS J1640–465 in the GeV energy band, together with a hard spectrum following a power law with a slope of $\Gamma = 1.8 \pm 0.1_{\text{stat}} \pm 0.2_{\text{syst}}$, which connects well with the H.E.S.S. spectrum. When compared with the analysis of Lemoine-Goumard et al. (2014), this analysis shows the strong impact of the confusion of sources on the spectral index result (see Section 2.3). Two origins of the multiwavelength emission are investigated: a PWN origin associated with the X-ray pulsar PSR J1640–4631 correlated with the radio SNR G338.3-0.0 and a scenario in which protons or electrons are accelerated in the SNR shock. Both hypotheses can be used to reproduce the detected emission, though the scenario involving electrons accelerated at the SNR shock requires a high electron-to-proton ratio. The H.E.S.S. PeVatron’s candidate HESS J1641–463 is detected as a point-like source and its GeV spectrum is similar to those seen from pulsars detected by Fermi-LAT. Upper limits at X-ray and radio wavelengths were also derived using Chandra X-ray and SGPS radio observations. The new data pose severe constraints on the scenarios investigated: assuming that a pulsar produces the GeV emission detected by the LAT, constraints on its spin-down power are obtained, as well as an upper limit on the magnetic field of its putative PWN detected at TeV energies. Alternative possibilities are that protons are accelerated up to hundreds of TeVs by direct acceleration of CR protons by the SNR G338.5+0.1 (HESS J1641–463) or by reacceleration of ambient CRs cannot be discarded. However, the very large ratio between the radio upper limit and the γ-ray emission casts doubts on the viability of such scenarios. Finally, the hard TeV spectrum detected by H.E.S.S. can be reproduced by the illumination of ambient dense gas by CRs that escaped from the SNR G338.3-0.0 coincident with HESS J1640–465. Future radio (pulsation search) and γ-ray observations (morphology and spectral curvature) by CTA are essential
to disentangle the nature of these two intriguing sources and VHE emitters.

The Fermi-LAT Collaboration acknowledges support for LAT development, operation, and data analysis from NASA and DOE (United States), CEA/Irfu and IN2P3/CNRS (France), ASI and INFN (Italy), MEXT, KEK, and JAXA (Japan), and the K.A. Wallenberg Foundation, the Swedish Research Council and the National Space Board (Sweden). Scientific analysis support in the operations phase from INAF (Italy) and CNES (France) is also gratefully acknowledged. This work was performed in part under DOE Contract DE-AC02-76SF00515.

Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Etudes Spatiales in France.

M.H.G., M.L.G., and S.G. acknowledge support from Agence Nationale de la Recherche (grant ANR-17-CE31-0014).

The contribution of J.D.G. is supported by the National Aeronautics and Space Administration (NASA) under grant No. NNX17AL74G issued through the NNH16ZDA001N Astrophysics Data Analysis Program (ADAP). J.D.G. is also supported by the NYU Abu Dhabi Research Enhancement Fund (REF) under grant RE022 and NYU Abu Dhabi Grant AD022.

S.G. acknowledges support from the Observatory of Paris (Action Fédératrice CTA).

We thank David A. Smith for useful discussions.

C.J.C. acknowledges support from the ERC under the European Union's Horizon 2020 research and innovation program (grant agreement No. 715051; Spiders).

This research has made use of NASA's Astrophysics Data System Bibliographic Services.

## Appendix

### Pulsar Wind Nebula Modeling of HESS J1640-465

As mentioned in Section 4.1.1, the best way to test a PWN origin for the γ-rays observed from HESS J1640-465 is to determine if the dynamical and spectral properties of this source can be reproduced by a model for the evolution of a PWN inside an SNR. We use a model similar to that described by Gelfand et al. (2009), which successfully reproduced the properties of other γ-ray emitting PWNe such as G54.1+0.3 (Gelfand et al. 2015) and Kes 75 (Gelfand et al. 2014). In this “one-zone” model, the magnetic field strength and particle distribution inside the PWN are assumed to be uniform (e.g., Reynolds & Chevalier 1984; Gelfand et al. 2009), and the SNR is expanding inside a medium with uniform number density $n_{ism}$. Furthermore, it assumes the density profile of the unshocked supernova ejecta consists of a uniform density ($\rho \propto r^0$) core surrounded by a $\rho \propto r^{-9}$ envelope (e.g., Blondin et al. 2001; Gelfand et al. 2009), where r is the distance from the explosion site (assumed to be the center of the SNR). While different supernova progenitors will likely result in ejecta with different density profiles, this does not substantially affect the evolution of the PWN (Chevalier 2005). In some cases, it is possible to assume the dependence between the ISM density $n_{ism}$, initial kinetic energy $E_{sn}$ and mass of the ejecta $M_{ej}$ as predicted by self-similar evolution of the SNR. However, this is not possible for HESS J1640-465 for two reasons:

1. We are attempting to reproduce the size of both the PWN and the SNR. Before the PWN collides with the SNR's reverse shock, its size is sensitive to the density of the innermost ejecta, whose dependence on $E_{sn}$ and $M_{ej}$ differs from the dependence of the SNR's forward shock on these parameters (and $n_{ism}$). Therefore, two different parameterizations are needed, which increases the number of model parameters.

2. The large size of the PWN relative to the SNR shell in this system ($\theta_{pwn} \sim 0.75\theta_{snr}$; Table 6) suggests that it might have already collided with the SNR reverse shock. In this case, the evolution of the PWN is sensitive to a number of different parameters (e.g., the time of the collision and the pressure of the reverse-shocked material; see Gelfand et al. 2009). While our model suggests this collision has not happened quite yet, there was no basis for rejecting this possibility a priori.

In this model, the pulsar injects energy into the PWN at a rate of $\dot{E}$, which evolves as (e.g., Goldreich & Julian 1969; Pacini & Salvati 1973)

$$ \dot{E}(t) = \dot{E}_0 \left(1 + \frac{t}{\tau_{sd}}\right)^{-\frac{p+1}{p-1}}, \quad (A1) $$

where $\dot{E}_0$ is the initial spin-down luminosity, p is the braking index, and $\tau_{sd}$ is the spin-down timescale (equivalent to the “characteristic age” at birth) of the associated pulsar PSR J1640-4631. As listed in Table 6, timing observations of this neutron star have determined its current braking index, characteristic age $t_{ch}$, and spin-down luminosity. Therefore, for a particular value of $\tau_{sd}$, we can determine the “true” age $t_{age}$,

$$ t_{age} = \frac{2t_{ch}}{p-1} - \tau_{sd} \quad (A2) $$

of this system as well as $\dot{E}_0$,

$$ \dot{E}_0 = \dot{E} \left( 1 + \frac{t_{age}}{\tau_{sd}} \right)^{\frac{p+1}{p-1}}, \quad (A3) $$

assuming that p is constant with time. This allows us to decrease the number of model parameters by two.

Furthermore, this model assumes that a constant fraction $\eta_B$ of this energy is injected into the PWN in the form of magnetic fields, while the remainder $(1 - \eta_B)$ is injected as relativistic $e^{\pm}$ whose spectrum is described by a broken power law

$$ \frac{d\dot{N}}{dE} = \begin{cases} \dot{N}_{\text{break}} \left(\frac{E}{E_{\text{break}}}\right)^{-p_1} & E_{\min} < E < E_{\text{break}} \\ \dot{N}_{\text{break}} \left(\frac{E}{E_{\text{break}}}\right)^{-p_2} & E_{\text{break}} < E < E_{\max} \end{cases}, \quad (A4) $$

where the five free parameters ($E_{min}$, $E_{break}$, $E_{max}$, $p_1$, and $p_2$) in Equation (A4) are assumed to be constant with time and the normalization $\dot{N}_{\text{break}}$ is calculated by requiring that

$$ (1 - \eta_B) \dot{E} = \int_{E_{\min}}^{E_{\max}} E \frac{d\dot{N}}{dE} dE \quad (A5) $$

at all times $t$. As the PWN evolves, the energy of these particles changes due to adiabatic expansion/contraction of the nebulae
<table><caption>**Table 6**

Observed Properties of HESS J1640-465, as Well as the Value Predicted by Our Model for the Combination of Input Parameters Given in Table 4</caption><thead><tr><th>Property</th><th>Observed Value</th><th>Model Prediction</th><th>Citation</th></tr></thead><tbody><tr><td>Current spin-down luminosity <i>˙E</i></td><td>4.4 × 10<sup>36</sup> erg s<sup>-1</sup></td><td>...</td><td>Gotthelf et al. (2014)</td></tr><tr><td>Characteristic age <i>t</i><sub>ch</sub></td><td>3350 yr</td><td>...</td><td>Gotthelf et al. (2014)</td></tr><tr><td>Braking index <i>p</i></td><td>3.15 ± 0.03</td><td>≡ 3.15</td><td>Archibald et al. (2016)</td></tr><tr><td>Radius θ<sub>snr</sub></td><td>4.45 ± 0.5</td><td>4.2</td><td>Shaver & Goss (1970)</td></tr><tr><td>Radius θ<sub>pwn</sub></td><td>3.3 ± 0.2</td><td>3.5</td><td>Lemiere et al. (2009)</td></tr><tr><td>660 MHz flux density <i>S</i><sub>660</sub></td><td>&lt; 690 ± 300 mJy</td><td>70 mJy</td><td>Castelletti et al. (2011)</td></tr><tr><td>2–25 keV unabsorbed flux</td><td>(1.7 ± 0.4) × 10<sup>-12</sup>(erg cm<sup>-2</sup> s<sup>-1</sup>)<sup>*</sup></td><td>1.67 × 10<sup>-12</sup> (erg cm<sup>-2</sup> s<sup>-1</sup>)<sup>*</sup></td><td>Gotthelf et al. (2014)</td></tr><tr><td>2–25 keV photon index</td><td>2.2<sup>+0.7</sup><sub>-0.4</sub></td><td>2.58</td><td>Gotthelf et al. (2014)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(0.35 GeV)</td><td>(1.32 ± 1.57) × 10<sup>-12</sup>*</td><td>0.76 × 10<sup>-13</sup>*</td><td>...</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(0.785 GeV)</td><td>1.38<sup>+2.19</sup><sub>-2.24</sub> × 10<sup>-12</sup>*</td><td>1.27 × 10<sup>-12</sup>*</td><td>...</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(1.68 GeV)</td><td>3.22<sup>+1.46</sup><sub>-1.48</sub> × 10<sup>-12</sup>*</td><td>1.99 × 10<sup>-12</sup>*</td><td>...</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(3.59 GeV)</td><td>3.64<sup>+1.02</sup><sub>-1.02</sub> × 10<sup>-12</sup>*</td><td>2.96 × 10<sup>-12</sup>*</td><td>...</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(7.67 GeV)</td><td>3.01<sup>+1.02</sup><sub>-1.02</sub> × 10<sup>-12</sup>*</td><td>4.15 × 10<sup>-12</sup>*</td><td>...</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(16.4 GeV)</td><td>3.54<sup>+1.03</sup><sub>-1.03</sub> × 10<sup>-12</sup>*</td><td>5.49 × 10<sup>-12</sup>*</td><td>...</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(35 GeV)</td><td>6.31<sup>+1.61</sup><sub>-1.61</sub> × 10<sup>-12</sup>*</td><td>6.77 × 10<sup>-12</sup>*</td><td>...</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(75 GeV)</td><td>8.81<sup>+2.75</sup><sub>-2.75</sub> × 10<sup>-12</sup>*</td><td>6.98 × 10<sup>-12</sup>*</td><td>...</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(160 GeV)</td><td>5.65<sup>+2.71</sup><sub>-2.71</sub> × 10<sup>-12</sup>*</td><td>7.15 × 10<sup>-12</sup>*</td><td>...</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(289 GeV)</td><td>5.34<sup>+2.27</sup><sub>-2.42</sub> × 10<sup>-12</sup>*</td><td>7.05 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(340 GeV)</td><td>&lt; 6.47 × 10<sup>-12</sup>*</td><td>6.38 × 10<sup>-12</sup>*</td><td>...</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(350 GeV)</td><td>5.92<sup>+1.01</sup><sub>-1.04</sub> × 10<sup>-12</sup>*</td><td>6.50 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(424 GeV)</td><td>6.99<sup>+0.78</sup><sub>-0.80</sub> × 10<sup>-12</sup>*</td><td>6.10 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(513 GeV)</td><td>5.57<sup>+0.60</sup><sub>-0.62</sub> × 10<sup>-12</sup>*</td><td>5.52 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(622 GeV)</td><td>5.66<sup>+0.57</sup><sub>-0.59</sub> × 10<sup>-12</sup>*</td><td>5.95 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(734 GeV)</td><td>5.08<sup>+6.36</sup><sub>-6.36</sub> × 10<sup>-12</sup>*</td><td>5.20 × 10<sup>-12</sup>*</td><td>...</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(753 GeV)</td><td>6.37<sup>+0.55</sup><sub>-0.57</sub> × 10<sup>-12</sup>*</td><td>5.34 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(913 GeV)</td><td>5.60<sup>+0.50</sup><sub>-0.52</sub> × 10<sup>-12</sup>*</td><td>5.04 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(1.11 TeV)</td><td>5.22<sup>+0.49</sup><sub>-0.50</sub> × 10<sup>-12</sup>*</td><td>4.38 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(1.34 TeV)</td><td>5.57<sup>+0.49</sup><sub>-0.51</sub> × 10<sup>-12</sup>*</td><td>4.20 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(1.62 TeV)</td><td>4.59<sup>+0.47</sup><sub>-0.49</sub> × 10<sup>-12</sup>*</td><td>4.28 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(1.97 TeV)</td><td>3.99<sup>+0.46</sup><sub>-0.48</sub> × 10<sup>-12</sup>*</td><td>3.91 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(2.38 TeV)</td><td>2.87<sup>+0.43</sup><sub>-0.45</sub> × 10<sup>-12</sup>*</td><td>3.44 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(2.89 TeV)</td><td>3.77<sup>+0.48</sup><sub>-0.50</sub> × 10<sup>-12</sup>*</td><td>3.20 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(3.50 GeV)</td><td>3.24<sup>+0.48</sup><sub>-0.50</sub> × 10<sup>-12</sup>*</td><td>2.73 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(4.24 TeV)</td><td>2.65<sup>+0.48</sup><sub>-0.51</sub> × 10<sup>-12</sup>*</td><td>2.56 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(5.13 TeV)</td><td>1.96<sup>+0.45</sup><sub>-0.48</sub> × 10<sup>-12</sup>*</td><td>2.63 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(6.22 TeV)</td><td>1.70<sup>+0.45</sup><sub>-0.49</sub> × 10<sup>-12</sup>*</td><td>2.14 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(7.53 TeV)</td><td>1.95<sup>+0.50</sup><sub>-0.55</sub> × 10<sup>-12</sup>*</td><td>2.04 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr><tr><td>ν<i>F</i><sub>ν</sub>(10.2 TeV)</td><td>8.77<sup>+3.21</sup><sub>-3.53</sub> × 10<sup>-13</sup>*</td><td>1.55 × 10<sup>-12</sup>*</td><td>Abramowski et al. (2014a)</td></tr></tbody></table>

Note. The errors for the values of $\nu F_{\nu}$ in the $\gamma$-ray regime are the square root of the sum of the squares of the statistical and systematic errors. The ... indicate quantities whose values were fixed during the modeling and the * indicates that the units are in erg cm$^{-2}$ s$^{-1}$.

(e.g., Gelfand et al. 2009)

$$ \dot{E}_{\text{ad}}(t, E) = -\frac{\dot{R}_{\text{pwn}}(t)}{R_{\text{pwn}}(t)}, \qquad (\text{A6}) $$

where $R_{\text{pwn}}$ is the radius of the (assumed spherical) PWN, calculated using the prescription described by Gelfand et al. (2009), and radiative losses resulting from the synchrotron and IC emission of this particles. The synchrotron emission is calculated assuming that, inside the PWN, the angle between a particle's velocity and local magnetic field is randomly distributed, and the evolution of the PWN's magnetic field is calculated assuming magnetic flux conservation as described

by Gelfand et al. (2009). For IC emission, we assume the particles are scattering photons from the cosmic microwave background (temperature $T_{\text{cmb}} = 2.7255$ K; Fixsen 2009) as well as additional background fields with a blackbody spectrum and normalization $K_{\text{IC}}$, such that this photon field has an energy density of

$$ u_{\text{ic}} = K_{\text{ic}} a_{\text{bb}} T_{\text{ic}}^{4}, \qquad (\text{A7}) $$

where $a_{\text{bb}} \approx 7.5657 \times 10^{-15}$ erg cm$^{-3}$ K$^{-4}$.

Similar to the process described by Gelfand et al. (2015), we used a Metropolis Markoff Chain Monte Carlo (e.g., Metropolis 1985) to determine the combination of the 15 model input parameters $\Theta$, listed in Table 4, the best reproduces the 35 observed properties $D$ listed in Table 6. This was done by
identifying the combination $\theta$ whose predicted observable quantities $M$ resulted in the lowest $\chi^2$, defined as

$$ \chi^2 = \sum_{i=1}^{i=35} \left( \frac{D_i - M_i}{\sigma_i} \right)^2, \qquad (A8) $$

where $\sigma_i$ is the error on observed quantity $D_i$ listed in Table 6. While such a procedure can be used to identify degeneracies between the model parameters as well as estimate their errors, the necessary mapping of the possible parameter space requires $\gtrsim 10^6$ trials (e.g., Gelfand et al. 2015), which is beyond the scope of this work.

## ORCID iDs

A. Mares https://orcid.org/0000-0003-3345-7800

M. Lemoine-Goumad https://orcid.org/0000-0002-4462-3686

C. J. Clark https://orcid.org/0000-0003-4355-3572

J. D. Gelfand https://orcid.org/0000-0003-4679-1058

D. A. Green https://orcid.org/0000-0003-3189-9998

M.-H. Grondin https://orcid.org/0000-0002-8383-251X

## References

Abdalla, H., Abramowski, A., Aharonian, F., et al. 2018a, A&A, 612, A1

Abdalla, H., Abramowski, A., Aharonian, F., et al. 2018b, A&A, 612, A2

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, PhRvD, 80, 122004

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011, ApJ, 734, 28

Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33

Abramowski, A., Aharonian, F., Benkhali, F. A., et al. 2014a, MNRAS, 439, 2828

Abramowski, A., Aharonian, F., Benkhali, F. A., et al. 2014b, ApJL, 794, L1

Acero, F., Ackermann, M., Ajello, M., et al. 2016, ApJS, 224, 8

Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, ApJ, 636, 777

Aharonian, F., Akhperjanian, A. G., de Almeida, U. B., et al. 2009, ApJ, 692, 1500

Ajello, M., Angioni, R., Axelsson, M., et al. 2020, ApJ, 892, 105

Archibald, R. F., Gotthelf, E. V., Ferdman, R. D., et al. 2016, ApJL, 819, L16

Atwood, W., Albert, A., Baldini, L., et al. 2013, arXiv:1303.3514

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071

Baade, W., & Zwicky, F. 1934, PhRv, 46, 76

Beilicke, M., de Naurois, M., Raue, M., Rowell, G., & Schlenker, S. 2005, in Proc of the 22nd Texas Symp. on Relativistic Astrophysics, ed. P. Chen et al. (Stanford, CA: Stanford Univ.), 654

Blondin, J. M., Chevalier, R. A., & Frierson, D. M. 2001, ApJ, 563, 806

Bruel, P., Burnett, T. H., Digel, S. W., et al. 2018, arXiv:1810.11394

Bucciantini, N., Arons, J., & Amato, E. 2011, MNRAS, 410, 381

Castelletti, G., Giacani, E., Dubner, G., et al. 2011, A&A, 536, A98

Celli, S., Morlino, G., Gabici, S., & Aharonian, F. A. 2019, MNRAS, 487, 3199

Chevalier, R. A. 2005, ApJ, 619, 839

Clark, C. J., Wu, J., Pletsch, H. J., et al. 2017, ApJ, 834, 106

Cristofari, P., Gabici, S., Casanova, S., Terrier, R., & Parizot, E. 2013, MNRAS, 434, 2748

Davies, B., de La Fuente, D., Najarro, F., et al. 2012, MNRAS, 419, 1860

de Wilt, P., Rowell, G., Walsh, A. J., et al. 2017, MNRAS, 468, 2093

Fazio, G. G., Helmken, H. F., O'Mongain, E., & Weekes, T. C. 1972, ApJL, 175, L117

Fixsen, D. J. 2009, ApJ, 707, 916

Funk, S., Hinton, J. A., Pühlhofer, G., et al. 2007, ApJ, 662, 517

Gabici, S., & Aharonian, F. A. 2014, MNRAS, 445, L70

Gabici, S., Aharonian, F. A., & Casanova, S. 2009, MNRAS, 396, 1629

Gaensler, B. M., & Slane, P. O. 2006, ARA&A, 44, 17

Gelfand, J. D. 2017, in Radiative Models of Pulsar Wind Nebulae, ed. D. F. Torres (Berlin: Springer International), 161

Gelfand, J. D., Slane, P. O., & Temim, T. 2014, AN, 335, 318

Gelfand, J. D., Slane, P. O., & Temim, T. 2015, ApJ, 807, 30

Gelfand, J. D., Slane, P. O., & Zhang, W. 2009, ApJ, 703, 2051

Goldreich, P., & Julian, W. H. 1969, ApJ, 157, 869

Gotthelf, E. V., Tomsick, J. A., Halpern, J. P., et al. 2014, ApJ, 788, 155

Grondin, M.-H., Funk, S., Lemoine-Goumad, M., et al. 2011, ApJ, 738, 42

Haverkorn, M., Gaensler, B. M., McClure-Griffiths, N. M., Dickey, J. M., & Green, A. J. 2006, ApJS, 167, 230

Inoue, T., Yamazaki, R., Inutsuka, S.-i., & Fukui, Y. 2012, ApJ, 744, 71

Johnston, S., Smith, D. A., Karastergiou, A., & Kramer, M. 2020, MNRAS, 497, 1957

Katz, B., & Waxman, E. 2008, JCAP, 01, 018

Kennel, C. F., & Coroniti, F. V. 1984, ApJ, 283, 694

Kerr, M. 2011, ApJ, 732, 38

Klepser, S., Carrigan, S., de Oña Wilhelmi, E., et al. 2013, in 33rd Int. Cosmic Ray Conf. (ICRC2013) (Trieste: PoS), 971

Klothes, R., & Dougherty, S. M. 2007, A&A, 468, 993

Lau, J. C., Rowell, G., Burton, M. G., et al. 2017, MNRAS, 464, 3757

Lemiere, A., Slane, P., Gaensler, B. M., & Murray, S. 2009, ApJ, 706, 1269

Lemoine-Goumad, M., Grondin, M.-H., Acero, F., et al. 2014, ApJL, 794, L16

Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396

Metropolis, N. 1985, in Monte Carlo: In the Beginning and Some Great Expectations, ed. R. Alcouffe et al., Vol. 240 (Berlin: Springer-Verlag), 62

Moskalenko, I., Johannesson, G., Orlando, E., Porter, T., & Strong, A. W. 2017, in 35th Int. Cosmic Ray Conf. (ICRC2017), 35 (Trieste: PoS), 279

Nishikawa, K. I., Hardee, P., Richardson, G., et al. 2005, ApJ, 622, 927

Pacini, F., & Salvati, M. 1973, ApJ, 186, 249

Pedaletti, G., de Oña Wilhelmi, E., Torres, D. F., & Natale, G. 2015, JHEAp, 5, 15

Pletsch, H. J., & Clark, C. J. 2014, ApJ, 795, 75

Porter, T. A., Moskalenko, I. V., & Strong, A. W. 2006, ApJL, 648, L29

Reynolds, S. P., & Chevalier, R. A. 1984, ApJ, 278, 630

Romani, R. W., & Watters, K. P. 2010, ApJ, 714, 810

Shaver, P. A., & Goss, W. M. 1970, AuJPA, 14, 77

Slane, P., Castro, D., Funk, S., et al. 2010, ApJ, 720, 266

Smith, D. A., Bruel, P., Cognard, I., et al. 2019, ApJ, 871, 78

Spitkovsky, A. 2008, ApJL, 682, L5

Supan, L., Supanitsky, A. D., & Castelletti, G. 2016, A&A, 589, A51

Tanaka, S. J., & Takahara, F. 2011, ApJ, 741, 40

Tang, Y., Yang, C., Zhang, L., & Wang, J. 2015, ApJ, 812, 32

Torres, D. F., Cillis, A., Martín, J., & de Oña Wilhelmi, E. 2014, JHEAp, 1, 31

Uchiyama, Y., Blandford, R. D., Funk, S., Tajima, H., & Tanaka, T. 2010, ApJL, 723, L122

Whiteoak, J. B. Z., & Green, A. J. 1996, A&AS, 118, 329

Wood, M., Caputo, R., Charles, E., et al. 2017, in 35th Int. Cosmis Ray Conf. (ICRC2017) (Trieste: PoS), 824

Xin, Y.-L., Liao, N.-H., Guo, X.-L., et al. 2018, ApJ, 867, 55

Zhu, B.-T., Zhang, L., & Fang, J. 2018, A&A, 609, A110

Zirakashvili, V. N., & Aharonian, F. A. 2010, ApJ, 708, 965