# MeV-GeV Gamma-Ray Emission from SNR G327.1-1.1 Discovered by the Fermi-LAT

Jordan Eagle, Daniel Castro, Tea Temim, Jean Ballet, Patrick Slane, Joseph

Gelfand, Matthew Kerr, Marco Ajello

► To cite this version:

Jordan Eagle, Daniel Castro, Tea Temim, Jean Ballet, Patrick Slane, et al. MeV-GeV Gamma-Ray Emission from SNR G327.1-1.1 Discovered by the Fermi-LAT. Astrophys.J., 2022, 940 (2), pp.143. 10.3847/1538-4357/ac9eb4. hal-03882945

HAL Id: hal-03882945

https://hal.science/hal-03882945

Submitted on 3 Dec 2022

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
OPEN ACCESS

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

CrossMark

# MeV-GeV Gamma-Ray Emission from SNR G327.1-1.1 Discovered by the Fermi-LAT

Jordan Eagle<fnref n="1" />,<fnref n="2" /> ID, Daniel Castro<fnref n="1" /> ID, Tea Temim<fnref n="3" /> ID, Jean Ballet<fnref n="4" /> ID, Patrick Slane<fnref n="1" /> ID, Joseph Gelfand<fnref n="5" /> ID, Matthew Kerr<fnref n="6" /> ID, and Marco Ajello<fnref n="2" /> ID

Received 2022 August 26; revised 2022 October 14; accepted 2022 October 28; published 2022 November 30

## Abstract

We report the discovery of MeV-GeV γ-ray emission by the Fermi-LAT positionally coincident with the TeV pulsar wind nebula (PWN) HESS J1554–550 within the host supernova remnant (SNR) G327.1–1.1. The γ-ray emission is point-like and faint but significant (>4σ) in the 300 MeV–2 TeV energy range. We report here the Fermi-LAT analysis of the observed γ-ray emission followed by a detailed multiwavelength investigation to understand the nature of the emission. The central pulsar powering the PWN within G327.1–1.1 has not been detected in any wave band; however, it is likely embedded within the X-ray nebula, which is displaced from the center of the radio nebula. The γ-ray emission is faint and therefore a pulsation search to determine if the pulsar may be contributing is not feasible. Prior detailed multiwavelength reports revealed an SNR system that is old, τ ~ 18,000 yr, where the interaction of the reverse shock with the PWN is underway or has recently occurred. We find that the γ-ray emission agrees remarkably well with a detailed broadband model constructed in a prior report based on independent hydrodynamical and semianalytic simulations of an evolved PWN. We further investigate the physical implications of the model for the PWN evolutionary stage incorporating the new Fermi-LAT data and attempt to model the distinct particle components based on a spatial separation analysis of the displaced PWN counterparts.

*Unified Astronomy Thesaurus concepts: Gamma-ray astronomy (628); Pulsar wind nebulae (2215); Supernova remnants (1667); X-ray astronomy (1810); Spectral energy distribution (2129); Radiative processes (2055)*

## 1. Introduction

Pulsar wind nebulae (PWNe) are descendants of core-collapse supernovae (CC SNe), and are powered by an energetic, rapidly rotating neutron star. As the neutron star loses spin, rotational energy is transferred into an ultra-relativistic particle wind, made up of mostly electrons and positrons. A shock wave, called the termination shock, is generated within the nebula where the ram pressure of the pulsar wind is balanced by the pressure of the highly relativistic plasma (Gaensler & Slane 2006; Slane 2017). The termination shock is the location where particles injected by the central pulsar are accelerated as they enter the non-thermal pool of relativistic particles in the nebula. Synchrotron emission from the relativistic electrons is observed from most PWNe, from radio wavelengths to hard X-rays, while the same relativistic electrons are thought to scatter off of local photon fields, resulting in inverse Compton (IC) emission at γ-ray energies. Most PWNe have been discovered at radio or X-ray bands with an increasing number of discoveries in TeV. Indeed, the majority of the TeV Galactic source population is found to originate from PWNe as observed by Imaging Air Cherenkov Telescopes (e.g., Acero et al. 2013).

The PWN broadband spectrum depends both on the particle spectrum that was initially injected at the termination shock by

the pulsar and how it was altered throughout the evolution of the PWN inside its supernova remnant (SNR; Reynolds & Chevalier 1984; Gelfand et al. 2009). The evolution of a PWN is connected to the evolution of the central pulsar, host SNR, and the structure of the surrounding interstellar medium (ISM; Gaensler & Slane 2006). Eventually, the relativistic particle population will be injected into the ISM to contribute to the Galactic cosmic ray (CR) electron–positron population. Current predictions from semianalytical models all indicate this occurs in the late-phase of PWN evolution (Blondin et al. 2001; Gelfand et al. 2009; Temim et al. 2015). Several detailed evolutionary models (e.g., Blondin et al. 2001; van der Swaluw et al. 2004; Gelfand et al. 2009) predict the spatial emission properties and dynamical and radiative evolution for PWNe inside a non-radiative SNR using magnetohydrodynamic and semianalytic modeling methods. There has been significant improvement in model performance when considering the pressure produced by both the particles and the magnetic field of the PWN, which influences the compression evolutionary phase that is initiated by the return of the reverse shock (e.g., Bandiera et al. 2020). Additionally, the PWN and SNR shock structures can be accurately described when accounting for the dependence between the properties of the host SNR and of the PWN (e.g., Bandiera et al. 2021). Other recent works have incorporated multi-dimensional, time-dependent leptonic transport codes that simulate PWN evolution from birth to present day (van Rensburg et al. 2020).

At later evolutionary stages, as is observed with G327.1–1.1, the interaction between the PWN and the reverse shock is a critical stage that initiates rapid synchrotron losses through the

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

Original content from this work may be used under the terms
of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further
distribution of this work must maintain attribution to the author(s) and the title
of the work, journal citation and DOI.

<footnote n="1">¹ Harvard-Smithsonian Center for Astrophysics Cambridge, MA 02138, USA; jeagle@clemson.edu</footnote>

<footnote n="2">² Department of Physics & Astronomy Clemson University Clemson, SC 29634, USA</footnote>

<footnote n="3">³ Department of Astrophysical Sciences Princeton University Princeton, NJ, 08544, USA</footnote>

<footnote n="4">⁴ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="5">⁵ New York University Abu Dhabi, P.O. Box 129188, Abu Dhabi, United Arab Emirates</footnote>

