The Eel Pulsar Wind Nebula: a PeVatron-Candidate Origin for HAWC J1826-128 and HESS J1826-130

HAWC J1826-128 is one of the brightest Galactic TeV gamma-ray sources detected by the High Altitude Water Cherenkov (HAWC) Observatory, with photon energies extending up to nearly $\sim$100 TeV. This HAWC source spatially coincides with the H.E.S.S. TeV source HESS J1826-130 and the"Eel"pulsar wind nebula (PWN), which is associated with the GeV pulsar PSR J1826-1256. In the X-ray band, Chandra and XMM-Newton revealed that the Eel PWN is composed of both a compact nebula ($\sim$15") and diffuse X-ray emission ($\sim$6'$\times$2') extending away from the pulsar. Our NuSTAR observation detected hard X-ray emission from the compact PWN up to $\sim$20 keV and evidence of the synchrotron burn-off effect. In addition to the spatial coincidence between HESS J1826-130 and the diffuse X-ray PWN, our multi-wavelength spectral energy distribution (SED) analysis using X-ray and gamma-ray data establishes a leptonic origin of the TeV emission associated with the Eel PWN. Furthermore, our evolutionary PWN SED model suggests (1) a low PWN B-field of $\sim$1 $\mu$G, (2) a significantly younger pulsar age ($t \sim5.7$ kyr) than the characteristic age ($\tau= 14.4$ kyr) and (3) a maximum electron energy of $E_{max} = 2$ PeV. The low B-field as well as the putative supersonic motion of the pulsar may account for the asymmetric morphology of the diffuse X-ray emission. Our results suggest that the Eel PWN may be a leptonic PeVatron particle accelerator powered by the $\sim$6-kyr-old pulsar PSR J1826-1256 with a spin-down power of $3.6 \times 10^{36}$ erg s$^{-1}$.


INTRODUCTION
The origin of extremely energetic cosmic rays (CRs) above ∼1 PeV is a long-standing problem in high energy astrophysics. Since the directional information of each CR is washed out by interstellar magnetic fields, the energetic CR acceleration in our Galaxy can be probed indirectly through their candidate accelerators, which are frequently associated with TeV γ-ray sources.
The High-Altitude Water Cherenkov (HAWC) experiment has opened a new window to the high energy γray sky. Equipped with 300 water Cherenkov detectors (WCDs), HAWC directly detects air shower particles produced by very high energy (VHE, E ≥ 100 GeV) up through ultra high energy (UHE, E ≥ 100 TeV) γ-rays in the upper atmosphere (Smith & HAWC Collaboration 2015). The HAWC observatory is distinct from imaging air Cherenkov telescopes such as H.E.S.S., VERITAS and MAGIC due to its sensitivity at energies above ∼10 TeV, a result of employing WCDs and collecting data from γ-ray sources continuously under all weather conditions. Alongside the LHAASO (Large High Altitude Air Shower Observatory) experiment, HAWC explores astrophysical sources in the highest energy γ-ray band beyond ∼100 TeV. A handful of HAWC sources detected at and beyond ∼100 TeV are of particular interest as they may be produced by one of the most extreme classes of cosmic particle accelerator -a Galactic PeVatron (Abeysekara et al. 2020). Primary Galactic PeVatron candi-dates that may accelerate electrons and protons to PeV energies are pulsar wind nebulae (Arons 2012), γ-ray superbubbles (e.g. the Cygnus cocoon; Abeysekara et al. 2021) and the supermassive black hole Sgr A* (HESS Collaboration et al. 2016).
While the number of Galactic PeVatron candidates detected by HAWC is growing, it remains uncertain whether leptonic or hadronic accelerators are primarily powering their UHE TeV emission. We now have strong indications that CR electrons and positrons are generated by leptonic processes in pulsar wind nebulae (PWNe), are able to exit these PWNe, and may diffuse to great distances (Cholis et al. 2018), in part due to the recent discovery of TeV γ-ray halos around the Geminga and Monogem pulsars (Abeysekara et al. 2017;Linden et al. 2017). The existence of γ-ray halos and successful leptonic model fits to UHE HAWC sources (e.g. Sudoh et al. 2021) both support a leptonic PWN hypothesis for these UHE HAWC sources and, by association, for Galactic PeVatrons and the origin of CRs.
PWNe are composed of highly relativistic particle winds that are injected by energetic rotation-powered pulsars into the surrounding SNR ejecta or interstellar medium. High energy observations have detected synchrotron and inverse Compton scattering (ICS) emission from numerous PWNe in the X-ray and TeV bands respectively, suggesting that non-thermal electrons are accelerated to GeV-PeV energies within PWNe. Given the complex emission mechanisms and morphology of TeV sources, the origins of TeV emission are best explored through multi-wavelength observations from the radio, X-ray, GeV, and up to TeV energy bands. For example, PWNe are often found in the vicinity of supernova remnants (SNRs) and molecular clouds, and therefore can be difficult to spatially resolve with γ-ray observations alone.
Over the past two decades, X-ray and TeV γ-ray observations have bolstered our understanding of PWN astrophysics and electron acceleration in the Galaxy. In a typical PWN, synchrotron radiation dominates from radio to keV/MeV photon energies, while the TeV emission is produced by ICS of accelerated electrons off CMB photons or local IR/optical/UV radiation fields from nearby dust, molecular clouds and star clusters (e.g. Slane 2017). As of today, more than 70 X-ray emitting PWNe have been discovered and nearly half of them have been detected in the TeV band (Kargaltsev et al. 2013). Chandra observations of PWNe have revealed detailed X-ray features such as tori, jets, and bow-shocks (Kargaltsev et al. 2015). On the other hand, NuSTAR, operating in its broad bandwidth of 3 to 79 keV, has detected hard X-ray synchrotron emission from some of the highest energy electron populations in our Galaxy (E e > 100 TeV), including PWNe and X-ray filaments in the Galactic Center (Zhang et al. 2014;Mori et al. 2015;Reynolds et al. 2017).
The highly extended UHE source eHWC J1825−134 belongs to a small group of nine γ-ray sources detected above 56 TeV, registered in the energetic HAWC (eHWC) source catalog (Abeysekara et al. 2020). This eHWC detection corresponds with LHAASO J1825−1326, a new UHE γ-ray detection over 1 • in extension with a photon E max of 0.42 ± 0.16 PeV (Cao et al. 2021). A more recent HAWC analysis was able to resolve eHWC J1825−134 into three separate high energy TeV sources, one of which is HAWC J1826−128 (Albert et al. 2021). This HAWC source is spatially coincident with the known TeV source HESS J1826−130, which may imply that HAWC J1826−128 is a separate detection of this TeV accelerator, providing a higher energy component to the original H.E.S.S. source. Albert et al. (2021) found that the HAWC J1826−128 spectrum is hard in TeV γ-rays, reaching photon energies of nearly 100 TeV. HESS J1826−130 also features a relatively hard spectrum (H. E. S. S. Collaboration et al. 2020) and the similarly hard spectral indices of both sources (Γ < 2) further supports their association. At larger transverse distances, HAWC J1826−128 and HESS J1826−130 are adjacent to the extended and softer TeV emission from HAWC J1825−138 in the south (the likely counterpart to HESS J1825−137) and HAWC J1825−134, a newly-discovered γ-ray source (detected above ∼200 TeV) likely powered by hadronic particle acceleration (see Figure 1 and Albert et al. 2021).
A GeV pulsar (PSR J1826−1256) of high spin-down luminosityĖ and its diffuse PWN are roughly coincident with HAWC J1826−128 (pulsar 5 .4 offset from γ-ray centroid) and HESS J1826−130 (pulsar 4 .6 offset from the γ-ray centroid). PSR J1826−1256 is one of the brightest GeV γ-ray pulsars in the Fermi 4-year catalog 3FGL (Acero et al. 2015), and has been categorized as "Vela-like" (Karpova et al. 2019). Located at R.A.(J2000) = 18 h 26 m 08.5 s decl.(J2000) = -12 • 56 33.36 , PSR J1826−1256 has a characteristic age of 14.4 kyr, spin-down luminosityĖ of 3.6 × 10 36 erg s −1 , and surface magnetic field B = 3.7 × 10 12 G . XMM-Newton and Chandra X-ray observations have characterized PSR J1826−1256 and the surrounding PWN (known as the Eel PWN; Roberts et al. 2007) in the X-ray band below 10 keV. Previous analyses of these XMM-Newton and Chandra observations and the various X-ray and TeV γ-ray properties of the Eel PWN include Karpova et al. (2019) and Duvidovich et al. (2019). In the X-ray band, the Eel PWN is char-  acterised by a tail-like extended emission over 6 (see Figure 2) and a more compact nebula within r ∼ 15 from the pulsar (Karpova et al. 2019). Hereafter we refer to them as "diffuse PWN" and "compact PWN", respectively. Additionally, Karpova et al. (2019) used the relation between distance and interstellar reddening to derive a new PSR J1826−1256 distance estimate of 3.5 kpc, which we adopt throughout this paper for all transverse distance calculations (e.g. 1 ∼ 1D 3.5kpc pc). Motivated by a potential association between the HAWC and H.E.S.S. sources and the Eel PWN, we performed a NuSTAR observation of the PWN in order to obtain broad-band X-ray data and combine it with the existing XMM data, as well as HAWC/H.E.S.S. TeV data. Our primary goal is to explore the hypothesis that the Eel PWN may be the PeVatron candidate HAWC J1826−128/HESS J1826−130 using multiwavelength spectral energy distribution (SED) fitting and an analysis of the PWN characteristics and morphology in the X-ray band. We begin by reviewing previous observations of the HAWC and H.E.S.S. sources ( §2), the status of Fermi-LAT detection of the pulsar and HESS J1826−130 ( §3), the X-ray observations performed by NuSTAR and XMM-Newton ( §4), and a 90 cm (330 MHz) radio image of the Eel PWN region ( §5). In §6 we perform analysis on the NuSTAR X-ray images, and in §7 we perform spectral analyses across the three X-ray telescopes. In §8 we investigate scenarios that may describe the observed emission from HAWC J1826−128 and HESS J1826−130 through SED fitting, and in §9 we discuss what is powering HAWC J1826−128 and HESS J1826−130, the implications of our fitting results, our conclusions about the Eel PWN morphology, and areas that require further study.  The centroids of HAWC J1826−128 and HESS J1826−130 are offset by approximately 10 , and both source extensions overlap with the diffuse Eel PWN (see Figure 1). Because of this spatial coincidence, along with the spectral consistency between the two sources within uncertainties, we proceed under the assumption that these are two separate detections of the same γ-ray source, as suggested by Albert et al. (2021).

Fermi-LAT OBSERVATIONS AND DATA ANALYSIS
The pulsar associated with the Eel PWN (PSR J1826−1256) is registered in both the 3FHL (Ajello et al. 2017) and 4FGL (Abdollahi et al. 2020) catalogs, as 3FHL/4FGL J1826.1−1256. To investigate whether Fermi-LAT detected any counts from the location of HESS J1826−130 that could not be attributed to the pulsar, we analyzed 13 years of Fermi Large Area Telescope (LAT) data collected between August 4, 2008, and August 3, 2021 (Fermi mission elapsed time 239557417s − 649724483s). We used the latest Pass 8 version of the data with "Source" class (evclass=128) "FRONT+BACK" type (evtype=3) events in the energy range from 100 MeV to 500 GeV. We used the instrument response function P8R3 SOURCE V3 and the filter expression (DATA QUAL>0)&&(LAT CONFIG==1).
To limit the γ-ray contamination from Earth's limb, the maximum zenith angle was set to 90 • . The region of interest (ROI) is a square area of 10 • ×10 • centered at the position of 4FGL J1826.1−1256. All point-like and extended sources within a square area of 30 • ×30 • from the center of the ROI listed in LAT 10-year Source Catalog (4FGL-DR2) were included in the model. The Galactic diffuse emission and extragalactic isotropic emission were modeled using gll iem v07.fits and iso P8R3 SOURCE V3 v1.txt, respectively 1 . We undertook the binned likelihood analysis with Fermipy  PSR J1826−1256 is a young, energetic pulsar with emission detected above 20 GeV (Ajello et al. 2017). Therefore, we generated test statistics (TS) maps above 30 GeV, 40 GeV, and 50 GeV to detect any emission that could not be attributed to the pulsar. All TS maps showed less than 1σ TS values in the region containing PSR J1826−1256 and HESS J1826−130. This non-detection of HESS J1826−130 is in agreement with Principe et al. (2020), who used off-pulse data from 1 GeV to 1 TeV and did not find evidence of residuals. We derived GeV flux upper limits for 4FGL J1826.1−1256 using the sed() method of Fermipy. In this method, the normalization is fit in each energy bin using a power-law model with the index = 2. The 95% confidence level upper limits were evaluated from the profile likelihood, which yielded a 30-52 GeV upper limit of 8.49 × 10 −13 erg cm −2 s −1 and a 52-91 GeV upper limit of 1.51 × 10 −12 erg cm −2 s −1 . We use these flux upper limits in our spectral energy distribution (SED) fits in §8- §8.3.
4. X-RAY OBSERVATIONS AND DATA REDUCTION PSR J1826−1256 and the surrounding Eel pulsar wind nebula have been observed by X-ray telescopes including NuSTAR and XMM-Newton.

NuSTAR
NuSTAR performed observations of PSR J1826−1256 on October 25, 2019 (ObsID 40501006002). The total exposure time was 76.5 ksec. NuSTAR's field of view (FOV) of 13 ×13 frames both the pulsar and PWN on the sub-detector DET0 for both focal plane modules (FPMA and FPMB). The NuSTAR data were processed and analyzed with NuSTARDAS Version 1.8.0 which is contained within HEASOFT v6.27.2, and we used the NuSTAR Calibration Database (CALDB) files v20191219.

XMM-Newton
XMM-Newton performed a 140 ksec observation of PSR J1826−1256 on October 11, 2014 (ObsID 0744420101). Data reduction and spectral extraction were performed in accordance with XMM-SAS v18.0 procedures. This observation was performed with the EPIC-MOS (European Photon Imaging Camera Metal Oxide Semiconductor) CCD detectors MOS1 and MOS2 in full-frame mode with the medium optical filter, which frames the entire diffuse PWN. The pn CCD detector was operated in the Small Window mode (FOV of 4'×4'), which is sufficient for pulsar and compact PWN analysis. Using calibrated event files, the data were reduced with the emchain and epchain tools. Using espfilt, good time intervals (GTIs) were found that excluded soft proton flares, which resulted in reduced exposure times of 77.4 ksec for MOS1, 80.8 ksec for MOS2, and 49.8 ksec for the pn. Canned FLAG values were used for MOS reduction, while FLAG was set to 0 for the pn. The corrected 0.5-10 keV images shown in Figure 2 were produced using the tasks eimageget and eimagecombine, with adaptive smoothing set to a radius between 1 and 4 pixels (1.1 per pixel). Figure 3 shows a 90 cm VLA field of the Eel PWN, containing numerous unidentified radio sources and large-scale emission across the region. The flux appears slightly increased near PSR J1826−1256 and the diffuse Eel PWN, however due to the large-scale diffuse emission we were unable to obtain an accurate flux for the PWN alone. We derived an upper limit for the 90 cm PWN emission from the measured flux density in the diffuse PWN region (marked as a red ellipse in Figure 3). We use this upper limit of 1.89 Jy in our SED fits in §8- §8.3. We estimate (using our evolutionary PWN modeling, see §8.3 with input parameters in Table 3) that the diffuse Eel PWN's contribution to the 330 MHz flux density is 15.9 mJy, and at 1.4 GHz is 9.4 mJy. Whether the diffuse Eel PWN emission can be directly detected or isolated at 330 MHz or 1.4 GHz remains a question we will investigate in future work.  from the pulsar, indicating evidence of the synchrotron burn-off effect. In order to further study the energy dependence of the source extension, we took advantage of NuSTAR's broadband capabilities and analyzed images in the 3.0-10.0 (soft) keV and 10.0-20.0 keV (hard) bands. The background was found to dominate above 20.0 keV. We used the HEASoft extractor (v5.36) tool to filter the cleaned event files into the soft and hard bands. Exposure maps were generated using the nuexpomap tool from nupipeline (v0.4.7), accounting for vignetting at 6.5 keV (soft images) and 15 keV (hard images). The four event files (FPMA and FPMB in each energy band) and their exposure maps were imported into XIMAGE (v4.5.1). Both the images and the associated exposure maps were smoothed using a σ = 3.69 (= 1.5 pixels) Gaussian kernel. We then applied exposure correction by dividing each source image by its exposure map. Exposure corrected images from the FPMA and FPMB in each energy band were summed to generate mosaic images, after correcting for the relative offset between the FPMA and FPMB images by using the brightest pixel as a reference point. The processed images are shown in Figure 4.

RADIO OBSERVATIONS
To study source extension at different energies, radial profiles (0 ≤ r ≤ 80 ) of the soft and hard images were generated around the brightest X-ray spot corresponding to the pulsar location. Since the background was found to be spatially uniform near the source, the difference between the average flux and the average value of the NuSTAR point spread function (PSF) at 60 ≤ r ≤ 80 was subtracted from the radial profile to measure source extension. A comparison with the onaxis NuSTAR PSF reveals extended emission from the compact PWN in both energy bands ( Figure 5), with a significance > 5σ from 18 -35 in the soft band and a significance > 3σ from 10 -28 in the hard band. We also found that the soft image is more extended by > 3σ significance relative to the hard image from 25 to 35 (reaching a maximum difference of 5.4σ at 33 ), consistent with the synchrotron burn-off effect suggested by Duvidovich et al. (2019). The soft and hard band profiles follow each other up to r ∼ 25 likely because the X-ray emission in this region is dominated by the pulsar component. The dip in the hard band below the NuS-TAR PSF at 45 -60 has less than a 2σ significance and is likely due to random fluctuations.
We further investigated 1-D spatial profiles in the NuSTAR images along two 20 wide axes, corresponding to "jet" and "torus" features reported by Karpova et al. (2019) using Chandra data. We found an excess at 25 NE of the pulsar position in the soft band image (perhaps corresponding to the "jet" feature). However, this feature was detected only with a significance of 2σ likely because NuSTAR's PSF is too broad to resolve the PWN substructures detected by Chandra.

X-RAY SPECTRAL ANALYSIS
Here we extract and analyze spectra from the NuS-TAR data alongside archived XMM-Newton data. We fit the NuSTAR and XMM-Newton spectra jointly in order to constrain the pulsar emission from PSR J1826−1256 and the compact Eel PWN. We then investigate the faint, diffuse Eel PWN to obtain X-ray data for SED fitting (see §8).

The Compact Eel PWN
To perform hard X-ray spectral analysis on the compact Eel PWN, source extraction regions were determined with physical justification. NuSTAR's angular resolution (58 HPD, 18 full width half maximum) is limited in comparison to XMM-Newton's, therefore a larger extraction radius is required to account for the wider PSF. The bright compact PWN radius is 15 as measured by Chandra (Karpova et al. 2019), therefore we determined that a NuSTAR extraction radius of 40 maximized PWN counts while minimizing counts from the diffuse X-ray nebula and background. We did not find a significant difference between using 30 and 40 NuSTAR extraction radii for our source region in the compact PWN analysis.
Annuli were avoided for background regions due to the diffuse PWN flux and high background levels in the NuSTAR data. Instead, we used source-free background regions on the same detector chip, which can be seen in Figure 4. NuSTAR spectral extraction was performed with nuproducts and response files were generated for an extended source. The spectra were binned to 20 counts per bin. To obtain XMM spectra, we used a 15 radius source spectral extraction region (the size of the bright, compact PWN) centered on PSR J1826−1256's location for MOS1, MOS2, and pn. We followed the standard XMM-SAS v18.0 procedures, and auxiliary response files (ARFs) were generated for an extended source.
Before fitting, the spectra were cut off above and below energies where the background began to dominate. For NuSTAR this gave a range from 4 to 20 keV, while the MOS1/2 spectra were fit between 1 to 8 keV and pn between 1 to 10 keV. For NuSTAR, this yielded 790 total and 562 net counts for FPMA, and 837 total and 475 net counts for FPMB. For XMM, this yielded 403 total and 367 net counts for MOS1, 422 total and 380 net counts for MOS2, and 816 total and 655 net counts for pn. Using XSPEC, we initially fit the NuSTAR + XMM-Newton spectra to an absorbed power-law model using a constant to normalize across the two telescopes: const*tbabs*(powerlaw+powerlaw). Because neither NuSTAR nor XMM can resolve the pulsar from the compact PWN, we added a second power-law to constrain the pulsar emission. Attempts to isolate the off-pulse emission in the NuSTAR observation using a timing solution were unsuccessful due to an unreliable 2σ pulsation detection. By freezing the pulsar's normalization and photon index, we were able to fit for the compact PWN's spectral parameters alone in both NuSTAR and XMM. For the pulsar power-law component, we froze Γ = 0.9 and a power-law normalization of 6.5 × 10 −6 , both obtained from the spatially-resolved Chandra spectral analysis in Karpova et al. (2019). This fit ( Figure 6) yielded a photon index Γ of 1.93 ± 0.19, a column density N H = (2.36 ± 0.34) × 10 22 cm −2 , and a reduced χ 2 = 1.07 for 151 degrees of freedom. The column density N H is consistent with Duvidovich et al. (2019) and Karpova et al. (2019) within error bars, while the photon index is slightly softer than the compact PWN photon index found by Karpova et al. (2019) using Chandra alone. All fitting parameters can be found in Table 1. We did not find any obvious evidence of a spectral break in the NuSTAR data between 10-20 keV. The hard Xray band spectra obtained by NuSTAR are shown in Figure 6, with no significant residuals when fit to a powerlaw model (note the NuSTAR-only fit parameters can be found in Table 1).

The Diffuse Eel PWN
To investigate the properties of the diffuse, cometary Eel PWN in X-rays, we analyzed XMM-Newton and NuSTAR spectra of the diffuse PWN region (marked by a red ellipse in Figures 1-3). XMM-Newton analyses of this region were first reported by Duvidovich et al. (2019) and Karpova et al. (2019). We followed the analysis procedures as used for the compact PWN in §4.2 and §7.1. The diffuse PWN region is defined as an ellipse of radii 3.5 × 1.6 , angled to avoid pulsed emission from PSR J1826−1256. The XMM-SAS v18.0 procedures were followed as before, and ARFs were generated for an extended source. The background was found to dominate outside of 1-6 keV, and fitting in this energy range yielded 8052 total and 4485 net counts for MOS1, and 8536 total and 4584 net counts for MOS2. The pn camera was operated in small-window mode for this observation and did not record counts from the entire diffuse PWN region. Fitting an absorbed power-law tbabs*powerlaw yielded a lower N H of 1.33 ± 0.15 × 10 22 cm −2 , a photon index Γ = 2.03 ± 0.15, and a reduced χ 2 = 0.84 for 155 degrees of freedom. The XMM-Newton MOS spectra are presented in Figure 7 after rebinning the data with a minimum 3σ significance in each bin.
Since NuSTAR did not significantly detect the diffuse PWN, we determined a 68% confidence interval upper limit for the diffuse PWN emission above 10 keV. We adopted only the FPMB data, since the FPMA data suffer from higher stray-light background contamination. We extracted NuSTAR FPMB spectra from the diffuse PWN region used for all XMM analysis (the red ellipse in Figures 1-3) and fit an unabsorbed power-law model with Γ = 2.0, the best-fit photon index of the diffuse PWN obtained by our fit to the XMM-Newton data (Table 1). We raised the flux normalization to match the background flux levels measured by NuSTAR, yielding 10-20 keV and 20-30 keV flux upper limits of 1.0 ×10 −12 erg cm −2 s −1 and 0.6 ×10 −12 erg cm −2 s −1 , respectively. These NuSTAR flux upper limits above 10 keV are presented in the subsequent SED plots along with XMM-Newton spectra of the diffuse PWN.

BROADBAND PWN SPECTRAL MODELING
Here we present SED modeling of the TeV source data using three different models -(1) a pure leptonic model, (2) a pure hadronic model, and (3) an evolutionary leptonic PWN model. In all sections, we adopted the XMM-Newton MOS1+2 spectra and NuSTAR hard X-ray flux upper limits of the diffuse PWN (the elliptical region shown in Figure 2) as well as the Fermi-LAT upper limits at the TeV source position ( §3) and VLA 90 cm radio upper limit of the diffuse PWN region ( §5). As shown in Table 1, the diffuse PWN dominates in total X-ray flux, while the compact PWN's contribution is less than 10% the diffuse PWN's flux (or more precisely ∼4% if the pulsar component is excluded; Karpova et al. 2019). In addition, the diffuse PWN is spatially coincident with the H.E.S.S. source, with the peak 5σ flux of HESS J1826−130 aligned with a region of the diffuse PWN ∼6 offset from PSR J1826−1256 (see H.E.S.S. contours in Figure 2). For these reasons, we use diffuse PWN X-ray spectra in the following SED analysis.
HAWC J1826−128 and HESS J1826−130 are offset in flux by a factor of ∼1.8-2, which could be explained by some combination of statistical and systematic errors, including both telescopes' modeling of the Galactic diffuse TeV emission (Abdalla et al. 2021). For this reason, and to avoid arbitrary scaling of one source's flux to the other, we fit our leptonic, hadronic, and evolutionary PWN models to the HAWC J1826−128 and HESS J1826−130 TeV datasets separately in order to explore and compare the resulting parameter spaces.
In all fits, we use HAWC J1826−128 and HESS J1826−130 spectral data from Albert et al. (2021) and H. E. S. S. Collaboration et al. (2020), respectively. In §8.1, we apply the Naima package for modeling SEDs (Zabalza 2015) to a pure leptonic scenario, considering the VLA 90 cm, XMM-Newton, NuSTAR, Fermi-LAT   (155)  GeV, and HAWC/H.E.S.S TeV data. In §8.2, we apply a pure hadronic model (also available in Naima) to the TeV data only, not using the X-ray data on account of their association with the (leptonic) diffuse PWN. In §8.3, after ruling out the hadronic case and favoring a leptonic interpretation, we adopted a more sophisticated leptonic SED model that accounts for time evolution of the multi-wavelength PWN emission ), jointly fitting VLA 90 cm, XMM, NuSTAR, Fermi-LAT GeV, and HAWC/H.E.S.S TeV data.

