OPEN ACCESS

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

CrossMark

# Fermi-LAT Gamma-Ray Emission Discovered from the Composite Supernova Remnant B0453-685 in the Large Magellanic Cloud

Jordan Eagle<fnref n="1" />,<fnref n="2" /> ID, Daniel Castro<fnref n="1" /> ID, Peter Mahhov<fnref n="3" />, Joseph Gelfand<fnref n="3" /> ID, Matthew Kerr<fnref n="4" /> ID, Patrick Slane<fnref n="1" /> ID, Jean Ballet<fnref n="5" /> ID, Fabio Acero<fnref n="5" /> ID, Samayra Straal<fnref n="3" /> ID, and Marco Ajello<fnref n="2" /> ID

Received 2022 October 10; revised 2023 January 30; accepted 2023 February 1; published 2023 March 1

## Abstract

We report the second extragalactic pulsar wind nebula (PWN) to be detected in the megaelectronvolt–gigaelectronvolt band by the Fermi-LAT, located within the Large Magellanic Cloud. The only other known PWN to emit in the Fermi band outside of the Milky Way is N157B, which lies to the west of the newly detected gamma-ray emission at an angular distance of 4°. Faint, pointlike gamma-ray emission is discovered at the location of the composite supernova remnant (SNR) B0453-685 with a ~4σ significance with energies ranging from 300 MeV–2 TeV. We present the Fermi-LAT data analysis of the new gamma-ray source, coupled with a detailed multiwavelength investigation to understand the nature of the observed emission. Combining the observed characteristics of the SNR and the physical implications from broadband modeling, we argue it is unlikely that the SNR is responsible for the gamma-ray emission. While the gamma-ray emission is too faint for a pulsation search, we try to distinguish between any pulsar and PWN component of SNR B0453-685 that could be responsible for the observed gamma-ray emission using semi-analytic models. We determine the most likely scenario is that the old PWN (τ ~ 14,000 yr) within B0453-685 has been impacted by the return of the SNR reverse shock with a possible substantial pulsar component below 5 GeV.

*Unified Astronomy Thesaurus concepts: Gamma-ray astronomy (628); Large Magellanic Cloud (903); X-ray astronomy (1810); Supernova remnants (1667); Pulsar wind nebulae (2215)*

## 1. Introduction

Pulsar wind nebulae (PWNe) are descendants of core-collapse supernovae, each powered by an energetic, rapidly rotating neutron star. As the neutron star spins down, rotational energy is translated into a relativistic particle wind, made up of mostly electrons and positrons (Slane 2017). The evolution of a PWN is connected to the evolution of the central pulsar, host supernova remnant (SNR), and the structure of the surrounding interstellar medium (ISM; Gaensler & Slane 2006). Eventually, the relativistic particle population will be injected into the ISM of the host galaxy and may contribute to the cosmic ray (CR) electron–positron population (Malyshev et al. 2009).

Synchrotron emission from relativistic electrons is observed from the majority of PWNe, from radio wavelengths to hard X-rays. Moreover, CR electrons are expected to scatter off of local photon fields, resulting in inverse Compton (IC) emission at gamma-ray energies (Gaensler & Slane 2006). Accordingly, the majority of PWNe have been discovered in the radio or X-ray bands and an increasing number of discoveries are being made in teraelectronvolt gamma rays. In fact, the majority of the Galactic teraelectronvolt source population is found to be PWNe as observed by Cherenkov Telescopes (~37, e.g., Wakely & Horan 2008; Acero et al. 2013). On the other hand, only 11 PWNe have been firmly identified in the megaelectronvolt–gigaelectronvolt band by Fermi-LAT (Atwood et al. 2009).

However, upgrades in the event processing of the Fermi-LAT data have significantly improved the spatial resolution and sensitivity of the instrument (Pass 8; Atwood et al. 2013). Taking advantage of the upgrade and using Fermi-LAT observations with ~11.5 yr of data, we have discovered a new Fermi-LAT gamma-ray source located in the Large Magellanic Cloud (LMC) that belongs to the composite SNR B0453-685. We combine the new gamma-ray measurements with available multiwavelength data for the region to determine that the PWN is the most likely origin of the gamma rays and that a pulsar may contribute to the lower-energy gamma-ray emission.

The broadband spectrum of a PWN depends on both the particle spectrum that was initially injected by the pulsar and how it was altered throughout the evolution of the PWN inside its SNR (Reynolds & Chevalier 1984; Gelfand et al. 2009). In order to rigorously explore the characteristics of the underlying particle population(s), we present a semi-analytic simulation for the dynamical and radiative evolution of a PWN inside an SNR.

In Section 2, we describe the SNR B0453-685 system. We present a multiwavelength analysis in Section 3, describing the X-ray analysis using archival Chandra observations in Section 3.2 and the Fermi-LAT data analysis in Section 3.3. We present simple broadband models investigating the gamma-ray origin in Section 4.1. We further simulate a broadband spectral model using a semi-analytic model of PWN evolution, which incorporates known properties of the system and we report the resulting best-fit spectral energy distribution (SED) in Section 4.2. We discuss the implications of observations and

<figure><img src="image_2.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; jordan.l.eagle@nasa.gov</footnote>

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

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

<footnote n="4">⁴Space Science Division, US Naval Research Laboratory, Washington, DC 20375-5352, USA</footnote>

<footnote n="5">⁵Université Paris-Saclay & Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France</footnote>
<figure><img src="image_3.png" /></figure>