<footnote n="6">⁶ Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Supernova remnant SNR G327.1–1.1 as observed by the Molonglo Observatory Synthesis Telescope at 843 MHz (red), Chandra X-rays between 0.5 and 7 keV (green), and the High Energy Stereoscopic System (H.E.S.S.) TeV and Fermi-LAT best-fit positions indicated by their 95% uncertainty radii, $R_{95} = 0^{\circ}035$ and $R_{95} = 0^{\circ}079$ for the TeV and MeV–GeV positions, respectively. The faint radio SNR shell is outlined as the red-dashed circle.</figcaption></figure>

compression of the PWN. The relativistic particle population is significantly disturbed at this stage. Constraining the physical properties of the particles is difficult, requiring a combination of simulation tools and observable properties over the entire electromagnetic spectrum.

In this paper we report the detection of faint γ-ray emission coincident with the middle-aged SNR G327.1–1.1. In Section 2 we describe the SNR G327.1–1.1 system. A multiwavelength analysis is described in Section 3, including re-analyzing archival Chandra observations to spatially separate the X-ray nebula components and their X-ray spectra in Section 3.2. We present the Fermi-LAT data analysis in Section 3.3. Broadband modeling incorporating the spatially separated multiwavelength data and the resulting best-fit spectral energy distribution (SED) model is described in Section 4. We discuss implications of observations and modeling in Section 5 and present our final conclusions in Section 6.

## 2. SNR G327.1–1.1

SNR G327.1–1.1 is a relatively old remnant with an estimated age $\tau \sim 18,000$ yr based on the Sedov–Taylor solution for an evolved SNR, and a distance of 9 kpc (Sun et al. 1999; Temim et al. 2009). The SNR is visible as a faint radio shell at 843 MHz accompanied by a bright, asymmetric PWN; see Figure 1. A similar structure is observed in the X-rays with thermal emission ($kT \sim 0.3$ keV) present from the SNR shell and strong nonthermal emission (photon index $\Gamma_X \sim 2.1$) from the PWN (Sun et al. 1999; Bocchino & Bandiera 2003; Temim et al. 2009, 2015). Radio and X-ray observations suggest the reverse shock, generated from the pressurized ISM gas swept up from the SNR forward shock, has recently interacted with the PWN asymmetrically, compressing and displacing the PWN to the east. As displayed in Figure 1, the radio nebula is displaced from the X-ray

emission while the X-ray emission has a cometary morphology, suggesting the pulsar is the point source embedded within the X-ray emission peak. Temim et al. (2009) searched for X-ray pulsations within the Chandra data, but the PWN outshines the pulsar, accounting for 94% of the Chandra flux for the compact region of the X-ray nebula. The pulsar is likely moving at a high velocity ($\sim 400$ km s$^{-1}$) to the north based on its location with respect to the geometric center of the SNR. The pulsar’s velocity in combination with an asymmetric reverse shock interaction could explain the southeastern displacement of the PWN (Temim et al. 2009, 2015). This scenario is consistent with the morphology, since the X-ray cometary structure is pointed northwest and elongates to the southeast where it meets the radio relic PWN which is displaced to the east.

It is possible that the formation of the PWN is from the pulsar’s high velocity, creating a bow shock with a comet-shaped X-ray nebula. However, the presence of the X-ray prongs ahead of the pulsar, in addition to evidence for a torus-like structure embedded inside the cometary nebula, suggests the reverse shock is responsible for the observed morphology.

Acero et al. (2011) reported the first detection of point-like TeV γ-ray emission in the direction of SNR G327.1–1.1 observed by the High Energy Stereoscopic System (H.E.S.S.). The TeV source, HESS J1554–550, is uniquely positioned inside the radio SNR shell and overlaps with the radio and X-ray PWN locations, with an uncertainty radius $R < 0^{\circ}035$ (see Figure 1), confirming the high-energy nature of the PWN. MeV–GeV γ-ray emission from the PWN was recently detected by the Fermi-LAT (Xiang et al. 2021).

## 3. Multiwavelength Information

### 3.1. Radio

SNR G327.1–1.1 is well-known for its composite radio structure which was first revealed in the Molonglo Observatory Synthesis Telescope radio survey (Whiteoak & Green 1996). The SNR exhibits a faint non-thermal radio shell and a bright core of emission with a bulge extending out of the core in the northern direction, sometimes referred to as the “finger.” The bulge of radio emission also outlines the X-ray cometary structure and presumed pulsar, supporting the idea that the SNR reverse shock has crushed the PWN and the pulsar has begun exiting the nebula.

More recently, Ma et al. (2016) reported Australia Telescope Compact Array radio observations at 6 and 3 cm that resolve the complex structure for the PWN into two distinct components: a radio “body” which is the relic PWN and a radio “head” which is the radio “finger” structure that overlaps with the compact X-ray nebula (see Figure 2 and Table 2 of Ma et al. 2016). The radio “body” or relic component dominates the overall radio spectrum, which is consistent with the picture that the radio relic contains the oldest particles, whereas the radio “head” represents a younger nebula with fainter radio emission, but brighter X-ray emission. The total radio spectrum is typical for other observed PWNe with a spectral index $\alpha = -0.3 \pm 0.1$. Additionally, the PWN is found to be highly linearly polarized with polarization fractions between 15% and 20% at 6 and 3 cm. The implied magnetic field direction is aligned along the “finger,” with at least a partially azimuthal magnetic field structure in the relic nebula.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Unsmoothed Chandra ACIS-I image of the pulsar wind nebula (PWN) region (2′′6 across) with the “body” and “head” regions used for the spectral extraction, which are chosen to be the same regions of the radio relic (“body”) and radio finger (“head”) spectral extraction regions used in Ma et al. (2016; see their Figure 2, bottom right panel).</figcaption></figure>

## 3.2. X-ray

SNR G327.1–1.1 has been studied in detail in the X-ray band by the Einstein Observatory, ROSAT, ASCA, BeppoSAX, XMM-Newton, and Chandra (see, e.g., Lamb & Markert 1981; Seward et al. 1996; Sun et al. 1999; Bocchino & Bandiera 2003; Temim et al. 2009, 2015). In Figures 1 and 2, the PWN/SNR in X-rays as seen by Chandra is displayed between energies 0.5 and 7 keV. The X-ray morphology consists of the diffuse SNR thermal emission which overlaps with the fainter SNR shape in radio. A slender nonthermal component in the X-ray, which overlaps well with the radio “finger,” corresponds to the PWN and the presumed pulsar, which is most likely the point source embedded within the bright X-ray nebula. Two prong-like structures extend from the northwest of the X-ray peak. It has been suggested that these structures result from an interaction between the PWN and re-heated SN ejecta, though no thermal emission has been detected (Temim et al. 2009, 2015). While the origin of the X-ray prongs remains unclear, it could alternatively point to enhancements in the magnetic field (Temim et al. 2015). In the next section, we extract and measure the X-ray spectra corresponding to the radio “body” or relic and radio “head” or finger components.