Leptonic Model Fitting in Naima
A leptonic SED model generally consists of two primary components from radio to TeV energies. The lower-energy emission (up to ∼MeV) is created by synchrotron radiation of primary electrons accelerated in the PWN. The GeV and TeV γ-ray emission is produced by Inverse Compton scattering (ICS) of the CMB and ambient photons by the same electrons. In addition, the model contains a synchrotron self-Compton (SSC) component which takes into account synchrotron photons up-scattered by energetic electrons inside the PWN. However, the SSC component is usually less significant than the ICS component in all but the most compact nebulae.
We constructed our leptonic SED model using the InverseCompton and Synchrotron components in Naima. We assumed an exponential cut-off powerlaw (ECPL) model for the electron energy distribution (f (E) ∝ AE −α e −E/Ecut electrons cm −3 erg −1 ). Naima's SED data fitting is performed using a Markov Chain Monte Carlo (MCMC) algorithm which samples posterior distributions of a set of model parameters, iteratively optimizing their respective χ 2 values. The free parameters that we fit to VLA, XMM, NuSTAR, Fermi-LAT, and either HAWC or H.E.S.S. data are a powerlaw index α, cut-off energy, flux normalization of the PWN electron energy distribution, PWN B-field, and FIR photon energy density. Voisin et al. (2016) investigated the ISM in this region and analyzed CO lines and maser emission consistent with molecular clouds. Assuming the large optical depth limit for the 12 CO line, we derived a low excitation temperature range of T exc < 20 K from the antenna temperatures reported for the nearby CO emission region labeled as "R1c" Voisin et al. (2016). We assumed T = 15 K as the temperature of the local FIR photon field and left its energy density as a free parameter.
When we fit the HAWC data alongside the X-ray data, we found that the best-fit model yields a power-law index of α = 2.2 and a cut-off energy of E cut = 590 TeV. The best-fit B-field is low at 1.9 µG, and the regional FIR energy density is elevated at 4.4 eV cm −3 in order to match the TeV fluxes. The total electron energy injected into the PWN is W e = 1.8 × 10 46 erg. The best fit to the HESS J1826−130 TeV data yields a slightly softer power-law index α = 2.5 as well as the same cutoff energy E cut = 590 TeV. The B-field is slightly higher at 3.0 µG, and the resulting FIR energy density is lower at 2.9 eV cm −3 . The total electron energy required in 10 6 10 3 10 0 10 3 10 6 10 9 10 12 10 15 Energy  Figure 8. A Naima leptonic model fit to the broadband SED with radio, X-ray, GeV upper limits, and HAWC J1826−128 spectral data. The electron energy distribution is assumed to be an exponential-cutoff power-law model. The fit parameters can be found in Table 2. 10 6 10 3 10 0 10 3 10 6 10 9 10 12 10 15  Figure 9. A Naima leptonic model fit to the broadband SED with radio, X-ray, GeV upper limits, and HESS J1826−130 spectral data. Though the HAWC and H.E.S.S. observations likely detected the same source, their γ-ray fluxes are slightly offset from one another, potentially as a result of systematics or differing analyses between the two telescopes (Abdalla et al. 2021). The fit parameters can be found in Table 2. . this model is W e = 1.6 × 10 46 erg. Our SED data and best-fit models are shown in Figures 8-9, and all fit parameters are listed in Table 2. The similarity between SED fitting parameters combined with their spatial association strongly suggests that HAWC J1826−128 and HESS J1826−130 are two separate detections of the same TeV γ-ray accelerator, as previously noted by Albert et al. (2021). Given the successful fitting with this 10 6 10 3 10 0 10 3 10 6 10 9 10 12 10 15