**Figure 1.** Left: SAO DS9 image of the LMC in the Hα band from the Southern H-Alpha Sky Survey Atlas (SHASSA, supported by the National Science Foundation; Gaustad et al. 2001). The P1–P4 labels identify the four brightest 4FGL point sources in the LMC, following the naming convention used in Ackermann et al. (2016). P1 is the most energetic pulsar ever detected, PSR J0540-6919, which lies <0.5° from P2. P2 is the possible Fermi-LAT PWN N157B. P3 is a high-mass binary system and P4 is the SNR N132D located near the 30 Doradus region. The four extended templates used to describe the diffuse gamma-ray emission from the LMC (components E1–E4 in Ackermann et al. 2016) are indicated with green circles. The location of SNR B0453-685 is marked in white with radius r = 0.05. The two closest known radio pulsars near SNR B0453-685 are labeled as white diamonds. (We used the ATNF radio pulsar catalog https://www.atnf.csiro.au/research/pulsar/psrcat/) Both are located too far away from the SNR to reasonably be a central pulsar candidate. The coordinates are labeled and are in equatorial J2000 degrees throughout the paper unless otherwise noted. Middle: the SAO DS9 image of the 1.4 GHz radio emission observed from SNR B0453-685 (Gaensler et al. 2003). The white contours correspond to the central PWN and the outer SNR shell as observed in X-ray (right panel). Right: tricolor X-ray flux map generated in SAO DS9 of SNR B0453-685 (Gaensler et al. 2003). Red represents soft X-ray emission between 0.5 and 1.2 keV, green represents medium flux between 1.2 and 2 keV, and blue represents hard flux from 2–8 keV. Soft and medium X-ray emission outlines and fills the entire SNR while the hard X-ray emission is heavily concentrated toward the center of the SNR where the PWN is located. (Manchester et al. 2005).

modeling and we provide our final conclusions in Sections 5 and 6.

## 2. SNR B0453-685

SNR B0453-685 is located in the LMC with a distance $d \approx 50$ kpc (Clementini et al. 2003). The LMC has an angular size of nearly 6° in the sky where SNR B0453-685 (angular size $r < 0.05°$) is positioned on the western wall of Hα emission as shown in Figure 1, left panel. SNR B0453-685 was identified by Gaensler et al. (2003) as a middle-aged ($\tau \sim 13$ kyr) composite SNR hosting a bright, polarized central core based on observations at 1.4 and 2.4 GHz frequencies and in X-rays between 0.3 and 8.0 keV; see the middle and right panels of Figure 1. A thin, faint SNR shell is visible in both radio and X-ray (0.3–2.0 keV) with the softer, diffuse X-ray emission filling the SNR. Within the radio SNR shell is a much brighter, large, and polarized, central core: the PWN. The PWN also dominates the hard X-ray emission (2.0–8.0 keV, Figure 1). While the radio and X-ray observations reported by Gaensler et al. (2003) indicated the composite morphology of the SNR, no pulsations from a central pulsar have been detected.

Manchester et al. (2006) performed a deep radio pulsar search in both of the Magellanic Clouds with the Parkes 64 m radio telescope and reported 14 total pulsars, 11 of which were located within the LMC, but none were associated with SNR B0453-685. It is reported in later investigations (e.g., Lopez et al. 2011; McEntaffer et al. 2012) using the same Chandra X-ray observations as those in Gaensler et al. (2003) that an X-ray point source was detected inside the central PWN core using the `wavdetect` tool within the Chandra data reduction software package, Chandra Interactive Analysis of Observations (CIAO; Fruscione et al. 2006). This remains the most promising evidence for the central pulsar.

Displayed in Figure 1, left panel, are the few known sources within the LMC that emit gamma rays in the Fermi-LAT band, labeled P1–P4 following the convention of Ackermann et al. (2016). Only one LMC PWN, N157B (P2), is identified as a gigaelectronvolt (Ballet et al. 2020) and teraelectronvolt (H.E.S.S. Collaboration et al. 2012) gamma-ray source, and it is located on the opposite (eastern) wall of the LMC with respect to SNR B0453-685. N157B is located in a very crowded area, accompanied by two bright gamma-ray sources nearby, SNR N132D and PSR J0540-6919. SNR B0453-685, however, is conveniently located in a much less crowded region of the LMC, making its faint pointlike gamma-ray emission detectable even against the diffuse LMC background, diffuse Galactic foreground, and the isotropic background emissions.

## 3. Multiwavelength Information

### 3.1. Radio

Australia Telescope Compact Array (ATCA) observations at 1.4 and 2.4 GHz revealed the composite nature of SNR B0453-685, indicating the presence of a PWN (Gaensler et al. 2003). The PWN is visible as a bright central core that is surrounded by the SNR shell roughly 2' in diameter. Gaensler et al. (2003) measure the flux density of the radio core to be $46 \pm 2$ mJy at both 1.4 and 2.4 GHz. The PWN radio spectrum is flat, with $\alpha = -0.10 \pm 0.05$ (Gaensler et al. 2003). No central point source such as a pulsar is seen, but the authors place an upper limit on a point source of 3 mJy at 1.4 GHz and 0.4 mJy at 2.4 GHz at the location of the emission peak and suggest the PWN to be powered by a Vela-like pulsar with a spin period of $P \approx 100$ ms, a surface magnetic field $B \approx 3 \times 10^{12}$ Gauss, and a spin-down luminosity of $\dot{E} \approx 10^{37}$ erg s$^{-1}$.

Haberl et al. (2012) observed SNR B0453-685 with ATCA at 4.8 and 8.6 GHz, providing radio flux density measurements
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Tricolor X-ray flux map generated in SAO DS9 of B0453-685. Red denotes the range of 0.5–1.2 keV, green denotes the range of 1.2–2 keV, and blue denotes the range of 2–8 keV. The source and background regions used for spectral analysis are indicated. The yellow dashed circle corresponds to the PWN region, the white annulus corresponds to the SNR region, and the large green annulus excluding six bright X-ray point sources corresponds to the background region.</figcaption></figure>

of both the SNR and PWN. The authors measure a flat radio spectrum for the PWN, with $\alpha_{pwn} = -0.04 \pm 0.04$, along with significant polarization from the PWN core at 1.4, 2.4, 4.8, and 8.6 GHz frequencies. The outer SNR shell, excluding the PWN contribution, has a radio spectral index $\alpha_{shell} = -0.43 \pm 0.01$, which is a typical value for radio SNR shells.

## 3.2. X-Ray

### 3.2.1. Chandra X-Ray Data Analysis

SNR B0453-685 has been analyzed in X-rays in great detail (Gaensler et al. 2003; Lopez et al. 2009; Haberl et al. 2012; McEntaffer et al. 2012) with data from the XMM-Newton and Chandra X-ray telescopes. Thermal X-ray emission dominates the soft X-rays and is largely attributed to the SNR, while the hard X-ray emission is concentrated toward the center of the remnant where the PWN is located (see Figure 1, right). In order to understand the gamma-ray origin, we must combine the new Fermi-LAT data with available multiwavelength data for the region. Therefore, we reanalyzed archival Chandra X-ray observations (ObsID: 1990) taken with the Advanced CCD Imaging Spectrometer (ACIS) on board the Chandra X-ray Observatory. The observation exposure is 40 ks and was completed on 2001 December 18. The entire SNR is imaged on one back-illuminated chip (called “S3,” see Figure 2). Data reprocessing was conducted using the standard processing procedures in the CIAO v4.12 (Fruscione et al. 2006) software package. The cleaned spectra are then extracted and background subtracted using one large annulus-shaped region surrounding the remnant. We model both SNR and PWN emission components using data extracted from the regions indicated in Figure 2 and perform a spectral analysis. A spectrum for each component is extracted using the

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top: 0.5–7 keV X-ray data and best-fit models for the two source models. The green solid line represents the nonthermal component from the PWN and the solid orange and blue lines represent the first and second thermal components of the SNR spectrum, respectively. Bottom: the residuals of the difference in the best-fit model and data for the SNR spectral fit (black, top) and the PWN spectral fit (red, bottom) in units of sigma.</figcaption></figure>

specextract tool in CIAO and modeled using Sherpa within CIAO (Freeman et al. 2001).

### 3.2.2. Chandra X-Ray Data Analysis Results

Data between 0.5 and 7 keV are used to model observed emission and is binned to at least 20 counts per bin. We fit the two source regions for the SNR and PWN components simultaneously and the best-fit model is displayed in Figure 3. A two-component collisionally ionized plasma model (xsvapec) is found to best describe the emission from the SNR and one nonthermal powlaw1d model is preferred for the PWN component (similar to prior works, e.g., Haberl et al. 2012; McEntaffer et al. 2012). We account for interstellar absorption along the line of sight by including the tbabs hydrogen column density parameter, which uses the abundances estimated by Wilms et al. (2000).

The initial values of elemental abundances are set to those estimated for the LMC in Russell & Dopita (1992) and are allowed to vary one by 1 in each fit iteration. We keep the abundance of an element free if it significantly improves the fit; otherwise, the value remains frozen at the following abundances relative to solar: He 0.89, C 0.26, N 0.16, O 0.32, Ne 0.42, Mg 0.74, Si 1.7, S 0.27, Ar 0.49, Ca 0.33, Fe 0.50, and Ni 0.62. Aluminum is not well constrained (see Section 4.3 in Russell & Dopita 1992) so we freeze its value to 0.33.

The PWN spectrum is contaminated by two thermal components from the SNR emission in addition to a nonthermal component described best as a power law. Because SNR emission contaminates the PWN emission, we link the thermal parameters of the two models using the scale1d
<table><caption>**Table 1**
Summary of the Statistics and Best-fit Model for the SNR and PWN Components in the X-Ray Analysis</caption><thead><tr><th>Data Points</th><th>Dof<sup>a</sup></th><th>Reduced χ<sup>2</sup></th></tr></thead><tbody><tr><td>204</td><td>191</td><td>0.94</td></tr><tr><th>Component</th><th colspan="2">Model</th></tr><tr><td>SNR</td><td colspan="2">tbabs × (vapec<sub>1</sub>+vapec<sub>2</sub>)</td></tr><tr><td>PWN</td><td colspan="2">tbabs × [(c<sub>1</sub> × vapec<sub>1</sub>+c<sub>2</sub> × vapec<sub>2</sub>) + powlaw]</td></tr></tbody></table>

**Notes.** The thermal components of the PWN spectrum are linked to the SNR model with the free coefficients $c_1$ and $c_2$.

<sup>a</sup> Degrees of freedom

parameter in Sherpa (Table 1). We leave the amplitude, $C_0$, free to vary in the fit for both thermal components.

The hydrogen column density is $N_H = 3.7^{+1.1}_{-0.9} \times 10^{21} \text{ cm}^{-2}$, the PWN power-law index is $\Gamma_X = 1.74^{+0.20}_{-0.20}$, and the unabsorbed X-ray flux of the PWN component between 0.5 and 7 keV is $f_X = 2.68 \pm 0.59 \times 10^{-13} \text{ erg cm}^{-2} \text{ s}^{-1}$. The $N_H$ value is reasonable compared to that measured in the direction of the LMC,<sup>8</sup> $N_H = 2.2 \times 10^{21} \text{ cm}^{-2}$ (Blackburn 1995). The best-fit model is consistent with other X-ray analyses (Haberl et al. 2012; McEntaffer et al. 2012), with the largest differences being the elemental abundances, which can be explained by the use of the Wilms et al. (2000) and the Verner et al. (1996) cross sections in this work, in addition to slight differences in the choice of model components for the thermal emission and detector capabilities. In particular, Haberl et al. (2012) analyzed XMM-Newton observations of the entire SNR, but the PWN is not resolved and thus only one global spectrum was used to characterize any SNR and PWN emission. The SNR is much brighter than the PWN in X-rays so the nonthermal component from the PWN in the XMM-Newton X-ray spectrum is not well constrained. McEntaffer et al. (2012) used Anders & Grevesse (1989) abundances and Balucinska-Church & McCammon (1992) cross sections, and instead of two thermal equilibrium models, vapec, their best-fit model assumes a two-component structure from a vapec+vnei combination, where the vnei models the second thermal component without ionization equilibrium conditions.

The best-fit temperatures for the two-component thermal model used to describe SNR emission are $kT = 0.34^{+0.02}_{-0.05}$ keV and $kT = 0.16^{+0.01}_{-0.01}$ keV, similar to that reported in McEntaffer et al. (2012). The PWN spectrum is nonthermal and best fit with a power-law and photon index, $\Gamma_X = 1.74^{+0.20}_{-0.20}$. The PWN's spectral index is slightly harder than that reported in McEntaffer et al. (2012), where an index $\Gamma_X \sim 2$ across the PWN region is measured, but is still in agreement within the 90% C.L. uncertainties. No synchrotron component is attributed to the SNR, but we estimate the 0.5–7 keV 90% C. L. upper limit of the flux for a nonthermal component to the SNR spectrum to be $F_X < 5.5 \times 10^{-13} \text{ erg cm}^{-2} \text{ s}^{-1}$. The best-fit parameters are listed in Table 2 along with the corresponding 90% confidence intervals using the conf tool in Sherpa.

<table><caption>**Table 2**
Summary of the 90% C.L. Statistics and Parameters for the Best-fit Model for Each Component in the X-Ray Analysis</caption><thead><tr><th>SNR<br>Component</th><th>Parameter</th><th>Best-fit Value</th></tr></thead><tbody><tr><td>tbabs<sup>a</sup></td><td>N<sub>H</sub>(10<sup>22</sup> cm<sup>-2</sup>)</td><td>0.37<sup>+0.11</sup><sub>-0.09</sub></td></tr><tr><td>vapec<sub>1</sub></td><td>kT(keV)</td><td>0.34<sup>+0.02</sup><sub>-0.05</sub></td></tr><tr><td>vapec<sub>2</sub></td><td>Normalization</td><td>3.67<sup>+2.55</sup><sub>-0.97</sub> × 10<sup>-3</sup></td></tr><tr><td></td><td>kT (keV)</td><td>0.16<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td></td><td>O</td><td>0.35<sup>+0.26</sup><sub>-0.11</sub></td></tr><tr><td></td><td>Ne</td><td>0.39<sup>+0.32</sup><sub>-0.13</sub></td></tr><tr><td></td><td>Mg</td><td>0.56<sup>+0.50</sup><sub>-0.33</sub></td></tr><tr><td></td><td>Fe</td><td>&lt; 0.70</td></tr><tr><td></td><td>Normalization</td><td>0.05<sup>+0.06</sup><sub>-0.03</sub></td></tr><tr><th>PWN<br>Component</th><th>Parameter</th><th>Best-fit Value</th></tr><tr><td>c<sub>1</sub></td><td>C<sub>0</sub></td><td>0.07<sup>+0.02</sup><sub>-0.05</sub></td></tr><tr><td>c<sub>2</sub></td><td>C<sub>0</sub></td><td>0.14<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td>powlaw</td><td>Γ</td><td>1.74<sup>+0.20</sup><sub>-0.20</sub></td></tr><tr><td></td><td>Amplitude</td><td>5.28<sup>+1.18</sup><sub>-1.01</sub> × 10<sup>-5</sup></td></tr></tbody></table>

**Notes.** Metal abundances are reported in solar units.

<sup>a</sup> Absorption cross section set to Verner et al. (1996)

### 3.3. Gamma Ray

#### 3.3.1. Fermi-LAT

The Fermi Gamma-ray Space Telescope houses the Large Area Telescope (LAT; Atwood et al. 2009). The LAT instrument is sensitive to gamma rays with energies from 50 MeV to >300 GeV (Abdollahi et al. 2020) and has been continuously surveying the entire sky every 3 hr since beginning operation in 2008 August.

We use 11.5 yr (from 2008 August to 2020 January) of Pass 8 SOURCE class data (Atwood et al. 2013; Bruel et al. 2018) between 300 MeV and 2 TeV. Photons detected at zenith angles larger than 100° were excluded to limit the contamination from gamma rays generated by CR interactions in the upper layers of Earth's atmosphere.

#### 3.3.2. Fermi-LAT Data Analysis

We perform a binned likelihood analysis with the latest FermiTools package<sup>9</sup> (v2.0.8) and FermiPy Python 3 package (v1.0.1, Wood et al. 2017), utilizing the P8R3_SOURCE_V3 instrument response function (IRF) and account for energy dispersion, to perform data reduction and analysis. We organize the events by point spread function (PSF) type using evtype=4, 8, 16, 32 to represent the 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.<sup>10</sup> We fit the square 10° ROI centered on the PWN position in equatorial coordinates using a pixel bin size 0.05 and 10 bins per decade in energy (38 total bins). The gamma-ray sky for the ROI is modeled from the latest

<footnote n="8"><sup>8</sup> Using the nh tool from the HEASARC FTOOLS package http://heasarc.gsfc.nasa.gov/ftools.</footnote>

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

<footnote n="10"><sup>10</sup> See FermiPy documentation for details: https://fermipy.readthedocs.io/en/0.6.8/config.html.</footnote>
comprehensive Fermi-LAT source catalog based on 10 yr of data, 4FGL (data release 2 (DR2), Ballet et al. 2020) for point and extended sources<fnref n="11" /> that are within 15° of the ROI center, as well as the latest Galactic diffuse and isotropic diffuse templates (gll_iem_v07.fits and iso_P8R3_SOUR-CE_V3_v1.txt, respectively).<fnref n="12" />

Because B0453-685 is located in the LMC, we need to properly account for the diffuse emission from the LMC. We employ in the 4FGL source model four additional extended source components to reconstruct the emissivity model developed in Ackermann et al. (2016) to represent the diffuse LMC emission. The four additional sources are 4FGL J0500.9-6945e (LMC far west), 4FGL J0519.9-6845e (LMC galaxy), 4FGL J0530.0-6900e (30 Dor west), and 4FGL J0531.8-6639e (LMC north). These four extended templates along with the isotropic and Galactic diffuse templates define the total background for the ROI.

With the source model described above, we allow the background components and sources with test statistic (TS) ≥25 and distances from the ROI center ≤3°0 to vary in spectral index and normalization. We computed a series of diagnostic TS and count maps in order to search for and understand any residual gamma-ray emission. The TS value is defined to be the natural logarithm of the ratio of 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)$$