### 3.2.1. Chandra Data Analysis

We re-analyze deep archival X-ray observations of SNR G327.1–1.1 performed with the Advanced CCD Imaging Spectrometer, ACIS-I, on board the Chandra X-ray Observatory for a total exposure time of 350 ks. The standard data reduction and cleaning were performed using the Chandra Interactive Analysis of Observations (CIAO v.4.12; Fruscione et al. 2006) software package, resulting in a total clean exposure time of 337.5 ks. An identical data reduction procedure has already been performed and reported (see Temim et al. 2015 for details). We re-analyze the X-ray data incorporating two different source extraction regions, following

<table><caption>**Table 1**
Summary of the 90% CL</caption><thead><tr><th>Component</th><th>Parameter</th><th>Best-fit Value</th></tr></thead><tbody><tr><td>Body or relic</td><td>Reduced χ<sup>2</sup></td><td>1.06</td></tr><tr><td>tbabs</td><td>n<sub>H</sub>(10<sup>22</sup>) cm<sup>-2</sup></td><td>2.30<sup>+0.10</sup><sub>-0.10</sub></td></tr><tr><td>powlaw (PWN)</td><td>Γ</td><td>2.43<sup>+0.04</sup><sub>-0.05</sub></td></tr><tr><td rowspan="3">vnei (SNR)</td><td>Amplitude</td><td>1.60<sup>+0.10</sup><sub>-0.20</sub> × 10<sup>-3</sup></td></tr><tr><td>kT(keV)</td><td>0.22<sup>+0.04</sup><sub>-0.02</sub></td></tr><tr><td>τ (10<sup>11</sup> cm<sup>-3</sup>s)</td><td>2.00<sup>+2.00</sup><sub>-1.00</sub></td></tr><tr><td></td><td>Normalization</td><td>0.07<sup>+0.08</sup><sub>-0.02</sub></td></tr></tbody></table><table><thead><tr><th>Component</th><th>Parameter</th><th>Best-fit Value</th></tr></thead><tbody><tr><td>Head or finger</td><td>Reduced χ<sup>2</sup></td><td>1.10</td></tr><tr><td>tbabs</td><td>n<sub>H</sub>(10<sup>22</sup>) cm<sup>-2</sup></td><td>1.90<sup>+0.10</sup><sub>-0.10</sub></td></tr><tr><td>powlaw (PWN)</td><td>Γ</td><td>2.09<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td rowspan="3">vnei (SNR)</td><td>Amplitude</td><td>2.50<sup>+0.20</sup><sub>-0.10</sub> × 10<sup>-3</sup></td></tr><tr><td>kT(keV)</td><td>0.27<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>τ (10<sup>11</sup> cm<sup>-3</sup>s)</td><td>1.00<sup>+1.00</sup><sub>-0.60</sub></td></tr><tr><td></td><td>Normalization</td><td>0.04<sup>+0.05</sup><sub>-0.01</sub></td></tr></tbody></table>

**Note.** Statistics and parameters for the best-fit model of each component measured in the X-ray analysis.

those used to measure the distinct radio “finger” and “body” PWN components from Ma et al. (2016). Figure 2 shows the regions used for the spectral extraction. A spectrum for each component is extracted using the `specextract` tool in CIAO and modeled using SHERPA within CIAO (Freeman et al. 2001). We measure two X-ray source spectra corresponding to the spatially separated radio components and their radio spectra.

### 3.2.2. Chandra Data Analysis Results

In order to study the properties of the two distinct structures observed for the PWN, we fitted spectra from two different source regions, shown in Figure 2. Since the SNR shell covers most of the field of view, we used the ACIS blank-sky background files adopted to our observation<fnref n="7" /> to extract background spectra from each extraction region. We used the high-energy data from 10 to 12 keV to compare the particle background in each ACIS-I chip in our observations to the particle background in the blank-sky data. This procedure is identical to the one performed and outlined in Temim et al. (2015). To fit the source spectra from each extraction region, we first subtracted the corresponding blank-sky spectrum. Both regions were fitted by a power-law model describing the synchrotron emission from the PWN and a non-equilibrium ionization thermal model `xsvnei` with solar abundances of Wilms et al. (2000) to characterize any thermal emission from the SNR. The results for each spectral region are fully consistent with the findings of prior works (e.g., Temim et al. 2015, see their Table 2 for the “diffuse” and “relic” regions which roughly correspond to the regions employed here). The results of the fits are listed in Table 1. The spectra and best-fit models for the regions are shown in Figure 3.

It is apparent from the softer X-ray spectrum that the body is entering the synchrotron cutoff regime, but the X-ray spectrum observed from the “head” is much harder. This distinction

<footnote n="7">⁷ http://cxc.harvard.edu/ciao/threads/akisbackground/</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Left: data and best-fit model for the “body” PWN component illustrated in Figure 2. Right: data and best-fit model for the “head” PWN component illustrated in Figure 2. In both panels, the nonthermal spectrum originating from the PWN is plotted in green and the thermal spectrum originating from the SNR is plotted in gray. The total X-ray spectral model is plotted in black.</figcaption></figure>

indicates that the two particle populations are under different external conditions, such as different magnetic field strengths and structures that ultimately change the evolutionary history of each population. The average X-ray spectrum (Temim et al. 2015) for the entire PWN demonstrates that there are important contributions from both populations to the total X-ray spectrum. We explore the particle properties from the spatially separated radio and X-ray spectra in Section 4.

### 3.3. γ-Ray

#### 3.3.1. Fermi-LAT