Hadronic Model Fitting in Naima
In the hadronic case, relativistic proton collisions with the ISM produce both neutral and charged pions. Energetic neutral π 0 mesons can decay into photons above ∼100 TeV, making hadronic processes a plausible mechanism for some PeVatron candidates (such as SNR G106.3+2.7; Tibet ASγ Collaboration et al. 2021). Additionally, the charged pions decay into secondary electrons and positrons, which contribute to ICS γ-ray emission and synchrotron radiation at X-ray energies and below. Since the X-ray emission around the pulsar is most likely associated with leptonic processes in both the compact and diffuse PWN regions, we do not include the X-ray data obtained by XMM and NuSTAR in our hadronic model fit. Instead, the observed X-ray emission can serve as an upper limit for the diffuse synchrotron X-ray emission expected from hadronic secondary electrons or hadronic thermal X-ray emission.
To investigate the hadronic scenario, we fit the hadronic PionDecay model, available in the Naima package, to both the HAWC and H.E.S.S. TeV data. We assumed an ECPL model for the proton energy distribution (f (E) ∝ AE −α e −E/Ecut protons cm −3 erg −1 ). As a target for the hadronic collisions, we considered the CO emission region at 3.5 kpc distance Voisin et al. (2016) labeled as "R1c", a kinematic velocity component of region "R1" which is consistent with the distance to PSR J1826−1256 found by Karpova et al. (2019). We derived a mean hydrogen density n H = 600 cm −3 from the n H2 = 220 cm −3 virial mass of the CO line emission component and the relation n H = 2.8 n H2 (Voisin et al. 2016). The best-fit model using the HAWC J1826−128 data yields a power-law index of α = 1.5 and a cutoff energy of E cut = 100 TeV. The total proton energy required in this model is W P = 1.0 × 10 47 erg. The HESS J1826−130 TeV data can be fit with the same power-law index of α = 1.5 and a lower cut-off energy E cut = 40 TeV due to the lack of higher energy bins. The total proton energy required in this model is W P = 3.0 × 10 46 erg. Note that Naima does not model the secondary electron energy distribution self-consistently from the proton-proton collisions in the PionDecay model component. Hence, following the method of Mori et al. (2020), we added separate Synchrotron and InverseCompton components to account for the contribution of secondary electrons in the X-ray and TeV bands. We assumed the same energy spectrum for protons and secondary electrons (electron α = 1.5 and cut-off energy E cut = 100 TeV for HAWC J1826−128). We adjusted the flux normalizations of the synchrotron and ICS components by imposing that the total energy of secondary electrons (W e ) is 1/3 of the total energy of injected protons, which accounts for the branching ratios of charged pions produced from p-p collisions (Coerver et al. 2019;Mori et al. 2020).
As can be seen in Figure 10, we found that the predicted synchrotron emission component from the secondary electron population, obtained by fixing the W e /W p ratio to 1/3, vastly overpredicts the PWN Xray flux even if a low ISM B-field of 3 µG (Crutcher 2012) is assumed. In this model, any B-field > 1 µG overpredicts the X-ray flux, while B < 1 µG is unable to match the observed X-ray spectral shape. The lack of  Figure 11. An evolutionary PWN model fit of the broadband SED with HAWC J1826−128 TeV data (dark red) and Fermi-LAT upper limits (purple), along with XMM-Newton X-ray data, NuSTAR upper limits (green), and a VLA 90cm upper limit (red). The fit parameters can be found in Table  3. , XMM-Newton X-ray data, NuS-TAR upper limits (green), and a VLA 90cm upper limit (red). The fit parameters can be found in Table 3.
this predicted diffuse X-ray emission in the TeV source region (which would appear substantially higher in flux than the compact or diffuse PWN at typical ISM Bfields), combined with the presence of an energetic pulsar and diffuse PWN coincident with HAWC J1826−128 and HESS J1826−130, rules out the pure hadronic case as the sole origin of the TeV emission.