The TS value quantifies the significance of a source detection with a given set of location and spectral parameters and the significance of such detection can be estimated by taking the square root of the TS value for 1 dof (Mattox et al. 1996). TS values >25 correspond to a detection significance >4σ for 4 dof.

We generated the count and TS maps for 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="13" /> We inspected the TS maps for additional sources, finding a faint pointlike gamma-ray source coincident in location with B0453-685 and no known 4FGL counterpart.<fnref n="14" />

Figure 1, left panel, demonstrates the total source model used in the analysis (except the isotropic and Galactic diffuse templates). Three additional point sources are added to the source model that models residual emission in the field of view (PS1, PS2, and PS3 in the right panel of Figure 4). PS3 corresponds to 4FGL-DR3 source J0517.9-6930c. A count and TS map between energies 1 and 10 GeV are shown in Figure 4 where the TS map, right panel, is generated from the global source model, which has no associated source at the position of B0453-685. Faint gamma-ray emission is visible and coincident with the SNR B0453-685.

### 3.3.3. Results of Fermi-LAT Data Analysis

To model the gamma-ray emission coincident with B0453-685 we add a point source at the PWN location (R.A., decl.) J2000 = (73°,408, -68°,489) to the 300 MeV–2 TeV source model. With a fixed location, we set the spectrum to a power law characterized by a photon index $\Gamma = 2$,

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

$E_0$ is set to 1000 MeV. We then allow the spectral index and normalization to vary. The TS value for a point source with a power-law spectrum and photon index, $\Gamma = 2.3 \pm 0.2$, is 23. We investigate the spectral properties of the gamma-ray emission by changing the spectral model to a log parabola shape following the definition,<fnref n="15" />

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

We fix $E_b = 4.0$ GeV but allow $\alpha, \beta$, and $N_0$ to vary in the fit. The TS value of a point source at the PWN/SNR position with a log parabola spectrum is 27 and has $\alpha = 2.5 \pm 0.4$ and $\beta = 0.5 \pm 0.3$. We test the spectral parameters once more using a spectrum typically observed with megaelectronvolt–gigaelectronvolt pulsars, a power law with a super-exponential cutoff (PLEC):<fnref n="16" />

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

where $E_0$ is the scale (set to 1000 MeV). The TS value of a point source at the position of B0453-685 with a PLEC spectrum is 27 and has $\Gamma_1 = 0.8 \pm 0.8$, $\Gamma_2$ is fixed to 0.7, and exponential factor $a = 0.009 \pm 0.005$, which corresponds to an $E_c \sim 1$ GeV energy cutoff. See Table 3 for a summary of the spectral parameters for each point-source test.