The Fermi Gamma-ray Space Telescope houses both the Gamma-ray Burst Monitor (Meegan et al. 2009) and the Large Area Telescope (LAT; Atwood et al. 2009). The LAT instrument is sensitive to γ-rays with energies 50 MeV to 2 TeV (Atwood et al. 2009) and has been continuously surveying the entire sky every 3 hr since it began operation in 2008 August. The latest upgrade to the event reconstruction process and instrumental response functions (dubbed Pass 8; Atwood et al. 2013) has improved the effective area, point-spread function (PSF), and background rejection capabilities. These upgrades have enabled improvements for the effective energy and angular resolution. The PSF is energy-dependent, becoming large at lower energies but improving at increasing energies. The PSF reaches optimal resolution at $\gtrsim 10$ GeV with a PSF $\sim 0.1$ (68% containment). A major advantage to the Pass 8 event reconstruction upgrade is the classification of detected photon events based on the quality of the angular reconstruction. The events are divided into four categories: PSF0, PSF1, PSF2, or PSF3 types. PSF0 has the worst reconstruction quality and each type increases in reconstruction quality such that PSF3 contains the highest-quality events. For the analysis described in this report, we consider all events.

#### 3.3.2. Fermi-LAT Data Analysis

We re-analyze the Fermi-LAT data using an optimized technique in order to maximize the source detection and to most accurately measure the MeV–GeV γ-ray spectrum. We perform a binned likelihood analysis using 11.5 yr (from 2008 August to 2020 January) of P8R3_SOURCE_V3 photons with

energy between 300 MeV and 2 TeV from all events. We utilize the latest FermiTools package<fnref n="8" /> (v.2.0.8) and FermiPy Python 3 package (v.1.0.1, Wood et al. 2017) to perform data reduction and analysis. We organize the events by PSF type, as described above, using `evtype=4,8,16,32` to represent PSF0, PSF1, PSF2, and PSF3 components. A binned likelihood analysis is performed on each event type and then combined into a global likelihood function for the region of interest (ROI) to represent all events.<fnref n="9" /> Photons detected at zenith angles larger than 100° were excluded to limit the contamination from γ-rays generated by CR interactions in the upper layers of Earth’s atmosphere. The data were additionally filtered to remove time periods when the instrument was not online (e.g., when flying through the South Atlantic Anomaly). The γ-ray data are modeled with the latest comprehensive Fermi-LAT source catalog, 4FGL-DR2 (Ballet et al. 2020), the LAT extended source template archive for the 4FGL catalog,<fnref n="10" /> and the latest Galactic diffuse and isotropic diffuse templates (gll_iem_v07.fits and `iso_P8R3_SOURCE_V3_v1.txt`, respectively).<fnref n="11" />

We fit the 10° ROI with 4FGL point sources and extended sources that are within 15° of the ROI center along with the diffuse Galactic and isotropic emission backgrounds. We allow the background components and sources with TS $\ge 25$ and a distance from the ROI center (chosen to be the TeV PWN position) $\le 3.0°$ to vary in spectral index and normalization.

We computed a series of diagnostic test statistic (TS) and count maps in order to search for and understand any residual γ-ray emission. The TS is defined to be the natural logarithm of the difference in the likelihood of one hypothesis (e.g., presence of one additional source) and the likelihood for the null hypothesis (e.g., absence of source):

$$TS = 2 \times \log\left(\frac{\mathcal{L}_1}{\mathcal{L}_0}\right). \qquad (1)$$

<footnote n="8">⁸ https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/</footnote>

<footnote n="9">⁹ See FermiPy documentation for details: https://fermipy.readthedocs.io/en/0.6.8/config.html.</footnote>

<footnote n="10">¹⁰ https://fermi.gsfc.nasa.gov/ssc/data/access/lat/10yr_catalog</footnote>

<footnote n="11">¹¹ LAT background models and appropriate instrument response functions: https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Left: 2° × 2° test statistic (TS) map of all events from 300 MeV to 2 TeV. Right: 2° × 2° 1–10 GeV TS map of all events. The maximum TS at the PWN/ SNR position is ~18 and ~23 in the 300 MeV–2TeV and 1–10 GeV energy ranges, respectively. The source model used only excludes the point-source modeling emission at the location of G327.1–1.1. The 843 MHz (white) radio contours and H.E.S.S. (blue) and Fermi-LAT (black) 95% uncertainty regions are plotted in both panels.</figcaption></figure>

The TS quantifies how significantly a source is detected at a given location for a specified set of spectral parameters and the significance of such a detection can be estimated by taking the square root of the TS value for one degree of freedom (DOF; Mattox et al. 1996). TS values >25 correspond to a detection significance >4σ for four DOFs.

We generated the count and TS maps in the following energy ranges: 300 MeV–2 TeV, 1–10 GeV, 10–100 GeV, and 100 GeV–2 TeV. The motivation for increasing energy cuts stems from the improving PSF of the Fermi-LAT instrument with increasing energies.<fnref n="12" /> We inspect the maps for additional sources, promptly finding a faint point-like γ-ray source coincident with SNR G327.1–1.1 in a relatively uncrowded region with no known 4FGL counterpart nearby.<fnref n="13" />

TS maps between energies 300 MeV–2 TeV and 1–10 GeV are shown in Figure 4 where significant residual emission coincides with the position and size of SNR G327.1–1.1, with a maximum TS value ~18 in the 300 MeV–2 TeV energy range and TS ~23 between 1 and 10 GeV. It is apparent the residual emission is probably point-like, with little evidence for extension beyond the PSF of the Fermi-LAT.

### 3.3.3. Fermi-LAT Data Analysis Results

To model the γ-ray emission coincident to the PWN inside SNR G327.1–1.1, we add a point source to the 300 MeV–2 TeV source model in addition to known 4FGL sources and the background emission templates for the diffuse Galactic and isotropic components. We first fit the point source with a

power-law spectrum,<fnref n="14" />

$$ \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{-\Gamma} \qquad (2) $$

and photon index Γ = 2. We localize the point source using GTAnalysis localize to find the best-fit position is at (R.A., decl.) = (238°.595, -55°.110) (J2000), an offset of 0°.06 from the exact PWN location, and a 95% positional uncertainty radius r = 0°.079, which overlaps well with the radio and X-ray positions of the PWN. With the new position, we then allow the spectral index and normalization to vary. The TS for a localized point source with a power-law spectrum and photon index Γ = 2.45 ± 0.13 is 29.05.

We investigate the spectral properties of the γ-ray emission by changing the spectrum to a log parabola representation:

$$ \frac{dN}{dE} = N_0 \left( \frac{E}{E_b} \right)^{-(\alpha + \beta \log E / E_b)} \qquad (3) $$

We repeat the same method a third time, replacing the log parabola point source with one described by a power law with a super-exponential cutoff (PLEC),

$$ \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{-\Gamma} \exp(-aE^{\Gamma_2}) \qquad (4) $$