Evolutionary PWN Model
In order to estimate the physical parameters of the progenitor supernova, associated neutron star, and its pulsar wind in the system, we fit the observed properties of the diffuse Eel PWN with the predictions of a physical model for the dynamical and radiative evolution of a spherical PWN inside a SNR (see Gelfand et al. (2009) for a detailed description of this particular model, and Gelfand (2017) for a review of similar dynamical, radiative PWN models).
In order to obtain posterior minimum χ 2 parameters for the Eel PWN, we made certain assumptions within this model to account for uncertain or unconstrained features of the pulsar, SNR and ISM. Firstly, we assumed that the unshocked SN ejecta density profile features a uniform density core which is surrounded by an envelope with ρ ∝ r −9 . This is a common assumption, but may not be fully representative of the varied ejecta density profiles of various progenitor SNRs (Gelfand 2017). Next, we assumed that the supernova ejecta with initial kinetic energy E sn and mass M ej is expanding into a medium with uniform density n ism . Using the observed current spin-down luminosityĖ and characteristic age t ch of PSR J1826−1256 as inputs for the model (and ensuring that each trial combination of input parameters reproduced these observables), for a given braking index p and spin-down timescale τ sd we determined the system's true age t age and initial spin-down luminositẏ E 0 as follows: (2) Next, we assumed that the entire spin-down luminos-ityĖ of PSR J1826−1256 is injected into the PWN as either magnetic fieldsĖ B or the kinetic energyĖ p of relativistic leptons: holding η B constant with time. Next, we assumed that the spectrum of particles injected into the PWN can be described by a broken power-law: where we further assumed that the free parameters E min , E break , E max , p 1 and p 2 are constant with time. The power-law normalizationṄ break across all times t is calculated by the relation: Lastly, we assumed that synchrotron and inverse Compton scattering are the only two mechanisms for electron radiative losses within the PWN, and that the B-field across the diffuse Eel PWN has a uniform strength (see Gelfand et al. 2009).
The total electron energy, magnetic field strength, and radius of the SNR reverse shock R RS are output by the model. All input and output parameters for this model are listed in Table 3. No priors were used in our model fitting -the possible value distribution for all parameters is assumed to be unbounded and completely flat (with the exception of the wind magnetization, which as a fraction necessarily must be < 1). A Metropolis Markov Chain Monte Carlo (MCMC) algorithm (see Gelman et al. (2013) and Section 3.2 of Gelfand et al. (2015) for more details) was used to determine the combination of parameters which best reproduced the observed X-ray and γ-ray fluxes of the diffuse Eel PWN, the observed radius of the coincident radio SNR candidate G18.45−0.05 (7.0 ± 0.6 ; Anderson et al. 2017), and a circularly-approximated PWN radius. As a best estimate for the full size of the γ-ray emitting PWN, we used our measured 3σ extension of the H.E.S.S source at E > 10 TeV (Figures 1 and 2) of 6 or ∼6 pc at D = 3.5 kpc in all model fits, with ±1 pc uncertainty.
The combinations of input parameters which resulted in the lowest χ 2 (HAWC J1826−128: χ 2 = 22.85 for 22 degrees of freedom and HESS J1826−130: χ 2 = 33.6 for 27 degrees of freedom) are given in Table 3. As shown in Figures 11 and 12, for this set of model parameters, inverse Compton scattering CMB photons alone is sufficient to explain the observed γ-ray emission. The addition of a low energy density photon field is compatible with this set of model parameters but did not improve the quality of the fit.
Several SNR/PWN systems have been found with comparable non-canonical parameters (e.g. our predicted progenitor SN explosion energy and ejecta mass, see Table 3). For example, Guest et al. (2019) derived a SN explosion energy of 3 × 10 49 erg in their analysis of the G21.5−0.9 progenitor SNR, and Gotthelf et al. (2021) found a low ejecta mass of 0.51 M and a SN explosion energy of 1.26 × 10 50 erg in their application of Gelfand et al. (2009)'s evolutionary model to the SNR Kes 75. Additionally, it has been suggested that the Crab progenitor SN was of an unusually low energy (∼10 50 erg; Yang & Chevalier 2015). A pulsar braking index > 3 as suggested by our fit to HAWC J1826−128 is unusual, but such braking indices have been measured for anomalous pulsars such as PSR J1640−4631 and could be explained by invoking higher order multipoles (Archibald et al. 2016). The HESS J1826−130 fit braking index is of a typical magnitude, ≤ 3.
As detailed by Gelfand et al. (2015) in their analysis of G54.1+0.3, there can be significant degeneracies between the different parameters of this model, for example between the mass and initial kinetic energy of the supernova ejecta with the density of the surrounding ISM, and between the braking index and spin-down timescale of the pulsar. Further parameter exploration is needed to derive proper statistical bounds on the values of these quantities. Further investigation of the degeneracy between braking index and spin-down timescale may allow for lower braking indices when fitting the HAWC spectrum. The pulsar spin-down timescale in our modeling is higher than most other systems, with best-fit values of 8.4 and 9.8 kyr for the HAWC and H.E.S.S. sources, respectively. However, other time-evolution fits to PWNe have yielded high spin-down timescales -for example, Torres et al. (2014) modeled a spin-down age of 6.6 kyr in a fit to HESS J1356-645, an offset relic TeV PWN, with a pulsar of similarĖ and comparable PWN morphology to the Eel.
The ISM density, HESS J1826−130 fit pulsar braking index, and wind magnetization parameters all agree with previous applications of time-evolution PWNe models (Torres et al. 2014;Gelfand et al. 2015). The low and maximum energy particle indices are also in agreement with the single power-law index found by Naima. Furthermore, the spherically-approximated PWN radius (6 pc) and SNR (∼8-10 pc) predicted by this set of model parameters agree well with the observed PWN volume and SNR candidate radius.  Roberts et al. (2007). HAWC J1826−128 and HESS J1826−130 likely represent independent observations of the same TeV source (this paper; Albert et al. 2021) and its TeV emission (as resolved by H.E.S.S. above 10 TeV) spatially coincides with the diffuse X-ray PWN. Both the Naima 0.001 0.005 Min. injected particle energy 0.3 GeV 0.5 GeV Max. injected particle energy 2.3 PeV 1.6 PeV Break injected particle energy 28 TeV 3.4 TeV Low energy particle index 1.7 1.9 High energy particle index 2.6 2.8 ForĖ = 3.6 × 10 36 ergs/s and CMB-only photon field leptonic SED model (assuming a homogeneous population of relativistic electrons; §8.1) and the evolutionary PWN model ( §8.3) yield fit parameters that are consistent with the PWN hypothesis. However, in comparison with other middle-aged PWNe, there are several distinct features of the Eel PWN which may account for its energetic TeV emission as a possible PeVatron candidate. Based on the pulsar and PWN parameters determined by our SED analysis, we investigate the properties of the Eel PWN in this section.