Fermi-LAT pulsars are often characterized as either a power-law or a PLEC spectrum and typically cut off at energies <10 GeV (e.g., Abdo et al. 2013). While we cannot firmly rule out that the observed gamma-ray emission is from the still-undetected pulsar based on the best-fit spectral parameters, it seems unlikely given the majority of the emission is measured in 1–10 GeV. Between the three tested spectral models, the log parabola and PLEC are only marginally preferred (e.g., TS\_LogParabola = 2Δ(ln $\mathcal{L}$) = 5.9) and carry another degree of freedom with respect to the power-law spectral model. We therefore conclude that the best characterization for the gamma-ray emission coincident with SNR B0453-685 is a power-law spectrum. The corresponding gamma-ray SED is displayed in Figure 5.

We localize the point source modeled using a power-law spectrum with GTAnalysis.localize to find the best-fit position and uncertainty. The localized position for the new gamma-ray source is offset by 0°.01 from the exact position of B0453-685 and has R.A., decl. = 73°,39, -68°,49 (J2000). The corresponding 95% positional uncertainty radius is $r = 0°.12$.

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

<footnote n="12">¹² LAT background models and appropriate instrument response functions: https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html.</footnote>

<footnote n="13">¹³ 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="14">¹⁴ The closest 4FGL source is the probable unclassified blazar 4FGL J0511.4-6804 ~2° away.</footnote>

<footnote n="15">¹⁵ For a review of Fermi-LAT source spectral models see https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html.</footnote>

<footnote n="16">¹⁶ This follows the PLSuperExpCutoff2 form used for the 4FGL-DR2. Details can be found at https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html#PLSuperExpCutoff2.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Left: smoothed ($\sigma = 0^{\circ}1$) 5 $^{\circ} \times 5^{\circ}$ count map of PSF3 events between 1 and 10 GeV with the locations of 4FGL sources in the field of view labeled. The pixel size is 0$^{\circ}01$ pixel$^{-1}$. Right: 5$^{\circ} \times 5^{\circ}$ TS map between 1 and 10 GeV. The maximum TS value at the SNR position is ~28. The 95% positional uncertainty for the best-fit gamma-ray point source is in blue. In both panels, the location and approximate size of the composite SNR B0453-685 ($r = 0^{\circ}02$) is marked in white with radius $r = 0^{\circ}05$.</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 templates. Both of these extended templates assume a symmetric 2D shape with width parameters radius and sigma, respectively. We fix the position but keep spectral parameters free to vary when finding the best-fit spatial extension for both templates. The best-fit parameters for the extension tests are presented in Table 4. The faint gamma-ray source does not display significant extension, consistent with the size of B0453-685 if observed by Fermi. We also perform a variability analysis following the method in the 4FGL catalogs using 1 yr time bins. There is no significant variability observed ($TS_{var} < 2$). Finally, we search the new gamma-ray source's 95% uncertainty region for the spatial overlap with any other objects that may be able to explain the observed gamma-ray emission. There are more than 150 LMC stars within the confidence region, but SNR B0453-685 is the only nonstellar object.<fnref n="17" />

### 3.3.4. Systematic Error from the Choice of Interstellar Emission Model and IRF

We account for systematic uncertainties introduced by the choice of the interstellar emission model (IEM) and the IRFs, which mainly affect the spectrum of the measured gamma-ray emission. We have followed the prescription developed in de Palma et al. (2013) and Acero et al. (2016), based on generating 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 Galactic Extended Source Catalog (Ackermann et al. 2017). The gamma-ray point source coincident with SNR B0453-685 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 from the effective area<fnref n="18" /> 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 negligible for B0453-685, which is not surprising given the location of the LMC with respect to the bright diffuse gamma-ray emission along the Galactic plane.

### 3.3.5. Systematic Error from Choice of Diffuse LMC

We must also account for the systematic error that is introduced by having an additional diffuse background component. This third component is attributed to the CR population of the LMC interacting with the LMC ISM and there are limitations to the accuracy of the background templates used to model this emission, similar to the Galactic diffuse background. We can probe these limitations by employing the straightforward method described in Ackermann et al. (2016) to measure systematics from the diffuse LMC. This requires replacing the four extended sources that represent the diffuse LMC in this analysis (the emissivity model, Ackermann et al. 2016) with four different extended sources to represent an alternative template for the diffuse LMC (the analytic model, Ackermann et al. 2016). The gamma-ray point source coincident with SNR B0453-685 is then refit with the alternative diffuse LMC template to obtain a new spectral flux that we then compare with the results of the emissivity model following Equation (5) in Acero et al. (2016). The systematic error from the choice of the diffuse LMC template is largest in the two lowest-energy bins, but negligible

<footnote n="17">17 https://simbad.u-strasbg.fr/simbad/sim-fcoo</footnote>

<footnote n="18">¹⁸ https://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html</footnote>
<table><caption>**Table 3**
Summary of the Best-fit Parameters and the Associated 68% C.L. Statistics for all Point-source Models Tested</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>G<sub>E</sub> (MeV cm<sup>-2</sup> s<sup>-1</sup>)</th><th>E<sub>b</sub> or a</th><th>TS</th></tr></thead><tbody><tr><td>Power law</td><td>-505,673</td><td>2.3 ± 0.2</td><td>...</td><td>...</td><td>7.5(± 2.2) × 10<sup>-7</sup></td><td>...</td><td>23</td></tr><tr><td>Log parabola</td><td>-505,670</td><td>...</td><td>2.5 ± 0.4</td><td>0.5 ± 0.3</td><td>5.0(± 1.4) × 10<sup>-7</sup></td><td>4000</td><td>27</td></tr><tr><td>Power law with exponential cutoff</td><td>-505,673</td><td>...</td><td>0.8 ± 0.8</td><td>0.7 (fixed)</td><td>5.1(± 1.3) × 10<sup>-7</sup></td><td>0.009 ± 0.005</td><td>27</td></tr></tbody></table>

Note. $G_E$ is the integrated energy flux for energies 300 MeV–2 TeV. The units for $E_b$ are MeV. The units for the exponential factor $a$ are MeV$^{-\Gamma_2}$.

<table><caption>**Table 4**
Summary of the Best-fit Parameters and the Associated Statistics for Each Spatial Template Used in Our Analysis</caption><thead><tr><th>Spatial Template</th><th>TS</th><th>TS<sub>ext</sub></th><th>95% Radius Upper Limit (°)</th></tr></thead><tbody><tr><td>Point source</td><td>23</td><td>...</td><td>...</td></tr><tr><td>Radial disk</td><td>23</td><td>0.1</td><td>0.2</td></tr><tr><td>Radial Gaussian</td><td>23</td><td>0.1</td><td>0.2</td></tr></tbody></table>

in higher-energy bins. The corresponding systematic error is plotted in Figure 5 in black.

# 4. Broadband Modeling

## 4.1. Investigating Gamma-Ray Origin

For a gamma-ray source at $d = 50$ kpc, the 300 MeV–2 TeV gamma-ray luminosity is $L_γ = 2.6 \times 10^{35}$ erg s$^{-1}$. We compare this value and the best-fit spectral index $Γ_γ = 2.3$ to Figure 17 in Acero et al. (2016), which plots the gigaelectronvolt luminosity against the power-law index for Fermi-LAT detected SNRs. There is a correlation between the gigaelectronvolt properties and age of an SNR, in particular, the softest (i.e., oldest) SNRs have larger gigaelectronvolt luminosities than harder (i.e., younger) SNRs. Comparing the gigaelectronvolt luminosity found here to those shown in Figure 17 in Acero et al. (2016), we see that the gamma-ray source is in agreement with the evolved SNRs. This observed correlation is likely due to evolved SNRs interacting with dense material (Acero et al. 2016), yet the SNR shell associated with B0453-685 does not show compelling evidence for such an interaction (Figure 1). We also compare the gigaelectronvolt luminosity to those of Fermi-LAT detected pulsars and PWNe (Abdo et al. 2013; Acero et al. 2013), finding that the gigaelectronvolt luminosity is characteristic of both source classes. Moreover, the spectral index $Γ_γ = 2.3 \pm 0.2$ is in agreement with Fermi-LAT detected SNRs, PWNe, and pulsars: $Γ_{γ,SNRs} \approx 2.3$, $Γ_{γ,PWNe} \approx 2.1$, and $Γ_{γ,PSRs} \approx 2.3$ are the average power-law indices for SNRs, PWNe, and pulsars in the 4FGL-DR2 catalog, respectively (Ballet et al. 2020).

In order to investigate the origin of the observed gamma-ray emission, we use the NAIMA Python package (v0.10.0 Zabalza 2015), which computes the radiation from a single nonthermal relativistic particle population and performs a Markov Chain Monte Carlo (MCMC) sampling of the likelihood distributions (using the `emcee` package, Foreman-Mackey et al. 2013). For the particle distribution in energy, we assume a power-law shape with an exponential cutoff,

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

We then test a combination of free parameters (namely, the normalization $A$, index $Γ$, energy cutoff $E_c$, and magnetic field

<figure><img src="image_7.png" /><figcaption>**Figure 5.** The best-fit gamma-ray SED for B0453-685 with $1\sigma$ statistical uncertainties in red for TS > 1 and 95% C.L. upper limits otherwise. The systematic error from the choice of diffuse LMC model is plotted in black. TS values for each spectral bin are plotted as the green histogram. The data are best characterized as a power law with $Γ = 2.3 \pm 0.2$.</figcaption></figure>

B) that can best explain the broadband spectra for the SNR and PWN independently. We consider only one photon field in all IC scattering calculations in this section, the cosmic microwave background (CMB).