See Table 2 for a summary of the best-fit parameters for each point-source test. Given that all three spectral models provide similar statistical fits, we model the emission with the power-law spectrum since it is the simplest and has the fewest degrees of freedom. As shown in the γ-ray SED in Figure 5, almost all of the emission is observed between 1 and 10 GeV.

<footnote n="12">¹² See https://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm for a review on the dependence of PSF with energy for Pass 8 data.</footnote>

<footnote n="13">¹³ The closest 4FGL sources are more than 1° away in 4FGL-DR2.</footnote>

<footnote n="14">¹⁴ For a review of Fermi-LAT source spectral models see https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html.</footnote>
<table><caption>**Table 2**
Summary of Best-fit Parameters and Associated Statistics with $1\sigma$ Statistical Uncertainties for Each of the Fermi-LAT Point Source Tests in the 300 MeV–2 TeV Energy Band</caption><thead><tr><th>Spectral Model</th><th>log L</th><th>Γ</th><th>α or Γ<sub>1</sub></th><th>β or Γ<sub>2</sub></th><th>N<sub>0</sub> (MeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>E<sub>b</sub> (GeV) or a</th><th>TS</th></tr></thead><tbody><tr><td>Power law</td><td>2143243.14</td><td>2.45 ± 0.13</td><td>...</td><td>...</td><td>(8.12 ± 1.70) × 10<sup>-13</sup></td><td>...</td><td>29.05</td></tr><tr><td>Log parabola</td><td>2143245.52</td><td>...</td><td>2.29 ± 0.19</td><td>0.20 ± 2.85 × 10<sup>-4</sup></td><td>(1.86 ± 0.35) × 10<sup>-13</sup></td><td>2.00</td><td>32.95</td></tr><tr><td>Power law with super-exp. cutoff</td><td>2143245.54</td><td>...</td><td>1.04 ± 1.58</td><td>0.67 (fixed)</td><td>(2.19 ± 2.88) × 10<sup>-12</sup></td><td>0.01 ± 2.79 × 10<sup>-5</sup></td><td>32.97</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Best-fit γ-ray spectral energy distribution between 300 MeV and 1 TeV over seven equally spaced logarithmic energy bins with $1\sigma$ statistical uncertainties (red) and the quadratic sum of systematic errors (black; see Section 3.3.4). TS values for each spectral bin are plotted as the green histogram. The best-fit power-law spectrum is fit with a photon index $\Gamma_{\gamma} = 2.45 \pm 0.13$.</figcaption></figure>

We run extension tests for the best-fit point source in FermiPy utilizing GTAnalysis.extension and the two spatial templates supported in the FermiPy framework, the radial disk and radial Gaussian. Both of these extended templates assume a symmetric 2D shape with width parameters radius and sigma, respectively. We allow the position to vary when finding the best-fit spatial extension for both templates. The summary of the best-fit parameters for the extended templates are listed in Table 3. The γ-ray emission shows no evidence for extension.

The origin of the emission is unlikely to be the SNR, given the distinct location inside the SNR shell in addition to a general lack of evidence for an energetic interaction between the SNR and the ISM that would produce γ-rays. The 300 MeV–2 TeV energy flux $F_{\gamma} = 5.12(\pm 0.97) \times 10^{-12}$ ergs cm$^{-2}$ s$^{-1}$ gives a γ-ray luminosity $L_{\gamma} = 4.94 \times 10^{34} (\frac{d}{9kpc})^2$ ergs s$^{-1}$. Taking the predicted current spin-down power $\dot{E} \sim 3.1 \times 10^{36}$ ergs s$^{-1}$ of the pulsar from Temim et al. (2015), we measure the γ-ray conversion efficiency $\eta_{\gamma} \sim 0.016$, which is comparable to other GeV PWNe (Acero et al. 2013). However, we cannot exclude a pulsar contribution, as several γ-ray pulsars have similar values of $\eta$ (Abdo et al. 2013). The γ-ray emission is not bright enough to perform a reliable pulsation search and we cannot rule out the scenario from the γ-ray SED alone. The results reported here are consistent with those reported in Xiang et al. (2021), where the coincident γ-ray emission is detected just below the Fermi–LAT threshold of $4\sigma$, $TS = 22.94$ for four DOFs and characterized best with photon index $\Gamma = 2.35 \pm 0.24$.

### 3.3.4. Systematic Error

We account for systematic uncertainties introduced by the choice of the interstellar emission model (IEM) and the

<table><caption>**Table 3**
Summary of Best-fit Parameters and Associated Statistics for Each Spatial Template Tested in the 300 MeV–2 TeV Energy Band</caption><thead><tr><th>Spatial<br>Template</th><th>TS</th><th>TS<sub>ext</sub></th><th>R.A., Decl.<br>(J2000)</th><th>Best-fit<br>Radius or<br>Sigma (°)</th><th>95% Radius<br>Upper<br>Limit (°)</th></tr></thead><tbody><tr><td>Point source</td><td>29.05</td><td>...</td><td>238.59,<br>-55.11</td><td>...</td><td>...</td></tr><tr><td>Radial disk</td><td>29.29</td><td>0.25</td><td>238.65,<br>-55.07</td><td>0.004</td><td>0.083</td></tr><tr><td>Radial<br>Gaussian</td><td>29.30</td><td>0.26</td><td>238.65,<br>-55.07</td><td>0.003</td><td>0.077</td></tr></tbody></table>

effective area, which mainly affect the spectrum of the faint γ-ray emission. We have followed the prescription developed in de Palma et al. (2013) and Acero et al. (2016), which generated eight alternative IEMs using a different approach than the standard IEM (see Acero et al. 2016 for details). For this analysis, we employ the eight alternative IEMs (aIEMs) that were generated for use on Pass 8 data in the Fermi-LAT Galactic Extended Source Catalog (Ackermann et al. 2017). The γ-ray point source coincident with PWN G327.1–1.1 is refit with each aIEM to obtain a set of eight values for the spectral flux that we compare to the standard model following Equation (5) in Acero et al. (2016). We estimate the systematic uncertainties introduced by the instrument response function (IRF) uncertainties<fnref n="15" /> while enabling energy dispersion as follows: ±3% for $E < 100$ GeV, ±4.5% for $E = 175$ GeV, and ±8% for $E = 556$ GeV. Since the IEM and IRF systematic errors are taken to be independent, we can evaluate both and perform the quadratic sum for the total systematic error. We find that the systematic errors are most important in the first two energy bins of the γ-ray spectrum, which effectively convert the first data point to an upper limit, but are negligible in the remaining energy bins. The corresponding quadratic sum of the systematic errors in the first four energy bins is plotted in Figure 5 in black.

## 4. Broadband Modeling

A broadband model was derived for G327.1–1.1 from a semianalytic simulation that predicts the evolution of a PWN inside a nonradiative SNR and was reported in Temim et al. (2015). Assuming a broken power law (BPL) particle injection spectrum and given a set of input parameters, the model simultaneously predicts the PWN and SNR size, the PWN magnetic field, and the PWN and SNR expansion velocities as a function of age. The final particle and photon spectra are measured by evolving the BPL injection spectrum with time while accounting for the evolution of the PWN and host SNR. The particle spectrum evolves until it corresponds to the age of

<footnote n="15">¹⁵ https://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Best-fit broadband model fits for emission from the evolved PWN in G327.1–1.1 along with radio observations of the PWN (Ma et al. 2016), X-ray fluxes, and spectral indices (this work), new Fermi-LAT data (this work), and H.E.S.S. γ-ray emission (Acero et al. 2011). For comparison, we plot the broadband representation derived from the evolved broken power-law injection spectrum in Temim et al. (2015) in solid blue.</figcaption></figure>

the SNR at which the simulated PWN and SNR radii match the observed sizes (see Gelfand et al. 2009; Temim et al. 2015 for details). The corresponding broadband photon spectrum can reproduce well the radio, X-ray, and γ-ray data of G327.1–1.1 considered by Temim et al. (2015), and suggests the remnant is old, $\tau \sim 18,000$ yr, with distinct particle components apparent in the simulated SED (see Figure 6).

The resolved PWN morphology in radio and X-ray gives a unique opportunity to study the individual properties of the particle populations through their distinct broadband spectral features, first investigated by Temim et al. (2015). In this section we attempt to model the particle populations independently in an effort to characterize each component that contributes to the observed broadband emission, incorporating the new Fermi-LAT data, and compare to the model presented in Temim et al. (2015). We use the NAIMA Python package (Zabalza 2015), designed to compute the non-thermal radiation from relativistic particle populations. NAIMA is a simple radiation code, in contrast to the semianalytic method employed by Temim et al. (2015), which uses a constant particle distribution and hence it does not account for the evolution of the total particle energy and spectrum as the PWN and host SNR evolve. To explore the possible particle distributions, we assume a power law shape with a PLEC,

$$f(E) = A \left(\frac{E}{E_0}\right)^{-\Gamma} \exp\left(-\frac{E}{E_c}\right) \quad (5)$$

where A is in eV⁻¹. We then input similar values as those estimated in Temim et al. (2015), such as the particle index and energy break $E_{break}$ in the particle spectrum, and also consider the same three photon fields for the IC scattering (ICS) component. The photon fields consist of the cosmic microwave background (CMB) with energy density $\Omega_{CMB} = 0.261$ eV cm⁻³ and temperature of 2.72 K, and two warmer photon fields from heated dust and nearby starlight, with energy densities 1.5 $\Omega_{CMB}$ and 2.0 $\Omega_{CMB}$ and temperatures $T = 3.5 \times 10^3$ K and $T = 50$ K, respectively. We then test the combination of free

PLEC parameters that can best explain the broadband spectrum.

While a single particle distribution can generally reproduce the observed synchrotron and IC spectral shapes, it cannot reproduce the radio spectrum nor explain the MeV–GeV γ-rays adequately. The sum of two leptonic populations under the same ambient photon fields thus represents the broadband spectrum of PWN G327.1–1.1: one component belonging to the older, larger and displaced radio relic and the second component belonging to the compact, younger X-ray nebula. The observed displacement of the radio and X-ray structures, in addition to the estimated age for the nebula as a whole, strongly indicate the presence of distinct particle components under different physical conditions as a result of the asymmetric crushing of the PWN. As such, we fit the magnetic field associated with each leptonic population independently. This assumption is further supported by the distinct radio and X-ray spectra measured for each population. As shown in Figure 6, the lower-energy population clearly dominates the radio spectrum for the nebula, but is in the synchrotron cut-off regime in the X-ray band, while the higher-energy population is faint in radio but bright in X-ray. Observationally, the larger radio to γ-ray ratio in the relic implies that the magnetic field is larger than in the X-ray head.

Table 4 lists the best-fit parameters for the two-leptonic PLEC particle spectra investigated here along with the best-fit parameters from the evolved BPL particle spectrum in Temim et al. (2015) for comparison. The corresponding best-fit broadband spectrum for the two PLEC lepton populations is displayed in Figure 6. We find the low-energy (i.e., older) population is best characterized by a particle index before the energy cutoff $\Gamma = 1.61$. We find the high-energy (i.e., younger) population is best-fit with a particle index before the cutoff $\Gamma = 2.15$. The implied magnetic field is much higher for the lower-energy population $B = 101 \mu$ G, about 10 times higher than the higher-energy population $B = 10.8 \mu$ G. The low-energy cutoff is $E_c = 0.99$ TeV and the high-energy cutoff is $E_c = 41.69$ TeV.
<table><caption>**Table 4**
Summary of Best-fit Parameters of the Broadband Models in Figure 6 with Those of Temim et al. (2015) for Comparison</caption><thead><tr><th>Particle Spectrum</th><th>Γ<sub>1</sub></th><th>Γ<sub>2</sub></th><th>E<sub>c</sub> or E<sub>break</sub> (TeV)</th><th>B (μG)</th></tr></thead><tbody><tr><td>PLEC1 (low-energy pop.)</td><td>1.61</td><td>...</td><td>0.99</td><td>101.00</td></tr><tr><td>PLEC2 (high-energy pop.)</td><td>...</td><td>2.15</td><td>41.69</td><td>10.76</td></tr><tr><td>BPL (Temim et al. 2015)</td><td>1.48</td><td>2.20</td><td>0.31</td><td>10.76</td></tr></tbody></table>

The low- and high-energy particle indices reported here are similar to the values 1.48 and 2.20 found in Temim et al. (2015) which were motivated by the observed average radio and X-ray indices. Because the low-energy population clearly dominates the radio emission, the comparable low-energy index is not surprising. The X-ray emission is a combination of the two populations as shown in our Figure 6 such that the similarity for the high-energy index found here and in Temim et al. (2015) is expected. Moreover, the X-ray spectrum of each particle component shows two important distinctions: the older population is entering the synchrotron cutoff regime while the younger population does not show strong evidence for synchrotron cutoff in the 0.5–6 keV X-ray band. This suggests the highest-energy particles are not suffering the same synchrotron losses as the lowest-energy particles and is supported by the high magnetic field value inferred from the lower-energy population. The high-energy magnetic field component is much lower, but matches the magnetic field strength inferred from Temim et al. (2015).

It is clear that the X-ray spectrum from the relic is underpredicted by the low-energy population in the broadband model. This could be explained by the assumption that both particle populations are characterized by a power law with an exponential cutoff distribution. The particle spectrum of an evolved PWN such as G327.1–1.1 has been shown to deviate beyond a simple power-law distribution even if the initial injection spectrum is a power law (Gelfand et al. 2009), such that if the low-energy population truly represents the relic nebula, the particle spectrum may not be best described in power-law form. This limitation is not present in the Temim et al. (2015) model, however, since the particle spectrum evolution is considered. The synchrotron emission peaks at a flux about an order of magnitude higher than what is predicted in the Temim et al. (2015) model and is shifted to lower energy by roughly one order of magnitude. This is due to the much stronger magnetic field strength combined with the low-energy cutoff that characterizes the oldest particle population. The last most striking difference is that the Temim et al. (2015) model can accurately predict the MeV–GeV γ-ray spectrum, intersecting with all of the most constraining data points in the Fermi-LAT SED. The IC emission in the Temim et al. (2015) model has a two-peak spectral shape that manifests from the evolutionary history of the particles throughout the PWN’s lifetime, which accounts for the total energy input of the PWN decreasing with time once the pulsar exits the nebula.

The divergences in the model fits, especially for the oldest particles, demonstrate the importance in considering the evolutionary history for broadband studies of PWNe. Despite the use of simplified radiative models to develop a new broadband representation, both models depict a physical scenario where the reverse shock has compressed the relic PWN while the exiting pulsar generates a new X-ray nebula which may have minimal to no impact from the reverse shock

and subsequent compression. We discuss the physical implications of the models in the following section.

## 5. Discussion

We have presented a broadband model for the PWN emission in G327.1–1.1 that is characterized by the sum of two-lepton PLEC spectra and compared to the evolved particle spectrum derived in Temim et al. (2015). The best-fit model estimated from the evolved particle spectrum is based on the dynamical model for PWN evolution inside a non-radiative SNR, first described in Gelfand et al. (2009). We further investigate the physical implications of the presented model and of the Temim et al. (2015) model by comparing the predicted observable properties from the evolutionary phases expected for PWNe from Gelfand et al. (2009) with those observed from the PWN inside SNR G327.1–1.1.

The initial expansion stage is the first phase of PWN evolution and occurs from the PWN pressure being much greater than that of the surrounding SN ejecta. Throughout the PWN expansion, the pulsar radiates more energy than the PWN loses from adiabatic expansion and synchrotron losses. The expansion phase introduces a rapid decline in the PWN magnetic field, which results in a steep decline in synchrotron luminosity with time. As a result, adiabatic losses will dominate over radiative losses until the reverse shock crushes the PWN, which also marks the end of the first expansion stage. The particle spectrum is dominated by previously injected particles at all energies for the majority of the PWN expansion. Similarly, the corresponding photon spectrum is dominated by previously injected particles at all energies, forming a single peak structure for both synchrotron and IC emission. The initial expansion stage occurs within the first ~1–5 kyr of the PWN lifetime, but depends on the ambient ISM density and the host SNR evolution. The initial expansion phase ends when the SNR reverse shock first collides with the PWN.

The particle spectrum and consequent photon spectrum during the first contraction phase will see dramatic changes once collision with the reverse shock has begun. The synchrotron lifetime of the highest-energy particles will become significantly less than the PWN age and will introduce an energy break in the particle spectrum. The most recently injected particles dominate beyond this energy break. During the PWN compression, the synchrotron luminosity will increase due to the strengthening magnetic field. Meanwhile, the IC luminosity depends mostly on the highest-energy particles, which are those most recently injected by the central pulsar. The corresponding photon spectrum sees the synchrotron emission peak decrease in energy due to the overall energy decreasing with time from synchrotron losses. The IC emission peak is also expected to decline in the first contraction stage both in part due to the central pulsar input and to the short synchrotron lifetime of the particles, such that only the highest-energy particles contribute to the TeV emission.

After the PWN undergoes significant compression, the maximum particle energy and magnetic field strength will eventually begin to decrease. The corresponding photon spectrum will see another shift in the synchrotron and IC peaks such that the majority of the broadband emission occurs in the radio and soft and hard X-ray bands. The continuous injection of high-energy particles by the central pulsar will result in two distinct components visible in both the particle and photon spectra: a lower-energy population composed of
particles injected at earlier times and a higher-energy comp-
onent composed of recently injected particles. The corresp-
onding photon spectrum features the double-peaked structure in
both synchrotron and IC emission processes, akin to the Temim
et al. (2015) model for G327.1–1.1 and similar to the one
presented in Section 4. At this stage, the radio emission is
attributed to synchrotron radiation from the oldest particles and
the GeV emission from the same population ICS against
ambient photons, while the rest of the emission is dominated by
recently injected particles. In the case where the pulsar exits the
nebula and does not re-enter, the pulsar still continues to inject
high-energy particles into its surroundings, generating a new
high-energy PWN that is separate from the relic PWN it has left
behind.

Much of the depicted evolutionary cycle depends on the
unique physical properties of the system and surroundings such
as the progenitor characteristics, the ISM density, and the
pulsar spin-down. The observed morphology and spectral
characteristics of G327.1–1.1 suggest it to be an evolved
system that has already encountered or is currently encounter-
ing the passage of the reverse shock through the SNR interior,
crushing the PWN in its first contraction.

The evolutionary cycle described here is only valid for one-
dimensional, spherically symmetric, non-radiative host SNR
systems. The spherically symmetric condition introduces a
series of compression and re-expansion cycles to the nebula
that are initiated by the passage of the reverse shock. For
increasingly asymmetric systems such as G327.1–1.1, the
evolutionary cycle after first compression can differ signifi-
cantly (Gelfand et al. 2009; Temim et al. 2015; Bandiera et al.
2020). The thin-shell approximation is also used in the model
presented by Temim et al. (2015), where the PWN is assumed
to be surrounded by a thin shell of swept-up SN ejecta. If this
assumption is not realistic, the impact of the reverse shock and
subsequent compression can be over-predicted (e.g., Bandiera
et al. 2020; Fiori et al. 2022). Despite these limitations,
particularly in characterizing the initial compression stage, the
physical implications from the simplified model we present in
Section 4 and the one reported in Temim et al. (2015) are
consistent with the observational properties for the system.

The size of the SNR (~22 pc) at the estimated distance
~9 kpc and the displaced PWN radio and X-ray counterparts
support the scenario of an older system that has experienced an
interaction with the SNR reverse shock. It is likely the reverse
shock has displaced the PWN to the east and, combined with
the velocity of the presumed pulsar (Temim et al. 2009, 2015),
has enabled the pulsar to exit the PWN leaving behind a radio
relic PWN while continuously injecting high-energy particles
into the surroundings, generating the X-ray nebula in its wake.
The radio and X-ray spectra further validate the scenario
depicted here, where significant synchrotron losses are
apparent, especially for the oldest particles in the relic, and
likely being the direct result of the reverse shock interaction.
The observed morphology also implies that the highest-energy
emission, which is observed close to the putative pulsar, is near
the maximum particle energy, and indicates a steepening in the
particle spectrum (Temim et al. 2015).

6. Conclusions

We have reported the detection of faint γ-ray emission
discovered coincident in location with the evolved
SNR G327.1–1.1. The position and extent is consistent with

an origin from the high-energy PWN which is also known to
emit TeV γ-rays as HESS J1554–550. We perform a multi-
wavelength investigation to explore the broadband models
supported by observations, utilizing a one-dimensional radia-
tive modeling package NAIMA. We compare the results to an
in-depth semianalytic simulation that predicts the evolved
particle injection spectrum based on hydrodynamically derived
quantities such as supernova explosion energy and ejecta mass
and reported in Temim et al. (2015). We investigate the
broadband emission structure based on our understanding of
the PWN evolutionary sequence outlined in Gelfand et al.
(2009). We find that the broadband model presented here and
the one reported in Temim et al. (2015), although obtained
through very different techniques, still provide consistent and
plausible physical implications. The faint γ-ray emission
reported here may have a pulsar contribution, particularly for
energies *E* < 10 GeV. Comparison of the two broadband
models explored here also suggests that particle injection
history is required for accurate γ-ray emission models of an
evolved system like G327.1–1.1. The Cherenkov Telescope
Array<fnref n="16" /> will be ~5–10 times more sensitive than H.E.S.S. at
1 TeV with an angular resolution that could spatially resolve
the PWN TeV morphology to be compared to the radio and
X-ray components, thus revealing how the distinct particle
components contribute to the lower-energy γ-ray emission.

The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K.A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. 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'Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.

*Software:* CIAO (v4.12 Fruscione et al. 2006), FermiPy (v1.0.1 Wood et al. 2017), Fermitools: Fermi Science Tools (v2.0.8 Fermi Science Support Development Team 2019), NAIMA (Zabalza 2015).

ORCID iDs

Jordan Eagle https://orcid.org/0000-0001-9633-3165
Daniel Castro https://orcid.org/0000-0002-0394-3173
Tea Temim https://orcid.org/0000-0001-7380-3144
Jean Ballet https://orcid.org/0000-0002-8784-2977
Patrick Slane https://orcid.org/0000-0002-6986-6756
Joseph Gelfand https://orcid.org/0000-0003-4679-1058
Matthew Kerr https://orcid.org/0000-0002-0893-4073
Marco Ajello https://orcid.org/0000-0002-6584-1703

<footnote n="16">¹⁶ https://www.cta-observatory.org/</footnote>
References

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

Acero, F., Ackermann, M., Ajello, M., et al. 2013, ApJ, 773, 77

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

Acero, F., Djannati-Ataï, A., Förster, A., Gallant, Y., & Renaud, M. 2011, ICRC (Beijing), 7, 185

Ackermann, M., Ajello, M., Baldini, L., et al. 2017, ApJ, 843, 139

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

Atwood, W. B., Baldini, L., Bregeon, J., et al. 2013, ApJ, 774, 76

Ballett, J., Burnett, T. H., Digel, S. W., & Lott, B. 2020, arXiv:2005.11208

Bandiera, R., Bucciantini, N., Martin, J., Olmi, B., & Torres, D. F. 2020, MNRAS, 499, 2051

Bandiera, R., Bucciantini, N., Martin, J., Olmi, B., & Torres, D. F. 2021, MNRAS, 508, 3194

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

Bocchino, F., & Bandiera, R. 2003, A&A, 398, 195

de Palma, F., Brandt, T. J., Johannesson, G., & Tibaldo, L. 2013, arXiv:1304.1395

Fermi Science Support Development Team 2019, Fermitools: Fermi Science Tools, Astrophysics Source Code Library, ascl:1905.011

Fiori, M., Olmi, B., Amato, E., et al. 2022, MNRAS, 511, 1439

Freeman, P., Doe, S., & Siemiginowska, A. 2001, Proc. SPIE, 4477, 76

Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, Proc. SPIE, 6270, 62701V

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

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

Lamb, R. C., & Markert, T. H. 1981, ApJ, 244, 94

Ma, Y. K., Ng, C. Y., Bucciantini, N., et al. 2016, ApJ, 820, 100

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

Meegan, C., Lichti, G., Bhat, P. N., et al. 2009, ApJ, 702, 791

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

Seward, F. D., Kearns, K. E., & Rhode, K. L. 1996, in Roentgenstrahlung from the Universe, ed. H. U. Zimmermann, J. Trümper, & H. Yorke, 277

Slane, P. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin (Cham: Springer), 2159

Sun, M., Wang, Z.-r., & Chen, Y. 1999, ApJ, 511, 274

Temim, T., Slane, P., Gaensler, B. M., Hughes, J. P., & Van Der Swaluw, E. 2009, ApJ, 691, 895

Temim, T., Slane, P., Kolb, C., et al. 2015, ApJ, 808, 100

van der Swaluw, E., Downes, T. P., & Keegan, R. 2004, A&A, 420, 937

van Rensburg, C., Venter, C., Seyffert, A. S., & Harding, A. K. 2020, MNRAS, 492, 3091

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

Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914

Wood, M., Caputo, R., Charles, E., et al. 2017, ICRC (Busan), 301, 824

Xiang, Y., Xing, Y., & Jiang, Z. 2021, ApJ, 912, 117

Zabalza, V. 2015, ICRC (The Hague), 34, 922