PWN Energetics
While our independent leptonic SED models (Naima and Gelfand et al. 2009) yielded similar particle indices and low B-fields, they differ in the best-fit total electron energy and photon field energy density values. This is primarily because our evolutionary PWN model takes into account electron injection and cooling throughout the modeled pulsar true age, while the generic Naima model assumes no time evolution of the electron energy distribution. The total electron energy (W e = (1.6-1.8) × 10 46 ergs) obtained by the Naima model is an order of magnitude lower than that of the evolutionary PWN model (W e = (3.5-5.0) × 10 47 ergs). As a result, the ICS component in the Naima model requires an elevated FIR photon field energy density of 2.9-4.4 eV cm −3 (in addition to the CMB component) to fit the observed TeV flux.
However, the Naima results are likely incompatible with the observed TeV γ-ray extension of HESS J1826−130. Note that the electron cooling time due to pitch angle averaged synchrotron and inverse Compton energy losses in the Thomson regime can be expressed by: (Rybicki & Lightman 1979) where U ph and U B are the photon field and magnetic field energy densities, respectively, γ 7 = E/(10 7 m e c 2 ) is the particle Lorentz factor, and B 5 = B × 10 5 G is the magnetic field. Therefore, given the FIR energy density and PWN B-field (B = 1.9-3.0 µG) required by the Naima model, > ∼ 40 TeV electrons (which would emit 10 TeV ICS γ-rays in these conditions) will cool in < 2 kyr, preventing them from diffusing across the entire observed extension of HESS J1826−130 (see §9.3, (8) and (9)). On the other hand, in the case of the evolutionary PWN model, the cooling time of ∼60 TeV electrons (which would emit γ-rays above 10 TeV in a CMB-only photon field) is far greater than our modeled pulsar true age (4.6-6.7 kyr). This suggests that the injected electrons in our evolutionary PWN model should be able to emit TeV γ-rays throughout the true pulsar age. Dividing the total electron energy of (3.5-5.0) × 10 47 ergs by the spin-down power (Ė = 3.6 × 10 36 erg s −1 ) yields an injection timescale of 3.1-4.4 kyr which is comparable to our modeled true age. In the case of a cooling timescale > true age and injection timescale ∼ true age as suggested by the evolutionary PWN model parameters, 4.6-6.7 kyr also provides sufficient time for these electrons to diffuse across the observed TeV γ-ray extension of HESS J1826−130 ( §9.3). Therefore, in addition to the fact that the evolutionary PWN model is more physically motivated, we consider it plausible that the FIR photon field is not significant in the Eel PWN region and that the electron energy is on the order ∼10 47 ergs -these assumptions are adopted in the subsequent B-field analysis ( §9.3).
To find the total rotational energy released from the pulsar (E tot ) we integrated the spin-down power (Ė) over the true age suggested by our evolutionary model (t = 4.6-6.7 kyr). Next, dividing the total electron energy W e by E tot indicated the fraction of rotational energy that was spent injecting electrons during the PWN evolution. The evolutionary model fits to the HAWC J1826−128 and HESS J1826−130 data predicted an electron energy of (3.5-5.0) × 10 47 ergs, andĖ integrated over true age is (5.2-7.6) × 10 47 ergs, implying that approximately two-thirds of the total rotational energy has been injected into the PWN.
Given the pulsar's spin-down power (Ė = 3.6 × 10 36 erg s −1 ), we found that the 0.5-8 keV X-ray and 1-10 TeV TeV efficiencies of the Eel PWN are η keV = 5.5 × 10 −4 and η T eV = (2-4) × 10 −3 respectively, assuming a pulsar distance of 3.5 kpc. The TeV efficiency range reflects the different values for TeV luminosity obtained by the HAWC and H.E.S.S source data. These X-ray and TeV efficiencies are consistent with those of middle-aged PWNe (Kargaltsev et al. 2015). Additionally, Karpova et al. (2019) found that the TeV and X-ray luminosity and L T eV /L keV ratio versus system age were consistent with the luminosity and ratio-to-age distributions of other PWNe detected in the TeV band (H. E. S. S. Collaboration et al. 2018;Kargaltsev et al. 2013).