### 4.1.1. PWN as Gamma-Ray Origin

The radio spectrum considered together with the hard X-ray spectrum of the PWN strongly indicate the presence of more than one particle population, which is also indicated by the estimated age and evolutionary phase of the host SNR. Moreover, the observed X-ray morphology of the nebula displays features consistent with an evolved SNR where the reverse shock has impacted the PWN, compressing the population of previously injected particles while the central pulsar continues to inject new high-energy particles (e.g., Gaensler et al. 2003; Haberl et al. 2012). The return of the reverse shock would additionally explain the significantly enhanced abundances relative to the local ISM, indicating the PWN plasma is becoming ejecta dominated (McEntaffer et al. 2012). Based on this, we instead incorporate two-leptonic particle populations under the same conditions (nebular magnetic field and ambient photon fields) and combine them to represent a two-leptonic broadband model. A two-leptonic broadband model can describe well the PWN radio, X-ray, and gamma-ray data, where the lower-energy particles dominate the radio and gamma-ray emission while the higher-energy particles are losing more energy in synchrotron radiation than in IC radiation, and therefore dominate in X-ray. We allow Population 1, the lower-energy population, to constrain the magnetic field strength, as the oldest particles likely dominate the synchrotron emission (Gelfand et al. 2009). It is possible each population is interacting with magnetic field regions of varying strength, but for simplicity, we fix the magnetic field value to the best fit found from the lower-energy population's broadband
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The best-fit broadband models for the three scenarios investigated to understand gamma-ray origin. Top: two-leptonic populations are required to explain the broadband PWN emission. Bottom left: a single leptonic population describing SNR synchrotron emission combined with a single hadronic population describing the gamma-ray emission via pion decay from the SNR. Bottom Right: the case where the leptonic population dominates over the hadronic population in the SNR. Radio data of PWN (blue) and SNR (purple) are from Haberl et al. (2012), X-ray data of PWN (cyan) and SNR (yellow) are described in detail in Section 3.2.2, and gamma-ray data (gray/green) in Section 3.3.3. The uncertainties to the corresponding radio data are very small at this flux scale.</figcaption></figure>

<table><caption>**Table 5**
Summary of the Statistics and Best-fit Models for the PWN and SNR Broadband Models Displayed in Figure 6</caption><thead><tr><th></th><th colspan="2">Two-leptonic PWN</th><th>Hadronic-dominant SNR</th><th>Leptonic-dominant SNR</th></tr><tr><th></th><th>Population 1</th><th>Population 2</th><th>(Hadrons Only)</th><th>(Leptons Only)</th></tr></thead><tbody><tr><td>Maximum log-likelihood (dof)</td><td>-2.07 (13-4)</td><td>-8.71 (13-3)</td><td>-2.06 (13-3)</td><td>-1.67 (13-4)</td></tr><tr><td colspan="5">Maximum likelihood values</td></tr><tr><td>W<sub>e</sub> or W<sub>p</sub><sup>a</sup></td><td colspan="2">2.84 &times; 10<sup>49</sup> 1.43 &times; 10<sup>47</sup></td><td>3.89 &times; 10<sup>51</sup></td><td>2.71 &times; 10<sup>50</sup></td></tr><tr><td>Index</td><td colspan="2">0.88 &plusmn; 0.13 2.05 &plusmn; 0.62</td><td>1.95 (fixed)</td><td>1.95 &plusmn; 0.05</td></tr><tr><td>log<sub>10</sub> E<sub>c</sub><sup>b</sup></td><td colspan="2">-0.45 &mu;m 0.11 2.35 &plusmn; 0.26</td><td>-0.68 &plusmn; 0.41</td><td>-0.17 &plusmn; 0.15</td></tr><tr><td>B<sup>c</sup></td><td>8.18 &plusmn; 4.25</td><td>8.18 (fixed)</td><td>4.82 &plusmn; 0.12</td><td>1.47 &plusmn; 0.29</td></tr></tbody></table>

**Notes.** The maximum log-likelihood can be understood as $\chi^2 = -2\ln L$.

<sup>a</sup> The total particle energy $W_e$ or $W_p$ is in units of erg.

<sup>b</sup> Logarithm base 10 of the cutoff energy in units of teraelectronvolts.

<sup>c</sup> Magnetic field in units of $\mu$Gauss.

model when searching for a model fit for the higher-energy population, $B \sim 8\mu$ G. The best-fit parameters for the low-energy population are $\Gamma = 0.88 \pm 0.13$ and $E_c = 0.35 \pm 0.11$ TeV. The best-fit parameters for the high-energy population are $\Gamma = 2.05 \pm 0.62$ and $E_c = 224^{+183}_{-101}$ TeV. The best-fit two-leptonic broadband model for the PWN is displayed in the top panel of Figure 6 and the corresponding

best-fit parameters for both particle populations are listed in Table 5.

The two-leptonic broadband model for the PWN has an estimated total particle energy of $W_e = 2.86 \times 10^{49}$ erg. The lower-energy population is responsible for $W_e = 2.84 \times 10^{49}$ erg and the higher-energy population with the remainder, $W_e = 1.43 \times 10^{47}$ erg.
### 4.1.2. SNR as Gamma-Ray Origin

There are two possible scenarios for the SNR to be responsible for the gamma-ray emission. The first is a single leptonic population that is accelerated at the SNR shock front, generating both synchrotron emission at lower energies and IC emission at higher energies in gamma rays (e.g., Reynolds 2008). The second scenario is a single leptonic population emitting mostly synchrotron radiation at lower energies combined with a single hadronic population emitting gamma rays through pion decay. We describe both of these models and their implications below.

To model the lower-energy SNR emission together with the newly discovered Fermi-LAT emission using a single leptonic population (i.e., the leptonic-dominant scenario), we require a particle index $\Gamma = 1.95 \pm 0.05$, an energy cutoff at 671 GeV, and an inferred magnetic field $B = 1.47$ µG. For the hadronic-dominant scenario, we model the broadband SNR emission assuming a single leptonic and single hadronic population. We measure the magnetic field value to be $B = 4.82 \pm 0.12$ µG for the synchrotron component under the electron-to-proton ratio assumption $k_{ep} = 0.01$ (Castro et al. 2011) and characterizing the gamma-ray emission via pion decay through proton-proton collisions at the SNR shock front. The pre-shock proton density $n_0$ has been estimated to be $\sim 0.4$ cm$^{-3}$ from the SNR X-ray emission measured along the rim region (Gaensler et al. 2003). The post-shock proton density at the SNR forward shock $n_H$ could be about 4 times as high as $n_0$; thus, for a compression ratio $n_H/n_0 = 4$, $n_H \sim 1.6$ cm$^{-3}$ (e.g., Vink & Laming 2003). We fix the target proton density $n_H = 1.6$ cm$^{-3}$ at the default differential cross section (Pythia8, Zabalza 2015) while also fixing the proton particle index to $\Gamma = 1.95$. The latter choice is motivated by the particle index being well defined from the radio data in the leptonic population, but is not well constrained for the hadronic population. We measure an energy cutoff $E_c = 0.194^{+0.27}_{-0.11}$ TeV for the proton spectrum that can best reproduce the observed gamma-ray spectrum. The best-fit broadband models for the SNR are displayed in the lower panels of Figure 6 and the corresponding parameters are listed in Table 5.

The best-fit leptonic-dominant model for the SNR yields a total electron energy of $W_e = 2.71 \times 10^{50}$ erg. This implies, assuming $k_{ep} = 0.01$, the total proton energy from undetected pion decay emission is $W_p = W_e \times 100 = 2.71 \times 10^{52}$ erg, requiring roughly 20 times the canonical expectation $E = 10^{51}$ erg to be in total SNR CR energy alone and a very low target density $n_H = 0.05$ cm$^{-3}$. The best-fit hadronic-dominant model requires a total proton energy of $W_p = 3.89 \times 10^{51}$ erg, a factor of almost 4 times greater than the typical supernova explosion energy of $E = 10^{51}$ erg.

Furthermore, the inferred magnetic field $B = 1.47$ µG in the leptonic-dominant model is comparable to the coherent component of the LMC magnetic field $B \sim 1$ µG (Gaensler et al. 2005), which is weaker than one would expect at the SNR shock front, where shock compression can amplify the magnetic field strength 4–5 times the initial value (see, e.g., Castro et al. 2011, and references therein). In order to explain the observed gamma-ray emission via pion decay with a reasonable energy in accelerated protons ($E \sim 10^{50}$ erg), the SNR must be interacting with dense ambient material (e.g., similar to W44 and IC 443, Ackermann et al. 2013; Chen et al. 2014; Slane et al. 2015). The radio and X-ray observations of the SNR show a fainter, limb-brightened shell compared to the