Comparison with other middle-aged PWNe
While the PWN hypothesis is well-supported, the exact properties of the diffuse Eel PWN remain uncertain. These uncertainties include the mechanism responsible for the highly asymmetrical X-ray morphology and the magnitude of the ambient B-field inside this extended X-ray and γ-ray tail. One possible explanation could be that the Eel PWN is produced by a supersonic pulsar, which has been suggested due to the PWN's diffuse Xray morphology and its similarity in appearance to other supersonic PWNe (Roberts et al. 2007;Kargaltsev et al. 2017). Another possible explanation could be an ongoing collision with the SNR reverse shock, which would prevent the outflow of radiating particles near the shock front and allow all remaining particles to diffuse away freely in the opposite direction, creating an asymmetric PWN which is often found to be offset from the pulsar position (Blondin et al. 2001;Karpova et al. 2019). Table 4 shows a comparison of the Eel PWN to three other evolved PWNe which are also characterised by a compact X-ray nebula around the pulsar and diffuse Xray emission. MSH 15−56 and G327.1−1.1 are both thought to have interacted with the SNR reverse shock, resulting in soft X-ray thermal line emission and radio brightness, while the origin of CTB 87's asymmetric morphology remains uncertain (Guest et al. 2020). Unlike the other three evolved PWNe, the Eel PWN lacks any significant radio emission (a signature often seen in PWNe interacting with SNR reverse shocks; Blondin et al. 2001), with our estimated 1.4 GHz flux density just 0.04% of CTB 87 and even less for MSH 15−56 and G327.1−1.1. Considering our modeled true age of 4.6-6.7 kyr (Table 3) compared to 17-20 kyr in the other PWNe, the Eel may not have had sufficient time to undergo substantial reverse shock interaction.
In addition to the lack of radio emission around the Eel PWN, our one-zone evolutionary model in its spherical approximation of the PWN and SNR also suggests that the PWN-reverse shock interaction has either not yet begun, or is just starting (see Table 3). While the pulsar's proper motion may cause one side of the PWN to interact with the reverse shock first and eventually result in an asymmetrical morphology, we did not find any significant evidence of thermal emission in our compact PWN spectra (see Figure 6). Instead, the compact Xray nebula is characterised by non-thermal X-ray emission with evidence of the synchrotron burn-off effect. The lack of features such as thermal emission lines in the compact PWN, an X-ray bright SNR, and a bright relic radio PWN coincident with the diffuse Eel PWN (all features seen in the reverse shock-interacting PWNe MSH 15−56 and G327.1−1.1; Temim et al. 2013Temim et al. , 2015 alongside the uncertain results of our modeled reverse shock prediction suggest that the Eel PWN may not have yet commenced a substantial interaction with the reverse shock of its host SNR. Therefore, the supersonic PWN hypothesis remains a plausible explanation for the asymmetric, diffuse X-ray emission seen extending from PSR J1826−1256.