bright, compact central PWN, providing little evidence of the SNR forward shock colliding with ambient media.

In conclusion, the energetics inferred from the SNR models lead us to favor the two-leptonic PWN broadband model as the most likely explanation for the gamma-ray emission reported here. We explore the most accurate representation of the PWN broadband data while also exploring the likelihood of a pulsar contribution in the following section.

### 4.2. PWN Evolution through Semi-analytic Modeling

We have established in the previous section that modeling the nonthermal broadband SED suggests that it most likely originates from two populations of leptons with different energy spectra, similar to that expected for evolved PWNe once they have collided with the SNR reverse shock (see, e.g., Gelfand et al. 2009; Temim et al. 2015). To determine if this depicted scenario can explain the intrinsic properties of this system, we model the observed properties of the PWN, assuming it is responsible for the detected Fermi-LAT gamma-ray emission, as it evolves inside the composite SNR B0453-685.

We use the dynamical and radiative properties of a PWN predicted by an evolutionary model, similar to that described by Gelfand et al. (2009), to identify the combination of neutron star, pulsar wind, supernova explosion, and ISM properties that can best reproduce what is observed. The model is developed using an MCMC fitting procedure (see, e.g., Gelfand et al. 2015, for details) to find the combination of free parameters that can best represent the observations. The observed sizes of the SNR and PWN together with the radio, X-ray, and gamma-ray data are used to calculate the final broadband model at an age, $t_{age}$. The predicted dynamical and radiative properties of the PWN that correspond to the best representation of the broadband data are listed in Table 6. The parameters velpsr, etag, kpsr, gpsr, and ecut are fixed to zero.

The analysis performed here is similar to that previously reported for MSH 15-56 (Temim et al. 2013), G54.1+0.3 (Gelfand et al. 2015) G21.5-0.9 (Hattori et al. 2020), Kes 75 (Gotthelf et al. 2021), and HESS J1640-465 (Mares et al. 2021). For the characteristic age $t_{ch}$ of a pulsar (see Pacini & Salvati 1973; Gaensler & Slane 2006), the age $t_{age}$ is defined as

$$t_{age} = \frac{2t_{ch}}{p-1} - \tau_{sd}, \qquad (6)$$

and the spin-down luminosity $\dot{E}$ is defined as

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

and are chosen for a braking index $p$, initial spin-down luminosity $\dot{E}_0$, and spin-down timescale $\tau_{sd}$ to best reproduce the pulsar’s likely characteristic age and current spin-down luminosity. A fraction $\eta_\gamma$ of this luminosity is converted to gamma-ray emission from the neutron star’s magnetosphere, and the rest ($1 - \eta_\gamma$) is injected into the PWN in the form of a magnetized, highly relativistic outflow, i.e., the pulsar wind. The pulsar wind enters the PWN at the termination shock, where the rate of magnetic energy $\dot{E}_B$ and particle energy $\dot{E}_P$ injected into the PWN is expressed as

$$\dot{E}_B(t) \equiv (1 - \eta_\gamma) \eta_B \dot{E}(t), \qquad (8)$$
<table><caption>**Table 6**
Summary of the Input Parameters for the Evolutionary System and their Best-fit Values Considering PWN-only and PWN+PSR Contributions to the Fermi-LAT Emission</caption><thead><tr><th>Shorthand</th><th>Parameter</th><th>PWN Best Fit</th><th>PWN+PSR Best Fit</th><th>Units</th></tr></thead><tbody><tr><td>loglh</td><td>Log-likelihood of SED</td><td>-19.9</td><td>-17.6</td><td>...</td></tr><tr><td>esn</td><td>Initial kinetic energy of supernova ejecta</td><td>5.24</td><td>5.21</td><td>10<sup>50</sup> erg</td></tr><tr><td>mej</td><td>Mass of supernova ejecta</td><td>2.24</td><td>2.42</td><td>M<sub>☉</sub></td></tr><tr><td>nism</td><td>Number density of surrounding ISM</td><td>0.97</td><td>1.00</td><td>cm<sup>-3</sup></td></tr><tr><td>brakind</td><td>Pulsar braking index</td><td>2.89</td><td>2.83</td><td>...</td></tr><tr><td>tau</td><td>Pulsar spin-down timescale</td><td>172</td><td>166</td><td>yr</td></tr><tr><td>age</td><td>Age of system</td><td>13,900</td><td>14,300</td><td>yr</td></tr><tr><td>e0</td><td>Initial spin-down luminosity of pulsar</td><td>6.95</td><td>6.79</td><td>10<sup>39</sup> erg s<sup>-1</sup></td></tr><tr><td>etag</td><td>Fraction of spin-down luminosity lost as radiation</td><td>≡0</td><td>0.246</td><td>...</td></tr><tr><td>etab</td><td>Magnetization of the pulsar wind</td><td>0.0006</td><td>0.0007</td><td>...</td></tr><tr><td>emin</td><td>Minimum particle energy in pulsar wind</td><td>1.77</td><td>2.26</td><td>GeV</td></tr><tr><td>emax</td><td>Maximum particle energy in pulsar wind</td><td>0.90</td><td>0.73</td><td>PeV</td></tr><tr><td>ebreak</td><td>Break energy in pulsar wind</td><td>76</td><td>72</td><td>GeV</td></tr><tr><td>p1</td><td>Injection index below the break (dN/dE ~ E<sup>-p1</sup>)</td><td>1.47</td><td>1.34</td><td>...</td></tr><tr><td>p2</td><td>Injection index below the break (dN/dE ~ E<sup>-p2</sup>)</td><td>2.36</td><td>2.36</td><td>...</td></tr><tr><td>ictemp</td><td>Temperature of each background photon field</td><td>1.02</td><td>1.13</td><td>10<sup>6</sup> K</td></tr><tr><td>icnorm</td><td>Log normalization of each background photon field</td><td>-17.9</td><td>-18.0</td><td>...</td></tr><tr><td>gpsr</td><td>Photon index of the gamma rays produced directly by the pulsar</td><td>...</td><td>2.00</td><td>...</td></tr><tr><td>ecut</td><td>Cutoff energy from the power law of the pulsar contribution</td><td>...</td><td>3.21</td><td>GeV</td></tr></tbody></table>

$$ \dot{E}_{\mathrm{P}}(t) \equiv (1 - \eta_{\gamma}) \eta_{\mathrm{P}} \dot{E}(t), \qquad (9) $$

where $\eta_B$ is the magnetization of the wind and defined to be the fraction of the pulsar's spin-down luminosity injected into the PWN as magnetic fields and $\eta_P$ is the fraction of spin-down luminosity injected into the PWN as particles. We assume the PWN IC emission results from leptons scattering off of the CMB similar to the previous modeling section; however, the total particle energy and the properties of the background photon fields cannot be independently determined. Since the evolutionary model accounts for the decline in total particle energy from the adiabatic losses of early PWN evolution and the increase of synchrotron losses at later times from compression, where both likely have a significant effect on the oldest particles, a second photon field is hence required. The second, ambient photon field is defined by temperature $T_{\mathrm{IC}}$ and normalization $K_{\mathrm{IC}}$, such that the energy density of the photon field $u_{\mathrm{IC}}$ is

$$ u_{\mathrm{IC}} = K_{\mathrm{IC}} a_{\mathrm{BB}} T_{\mathrm{IC}}^{4}, \qquad (10) $$

where $a_{\mathrm{BB}} = 7.5657 \times 10^{-15}$ erg cm$^{-3}$ K$^{-4}$. Additionally, we assume the particle injection spectrum at the termination shock

is well described by a broken power-law distribution

$$ \frac{d\dot{N}_{e^{\pm}}(E)}{dE} = \begin{cases} \dot{N}_{\mathrm{break}} \left(\frac{E}{E_{\mathrm{break}}}\right)^{-p_1} & E_{\mathrm{min}} < E < E_{\mathrm{break}} \\ \dot{N}_{\mathrm{break}} \left(\frac{E}{E_{\mathrm{break}}}\right)^{-p_2} & E_{\mathrm{break}} < E < E_{\mathrm{max}} \end{cases}, \qquad (11) $$

where $\dot{N}_{e^{\pm}}$ is the rate that electrons and positrons are injected into the PWN, and $\dot{N}_{\mathrm{break}}$ is calculated using

$$ (1 - \eta_{\mathrm{B}}) \dot{E} = \int_{E_{\mathrm{min}}}^{E_{\mathrm{max}}} E \frac{d\dot{N}(E)}{dE} dE. \qquad (12) $$

We show the SED for PWN B0453-685 that can reasonably reproduce the observed spectrum shown in Figure 7.

To investigate the potential for a pulsar contribution to the Fermi-LAT data, we model the broadband spectrum again by adding a second emission component from the pulsar. Only the parameter *velpsr* is fixed to zero. In this case, we assume any Fermi-LAT pulsar flux can be described by a power law with an exponential cutoff of

$$ \frac{dN_{\gamma}}{dE} = N_0 E^{-\Gamma} \exp\left(-\frac{E}{E_{\mathrm{cut}}}\right), \qquad (13) $$
<figure><img src="image_9.png" /><figcaption>**Figure 7.** The best-fit SED assuming all Fermi-LAT emission is non-magnetospheric in origin (i.e., PWN only) obtained through the evolutionary model method described in Section 4.2. The colored points represent the values of observed data that the model used as comparison points for fitting and are the same values as those in the top panel of Figure 6. The small discontinuities in the SED between $\nu \sim 10^{25}$ and $10^{26}$ Hz are artifacts resulting from the specific numerical implementation of this code and are not physical.</figcaption></figure>

which is a common spectral characterization observed from gamma-ray pulsars (Abdo et al. 2013). We find that the pulsar together with its nebula can readily explain the lower-energy Fermi-LAT emission with a cutoff energy of $E_c = 3.21$ GeV and spectral index $\Gamma = 2.0$. The results are similar to the model presented for PWN Kes 75 and its central pulsar (Straal et al. 2023). Figure 8 displays both gamma-ray SEDs for the two considered cases where the Fermi-LAT emission is PWN only (left panel) and where there is both a PWN and pulsar contribution (right panel). If there is a pulsar contribution to the Fermi-LAT emission, it is likely to dominate for $E \lesssim 3$ GeV, whereas the PWN may only begin to dominate beyond this energy. We discuss the physical implications of the presented broadband models in the next section.

## 5. Discussion

The discovery of faint pointlike gamma-ray emission coincident with the SNR B0453-685 is presented. We can determine the physical properties of the host SNR and ambient medium from the broadband models presented in Sections 4.1 and 4.2 and compare them to the theoretical values expected for a middle-aged SNR in the Sedov–Taylor phase.

First, we can estimate the post-shock electron density assuming $\frac{n_e}{n_H} = 1.2$ and taking $n_H \sim 1.6 \text{ cm}^{-3}$ to find $n_e \sim 1.9 \text{ cm}^{-3}$. This result is consistent with prior works finding a range of values for a filling factor of $f, n_e/f \sim 1.5-8.0 \text{ cm}^{-3}$ (where a uniform density has $f=1$, Gaensler et al. 2003; Haberl et al. 2012; McEntaffer et al. 2012). The post-shock proton density $n_H = 1.6 \text{ cm}^{-3}$ is less than the average pre-shock LMC ISM density $n_0 \sim 2.0 \text{ cm}^{-3}$ (Kim et al. 2003). The total proton energy and the post-shock proton density characterizing pion decay emission are inversely proportional. If we assume $n_H$ is the expected shock-compressed LMC ISM density, then $n_H = 8.0 \text{ cm}^{-3}$. This would scale down the total energy in protons by a factor of $\frac{n_{h,LMC}}{n_{h,X-ray}} = 5$ to $W_p \sim 7.8 \times 10^{50}$ erg. This is a more reasonable particle energy, but both SNR models challenge the X-ray observations of the SNR shell,

which indicate an explosion energy as low as $E_{SN} \approx 10^{50}$ erg (Gaensler et al. 2003; Haberl et al. 2012).

The angular diameter of SNR B0453-685 in both radio and X-ray is $0.036$ (Figure 1), which corresponds to a shock radius of $R_s = 15.7$ pc at a distance of $d = 50$ kpc. We can evaluate the SNR age assuming it is in the Sedov–Taylor phase (Sedov 1959):

$$ \tau = \left( \frac{R_s}{2.3 \text{ pc}} \left( \frac{E}{10^{51} \text{ erg}} \right)^{\frac{1}{5}} \left( \frac{\rho_0}{10^{-24} \text{ g cm}^{-3}} \right)^{-\frac{1}{5}} \right)^{\frac{5}{2}} \text{100 yr.} \quad (14) $$

The SNR age estimates range between 13 kyr (Gaensler et al. 2003) using $E = 5 \times 10^{50}$ erg and $\rho_0 = m_H n_0 = 0.4 \times 10^{-24} \text{ g cm}^{-3}$ where $m_H$ is the mass of an H atom, and 15.2 kyr using $E = 7.6 \times 10^{50}$ erg and $\rho_0 = 0.3 \times 10^{-24} \text{ g cm}^{-3}$ (Haberl et al. 2012). McEntaffer et al. (2012) find the largest age estimates of $\tau \sim 17-23$ kyr using equilibrium shock velocity estimates $\sim 280-380 \text{ km s}^{-1}$. We adopt the SNR age reported in Gaensler et al. (2003), $\tau \sim 13$ kyr, which corresponds to a shock velocity $v_s = 478 \text{ km s}^{-1}$ from $v_s = \frac{2R_s}{5t}$. The age predicted from the evolutionary method in Section 4.2, $t \sim 14.3$ kyr, is in good agreement with prior work. The ambient proton density predicted in Section 4.2, $n_0 = 1.0 \text{ cm}^{-3}$, is somewhat higher than the values estimated in prior work (Gaensler et al. 2003; Haberl et al. 2012). In any case, the $n_0$ estimates are much lower than the average LMC ISM value $n_0 \sim 2 \text{ cm}^{-3}$ (Kim et al. 2003), and indicate that the ambient medium surrounding SNR B0453-685 may be less dense than the average LMC ISM. This is supported by Figure 1, left panel, where a possible density gradient decreasing from east to west is apparent. While $H\alpha$ emission is not a direct tracer for molecular material, it is a byproduct of SNRs interacting with molecular material (e.g., Winkler et al. 2014; Eagle et al. 2019). The lower ambient particle density estimate is also consistent with the observed faint SNR shell in radio and X-ray. It therefore seems unlikely for the SNR to be the gamma-ray origin, whether leptonic or hadronic.

We instead favor a model where the observed gamma rays are produced by an energetic neutron star and its resultant PWN, which can adequately describe the observed properties of this system as detailed in Section 4.2. The explosion energy predicted by the evolutionary model, $E = 5.2 \times 10^{50}$ erg, is very similar to that inferred by X-ray observations, $E \sim 5-7.6 \times 10^{50}$ erg (Gaensler et al. 2003; Haberl et al. 2012). Additionally, the magnetic field and total particle energy in the PWN from the evolutionary model are predicted to be $5.9 \mu\text{G}$ and $W_e = 5.4 \times 10^{48}$ erg, respectively, which is roughly consistent with the values implied by the NAIMA modeling in Section 4.1, $8.18 \mu\text{G}$ and $W_e = 2.9 \times 10^{49}$ erg. Lastly, one can estimate the gamma-ray efficiency $\eta = \frac{L_\gamma}{E}$ from the predicted current spin-down power of the central pulsar in the evolutionary model, $\dot{E} \sim 8.1 \times 10^{35}$ erg s$^{-1}$. For a 300 MeV–2 TeV gamma-ray source at $d = 50$ kpc, the gamma-ray luminosity is $L_\gamma = 2.6 \times 10^{35}$ erg s$^{-1}$, which corresponds to $\eta = 0.32$. This efficiency value is not uncommon for gamma-ray pulsars (e.g., Abdo et al. 2013), though it is a more compatible value to expect from evolved PWNe.

From the presented semi-analytic evolutionary models, we find the best representation of the data occurs with the supernova energy values $\sim 5 \times 10^{50}$ erg, $\sim 2.3 M_\odot$ for
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Left: the gamma-ray spectral evolutionary model assuming all Fermi-LAT emission is non-magnetospheric in origin (i.e., PWN only). Right: the gamma-ray spectral evolutionary model assuming magnetospheric contribution to the Fermi-LAT emission. The dotted line indicates the pulsar contribution and the dashed line indicates the PWN contribution. The colored points represent the values of observed data that the model used as comparison points for fitting and are the same values as those in the top panel of Figure 6. In both panels, the discontinuous spectral features beyond $\nu \sim 10^{25}$ Hz are numerical artifacts and can be ignored.</figcaption></figure>

supernova ejecta, and ~1.0 cm⁻³ for the density of the ISM (see Table 6). These values can then be used in combination with other models to survey the possible physical characteristics of the progenitor for SNR B0453-685. For example, a correlation reported in Ertl et al. (2020) found that the only supernovae that have an explosion energy of ~5 × 10⁵⁰ erg are those whose progenitors have a final helium core mass <3.5 M☉. Given an ejecta mass of ~2.3 M☉ from the presented evolutionary model, we calculate a neutron star mass $M_{NS} = 3.5-2.3 M☉ = 1.2 M☉$, which is reasonable (see, e.g., Kaper et al. 2006).

A core-collapse supernova progenitor cannot have an initial mass smaller than 8 M☉. We can use the known inverse correlation between the age and mass of a main-sequence star,

$$ \frac{t_{\text{MS}}}{t_{\text{Sun}}} \sim \left( \frac{M}{M_{\text{Sun}}} \right)^{-2.5} \qquad (15) $$