Magnetic Field of the Diffuse PWN
The magnetic field strength B inside the diffuse Eel PWN (the location of a majority of the TeV-emitting electrons) has not been previously well constrained. Duvidovich et al. (2019) suggested that B > 2 µG based on   2019), 0.5-10 keV flux in erg cm −2 s −1 ; e this paper, see §7.2; f this paper, estimate from evolutionary PWN model; g Ma et al. (2016); h modeled true age vs pulsar characteristic age; i relic PWN age the peak energy of the diffuse X-ray spectrum and the assumption that the synchrotron cooling time must be less than the pulsar's characteristic age. Karpova et al. (2019) suggested B ∼ 5 µG on the basis of particle diffusion lengths across the angular size of HESS J1826−130. Similarly, we find that low B-fields (∼0.7-3 µG) are required to jointly fit the X-ray and TeV data in both the Naima leptonic SED and evolutionary PWN models (see Tables 2 and 3). Low B-fields have previously been found in other TeV-detected PWNe: fitting a different evolutionary PWN model to multi-wavelength SED data of G292.2-0.5 and HESS J1356-645 yielded B = 4 µG and 3.1 µG, respectively (Torres et al. 2014). Low wind magnetization fractions (we obtained typical values in our evolutionary model fits of 0.1-0.5%, see Table 3) indicate that only a small fraction of the pulsar's spindown luminosity is injected into the PWN in the form of magnetic fields. Additionally, adiabatic expansion of the PWN decreases its magnetic field strength, as shown by Gelfand et al. (2009), Tanaka & Takahara (2010, and other similar works. A lower PWN B-field increases the synchrotron lifetime of the diffusing particles, which is consistent with the large angular extensions of the Eel PWN in the X-ray and TeV bands. Our Naima parameters are derived from generic nonthermal particle modeling that does not account for time-evolution or cooling. For this reason, below we investigate the cooling and diffusion of electrons using the PWN parameters obtained by the evolutionary PWN model. For all transverse distances and diffusion lengths, we assumed that the pulsar distance is 3.5 kpc. While synchrotron X-ray radiation of leptons in the diffuse PWN falls below XMM-Newton sensitivities at distances greater than 4-6 from the pulsar, this X-ray faint region is spatially associated with the brightest γrays from HESS J1826−130 (see the offset location of the HESS J1826−130 5σ peak flux, Figure 2). This off-set between the diffuse PWN detected by XMM-Newton and γ-ray centroid indicates that the true extension of the diffuse PWN is likely larger than 4-6 , with some PWN particles reaching the HESS J1826−130 peak 6-8 from the pulsar. This suggests long diffusion lengths for ICS leptons in the PWN. Additionally, the ASCA source catalog in 2001 reported a diffuse X-ray source extended over 8 from PSR J1826−1256 (Roberts et al. 2001). As a best estimate for the diffusion lengths required for particles to fill the observed TeV source at a given energy, we used our measured 3σ extension of the H.E.S.S source at E > 10 TeV (Figures 1 and 2) of 6 or ∼6 pc at D = 3.5 kpc. Given that the H.E.S.S source spectrum peaks at E ∼ 10 TeV, we assumed that the ICS γ-ray energy at this radius is typically E ICS ∼ 10 TeV.
CMB photons alone are sufficient to account for the TeV emission in our evolutionary PWN model. In a CMB-only photon field, 63 TeV electrons would emit 10 TeV γ-rays via ICS. These ∼60 TeV electrons have a cooling time (7) of 15 kyr for the derived B-field of ∼1 µG. Next, we calculated how far these energetic electrons can diffuse within their cooling time with this formula: where t d is the diffusion timescale after escape and R d is a lower limit for diffusion length. Assuming Bohm-type diffusion in the region surrounding PSR J1826−1256, the diffusion coefficient D B can be written as (9) Using our evolutionary PWN model parameters (Table 3), 63 TeV electrons within B = 0.6-1.5 µG have a cooling time (∼15 kyr) greater than our modeled true pulsar age of 4.6-6.7 kyr. Using this true age to calcu-late a diffusion radius instead of the cooling time, the 63 TeV electrons are able to diffuse to 6 pc, which is consistent with the HESS J1826−130 extension above 10 TeV.
If the pulsar is fast moving, its proper motion could modify the size of the X-ray and TeV emission region. However, we find that the pulsar will travel between 2-7 pc within our estimated true pulsar age, assuming a transverse velocity between 500-1,000 km s −1 . In either case, the estimated diffusion length (plus an additional distance due to the pulsar's proper motion) is comparable to the measured 3σ extension of HESS J1826−130 above 10 TeV. Therefore, our estimated PWN B-field of ∼1 µG is well justified due to its consistency with both the broadband SED data and the TeV source size.
Our estimated diffuse PWN B-field is also compatible with estimations of electron diffusion in the X-ray band. The peak of the diffuse PWN's synchrotron Xray spectrum (1.8 keV, see Figure 7) indicates that 210 TeV electrons will radiate 1.8 keV photons in a ∼1 µG B-field, with electron energy inversely proportional to B-field. These electrons will cool over a timescale longer than our estimated true pulsar age in a CMB-only environment, therefore using the true age (4.6-6.7 kyr) as a diffusion timescale suggests that electrons radiating at ∼2 keV will be able to diffuse > 6 pc in a B-field of ∼1 µG, which is consistent with both the broadband SED modeling and the observed diffuse PWN size in X-rays.

The Eel as a PeVatron PWN?
In fitting our evolutionary PWN model to the HAWC J1826−128 and HESS J1826−130 data over a large parameter space, we found that a maximum electron energy E max > ∼ 1 PeV is strongly preferred with the observed PSR J1826−1256 spin-down power ofĖ = 3.6 × 10 36 erg s −1 as a fixed input. This is because a lower E max requires a large fraction of the spin-down power to be injected into lower energy electrons, thus forcing our time-evolution SED model to over-predict the observed γ-ray fluxes at E < 10 TeV and the GeV upper limits. The observed TeV γ-ray fluxes and GeV upper limits combined with the pulsar's high spin-down power make a robust case for E max > ∼ 1 PeV in our multiwavelength evolutionary SED modeling of the Eel PWN, and parameter exploration with E max < 1 PeV did not yield a successful fit. Even without time dependence, the generic Naima leptonic model also suggests high electron energies, requiring a cut-off electron energy of E max ∼ 600 TeV. While there are some uncertainties associated with our model assumptions (e.g. representing an asymmetrical PWN as spherical), both the leptonic Naima model and our evolutionary PWN model sug-gest the acceleration of leptons to energies on the order of hundreds of TeV up to a maximum of ∼2 PeV.

CONCLUSION
Recent γ-ray and X-ray observations have suggested that the Eel PWN is a counterpart to the bright HAWC source HAWC J1826−128, also detected by H.E.S.S primarily above E ∼ 1 TeV. The HAWC detection up to nearly E ∼ 100 TeV suggests that the TeV source may be a PeVatron candidate, i.e., a long-sought, extreme particle accelerator in our galaxy (Albert et al. 2021). NuSTAR has revealed evidence of strong synchrotron burn-off in the compact X-ray nebula up to ∼20 keV, while a previous XMM observation detected tail-like diffuse X-ray emission extending ∼6 from the pulsar. The diffuse X-ray PWN coincides well with the H.E.S.S TeV source detected above 10 TeV, while the contribution of the compact PWN and the pulsar is less than 10% in the X-ray band.
Based on our multi-wavelength SED and morphology study, we ruled out a pure hadronic origin of the TeV emission due to the inability of a hadronic SED model to fit both the X-ray and TeV data (along with the spatial coincidence between the diffuse Eel PWN and HESS J1826−130 above 10 TeV). Through our time-evolution leptonic SED modeling as well as theoretical estimates of the electron cooling and diffusion timescales, we found that (1) the radiative efficiencies in the X-ray and TeV bands are comparable to those of other middle-aged PWNe; (2) the true pulsar age is ∼5.7 kyr (4.6-6.7 kyr) which is significantly smaller than the spin-down age of 14.4 kyr; (3) the PWN B-field is low at B ∼ 1 µG; (4) the asymmetric X-ray nebula morphology is likely due to supersonic motion of the pulsar; (5) the lack of radio emission is notable compared to other evolved PWNe with similar X-ray and TeV properties; (6) the maximum electron energy may extend above ∼1 PeV. Our modeled maximum electron energy reaching ∼1 PeV and beyond is particularly exciting, as it further suggests the Eel PWN may be a leptonic PeVatron.
Why the Eel PWN stands out as an extreme leptonic particle accelerator will remain an elusive question until other TeV-detected PWNe are further investigated in similar multi-wavelength approaches. These investigations are ongoing with seven other middle-aged PWNe (including several PeVatron candidates) detected by HAWC and LHAASO (Mori et al. 2021). In a larger context, it is still unclear whether evolved PWNe such as the Eel are dominant contributors to the Galactic Pe-Vatron population as well as the positron excess (Malyshev et al. 2009). As demonstrated in this paper, an extensive X-ray survey of middle-aged PWNe with NuS-TAR, XMM and other X-ray telescopes, along with use of next-generation TeV telescopes such as the Cherenkov Telescope Array (CTA) and the Southern Wide-field Gamma-ray Observatory (SWGO), will be essential to answering these fundamental questions in high-energy astrophysics.