to get a maximum possible lifetime of $\tau \sim 20 \ 10^6$ yr for any supernova progenitor. A map by Harris & Zaritsky (2009) of the LMC with age and metallicity data distributions provides the age and metallicity distributions for the LMC regions closest to the location for B0453-685. By compiling the data in Harris & Zaritsky (2009), we can see that there was possibly a burst of star formation in those regions around the maximum possible lifetime estimate, as it contains many stars that are from $10^{6.8}$ (~6.3 million) to $10^{7.4}$ (~25 million) yr old. From this, the progenitor would have had a main-sequence lifetime comparable to the maximum possible lifetime for us to observe the SNR today. We can use Equation (15) to estimate the mass of the precursor star of B0453-685 to be between 11 and 19 M☉. However, as said above, the presented model predicts a pre-explosion helium core of 3.5 M☉, which does not reach the 11–19 M☉ dictated by the above analysis. The similarity between the inferred final core mass $M_{NS} = 1.2 M☉$ suggested by the presented modeling and the predicted pre-explosion mass $M_{\text{pre-explosion}} = 3.5 M☉$ from Ertl et al. (2020) implies that the progenitor lost its envelope before exploding.

If the models presented are correct, then there are two plausible ways to explain the loss of ~7.5–15.5 M☉ of material before exploding: an isolated star could have lost mass by way

of stellar wind, while a star that is part of a binary system could have transferred some of its mass to the other star. To account for stellar wind quantitatively, we looked at the model presented in Sukhbold et al. (2016) where it is shown that normal ejecta mass for a 10–15 M☉ star is 8–10 M☉, respectively. However, stellar wind can only account for up to 3 M☉ in mass loss for stars more massive than 15 M☉. Additionally, it is known that low-metallicity stars experience less mass loss (Heger et al. 2003), and the young stars in the LMC region of B0453-685 all have metallicity of ~0.008 Z☉. In summary, it seems plausible that the progenitor for B0453-685 was part of a binary star system.

## 6. Conclusions

We have reported the discovery of faint, pointlike gamma-ray emission by the Fermi-LAT that is coincident with the composite SNR B0453-685, located within the LMC. We provide a detailed multiwavelength analysis that is combined with two different broadband modeling techniques to explore the most likely origin of the observed gamma-ray emission. We compare the observed gamma-ray emission to the physical properties of SNR B0453-685 to determine that the association is probable. We then compare the physical implications and energetics from the best-fit broadband models to the theoretically expected values for such a system and find that the most plausible origin is the PWN within the middle-aged SNR B0453-685 and possibly a substantial pulsar contribution to the low-energy gamma-ray emission below $E < 5$ GeV. Theoretical expectation based on observational constraints and the inferred values from the best-fit models are consistent, despite assumptions about the SNR kinematics and environment in the evolutionary modeling method such as a spherically symmetric expansion into a homogeneous ISM density. The megaelectronvolt–gigaelectronvolt detection is too faint to attempt a pulsation search and the gamma-ray SED cannot rule out a pulsar component. We attempt to model the gamma-ray emission assuming both PWN and pulsar contributions and the results indicate that any pulsar gamma-ray signal is likely to be prominent below $E \le 5$ GeV, if present. Further work should explore the gamma-ray data, particularly for energies $E < 10$ GeV to investigate the potential for a pulsar
contribution as well as the possibilities for PWN and/or pulsar emission in the megaelectronvolt band for future megaelectron-volt space missions such as COSI<fnref n="19" /> and AMEGO.<fnref n="20" /> The IC emission spectra reported here may be even better constrained when combined with teraelectronvolt data from ground-based Cherenkov telescopes such as H.E.S.S. or the upcoming Cherenkov Telescope Array.<fnref n="21" />

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 from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France is gratefully acknowledged. This work was performed in part under DOE Contract DE-AC02-76SF00515. The work at NRL is supported by NASA.

**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

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

Matthew Kerr https://orcid.org/0000-0002-0893-4073

Patrick Slane https://orcid.org/0000-0002-6986-6756

Jean Ballet https://orcid.org/0000-0002-8784-2977

Fabio Acero https://orcid.org/0000-0002-6606-2816

Samayra Straal https://orcid.org/0000-0003-4136-7848

Marco Ajello https://orcid.org/0000-0002-6584-1703

## References

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

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

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

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

Ackermann, M., Ajello, M., Allafort, A., et al. 2013, Sci, 339, 807

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

Ackermann, M., Albert, A., Atwood, W. B., et al. 2016, A&A, 586, A71

Anders, E., & Grevesse, N. 1989, GeCoA, 53, 197

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

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

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

Balucinska-Church, M., & McCammon, D. 1992, ApJ, 400, 699

Blackburn, J. K. 1995, in ASP Conf. Ser. 77, Astronomical Data Analysis Software and Systems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Hayes (San Francisco, CA: ASP), 367

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

Castro, D., Slane, P., Patnaude, D. J., & Ellison, D. C. 2011, ApJ, 734, 85

Chen, Y., Jiang, B., Zhou, P., et al. 2014, in IAU Symp. 296, Supernova Environmental Impacts, ed. A. Ray & R. A. McCray (Cambridge: Cambridge Univ. Press), 170

Clementini, G., Gratton, R., Bragglia, A., et al. 2003, AJ, 125, 1309

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

Eagle, J., Marchesi, S., Castro, D., et al. 2019, ApJ, 870, 35

Ertl, T., Woosley, S. E., Sukhbold, T., & Janka, H.-T. 2020, ApJ, 890, 51

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

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

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., Haverkorn, M., Staveley-Smith, L., et al. 2005, in The Magnetized Plasma in Galaxy Evolution, ed. K. Chyży, K. Otmianowska-Mazur, M. Soida, & R.-J. Dettmar, 209, arXiv:astro-ph/0503371

Gaensler, B. M., Hendrick, S. P., Reynolds, S. P., & Borkowski, K. J. 2003, ApJL, 594, L111

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

Gaustad, J. E., Rosing, W., McCullough, P., & Van Buren, D. 2001, in ASP Conf. Ser. IAU Colloq. 183: Small Telescope Astronomy on Global 246, ed. B. Paczynski, W.-P. Chen, & C. Lemme (San Francisco, CA: ASP), 75

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

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

Gotthelf, E. V., Safi-Harb, S., Straal, S. M., & Gelfand, J. D. 2021, ApJ, 908, 212

Haberl, F., Filipovic, M. D., Bozzetto, L. M., et al. 2012, A&A, 543, A154

Harris, J., & Zaritsky, D. 2009, AJ, 138, 1243

Hattori, S., Straal, S. M., Zhang, E., et al. 2020, ApJ, 904, 32

Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003, ApJ, 591, 288

H.E.S.S. Collaboration, Abramowski, A., Acero, F., et al. 2012, A&A, 545, L2

Kaper, L., van der Meer, A., van Kerkwijk, M., & van den Heuvel, E. 2006, Msngr, 126, 27

Kim, S., Staveley-Smith, L., Dopita, M. A., et al. 2003, ApJS, 148, 473

Lopez, L. A., Ramirez-Ruiz, E., Badenes, C., et al. 2009, ApJL, 706, L106

Lopez, L. A., Ramirez-Ruiz, E., Huppenkothen, D., Badenes, C., & Pooley, D. A. 2011, ApJ, 732, 114

Malyshev, D., Cholis, I., & Gelfand, J. 2009, PhRvD, 80, 063005

Manchester, R. N., Fan, G., Lyne, A. G., Kaspi, V. M., & Crawford, F. 2006, ApJ, 649, 235

Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993

Mares, A., Lemoine-Goumard, M., Acero, F., et al. 2021, ApJ, 912, 158

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

McEntaffer, R. L., Brantseg, T., & Presley, M. 2012, ApJ, 756, 17

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

Reynolds, S. P. 2008, ARA&A, 46, 89

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

Russell, S. C., & Dopita, M. A. 1992, ApJ, 384, 508

Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics (New York: Academic Press)

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

Slane, P., Bykov, A., Ellison, D. C., Dubner, G., & Castro, D. 2015, SSRv, 188, 187

Straal, S. M., Gelfand, J. D., & Eagle, J. L. 2023, ApJ, 942, 103

Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H.-T. 2016, ApJ, 821, 38

Temim, T., Slane, P., Castro, D., et al. 2013, ApJ, 768, 61

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

Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487

Vink, J., & Laming, J. M. 2003, ApJ, 584, 758

Wakely, S. P., & Horan, D. 2008, ICRC (Mérida), 3, 1341

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

Winkler, P. F., Williams, B. J., Reynolds, S. P., et al. 2014, ApJ, 781, 65

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

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

<footnote n="19">¹⁹ https://cosi.ssl.berkeley.edu/</footnote>

<footnote n="20">²⁰ https://asd.gsfc.nasa.gov/amego/index.html</footnote>

<footnote n="21">²¹ https://www.cta-observatory.org/</footnote>