# TeV cosmic ray nuclei acceleration in shell-type supernova remnants with
hard $\gamma$-ray spectra

Houdun Zeng Key Laboratory of Dark Matter and Space Astronomy, Purple Mountain
Observatory, Chinese Academy of Sciences Nanjing 210034, People’s Republic of
China  Yuliang Xin  School of Physical Science and Technology, Southwest
Jiaotong University, Chengdu 610031, People’s Republic of China  Shuinai Zhang
Siming Liu Key Laboratory of Dark Matter and Space Astronomy, Purple Mountain
Observatory, Chinese Academy of Sciences Nanjing 210034, People’s Republic of
China

###### Abstract

The emission mechanism for hard $\gamma$-ray spectra from supernova remnants
(SNRs) is still a matter of debate. Recent multi-wavelength observations of
TeV source HESS J1912+101 show that it is associated with an SNR with an age
of $\sim 100$ kyrs, making it unlikely produce the TeV $\gamma$-ray emission
via leptonic processes. We analyzed Fermi observations of it and found an
extended source with a hard spectrum. HESS J1912+101 may represent a peculiar
stage of SNR evolution that dominates the acceleration of TeV cosmic rays. By
fitting the multi-wavelength spectra of 13 SNRs with hard GeV $\gamma$-ray
spectra with simple emission models with a density ratio of GeV electrons to
protons of $\sim 10^{-2}$, we obtain reasonable mean densities and magnetic
fields with a total energy of $\sim 10^{50}$ ergs for relativistic ions in
each SNR. Among these sources, only two of them, namely SN 1006 and RCW 86,
favor a leptonic origin for the $\gamma$-ray emission. The magnetic field
energy is found to be comparable to that of the accelerated relativistic ions
and their ratio has a tendency of increase with the age of SNRs. These results
suggest that TeV cosmic rays mainly originate from SNRs with hard $\gamma$-ray
spectra.

Galactic cosmic rays (567); Gamma-ray sources (633); Non-thermal radiation
sources (1119); Gamma-ray astronomy (628); Supernova remnants (1667)

##  1 Introduction

Although it is generally accepted that soft $\gamma$-ray spectra of supernova
remnants interacting with molecular clouds (SNRs) result from decay of
$\pi^{0}$ produced via inelastic collisions of high-energy ions with nuclei in
the background due to evolution of SNR shocks in a high-density environment
(Abdo et al., 2010; Giuliani et al., 2011; Zeng et al., 2019), the nature of
hard $\gamma$-ray spectra has been a subject of extensive investigations (Yuan
et al., 2012; Gabici & Aharonian, 2014; Zhang & Chen, 2016; H. E. S. S.
Collaboration et al., 2018a; Celli et al., 2019). In the leptonic scenario for
the $\gamma$-ray emission, the model parameters are well constrained and
appear to be consistent with expectation of diffusive shock particle
acceleration mechanism (Zhang & Liu, 2019a). Hadronic models require stronger
magnetic fields and less efficient electron acceleration (Butt et al., 2008).
These results have profound implications on the origin of cosmic rays (CRs),
especially those with energies lower than the CR spectral knee energy of $\sim
1$ PeV (Zhang & Liu, 2020).

Considering the anomalous CR spectra discovered with space measurements
(Aguilar et al., 2015a, b), Zhang et al. (2017) proposes that GeV CRs are
mostly accelerated in SNRs interacting with molecular clouds with relatively
lower shock speeds, giving rise to softer spectra, while higher energy CRs may
be attributed to particle acceleration in relatively younger SNRs with higher
shock speeds, leading to slightly harder high-energy spectra. In this
scenario, high-energy electron acceleration efficiency in young SNRs can
indeed be lower than that for GeV electrons in old SNRs (Zhang & Liu, 2019c),
reminiscence of the hadronic scenario for hard $\gamma$-ray spectra (Butt et
al., 2008). One of the challenges facing hadronic models is that ion spectra
need to cut off at tens of TeVs to account for the observed $\gamma$-ray
spectral cutoffs, suggesting that SNRs are not PeVatrons. Recent CR proton
spectral measurement by the DAMPE shows that there appears to be a spectral
hump at tens of TeVs (DAMPE Collaboration et al., 2019), which has been
attributed to a nearby SNR, such as Geminga (Qiao et al., 2019), implying
indeed that shocks of SNRs can only accelerate protons up to a few tens of TeV
(Lagage & Cesarsky, 1983; Bell et al., 2013). Observations of young SNR Cas A
and $\gamma$-Cygni SNR also imply a high-energy cutoff of the ion distribution
in the TeV energy range (Zhang & Liu, 2019b; Abeysekara et al., 2020; MAGIC
Collaboration et al., 2020).

HESS J1912+101 was first discovered in 2008 by the H.E.S.S. collaboration with
a shell structure (Aharonian et al., 2008). However its radio counterpart had
not been identified until very recently with polarization measurement (Reich &
Sun, 2019), implying presence of large scale strong magnetic field.
Interestingly, observations of molecular clouds in the direction of HESS
J1912+101 suggest that it is associated with an SNR with an age of $70-200$
kilo-years (kyrs) (Su et al., 2017), which is consistent with the
characteristic age of 170 kyrs for pulsar J1913+1011 inside it (Morris et al.,
2002). SNR G279.0+1.1 has similar properties with $\gamma$-ray emission up to
0.5 TeV detected recently (Araya, 2020). A recent population study of SNRs
shows that higher energy particles are preferentially accelerated in younger
SNRs with higher shock speeds (Zeng et al., 2019). In the absence of
continuous TeV particle acceleration in old SNRs, it is very challenging to
produce the TeV emission from HESS J1912+101 via leptonic processes.

In section 2, we show challenges to reproduce $\gamma$-ray emission from HESS
J1912+101 and SNR G279.0+1.1 via leptonic processes. Section 3 is dedicated to
the study of G296.5+10.0 for its hard $\gamma$-ray spectrum similar to HESS
J1912+101. In section 4 we fit the multi-wavelength spectra of another 10 SNRs
with hard $\gamma$-ray spectra in the hadronic and leptonic scenarios and
discuss the model implications. Our conclusions are drawn in section 6.

##  2 HESS J1912+101

Recent radio observations of HESS J1912+101 revealed strongly polarized
emission at 6 cm from the northeast half of the shell (Reich & Sun, 2019). Due
to strong radio emission from the surrounding, a shell structure cannot be
identified in the total intensity maps. The total polarized flux density at 6
cm is about 0.5$\pm$0.2 Jy, which can be considered as a lower limit to the
total flux density. Assuming a polarization fraction of 20%, Reich & Sun
(2019) obtained a total flux density of 2.5$\pm$ 1.0 Jy, which will be treated
as an upper limit in the following. The strong polarization of radio emission
implies the presence of large scale magnetic fields, which should also be
stronger than the typical value of $3\mu$G for the interstellar medium.

Via CO and HI observations, Su et al. (2017) showed that HESS J1912+101 is
likely associated with an old SNR with an age of $0.7-2.0\times 10^{5}$ years
at a distance of $\sim$ 4.1 kpc. The good correlation between the TeV emission
and the disturbed gas revealed by these observations makes them suggest a
hadronic origin for the $\gamma$-ray emission.

Figure 1: Electron cooling time due to synchrotron and inverse Compton
emission. The magnetic fields and properties of the background soft photons
are indicated. The grey band indicates the age range obtained via molecular
cloud observations (Su et al., 2017).

Recently Zeng et al. (2019) studied a sample of $\gamma$-ray SNRs and found
that in general high-energy particles in these sources have a broken power-law
distribution with the break energy and the low-energy spectral index
decreasing with the increase of SNR age. These results imply that higher
energy particles are mostly accelerated in younger SNRs with relatively higher
shock speeds. The acceleration of TeV particles quenches for SNR with an age
greater than 10 kyrs. This challenges a leptonic origin for the TeV emission
from HESS J1912+101.

Figure 1 shows the energy loss timescale of electrons due to the synchrotron
and inverse Compton (IC) processes. It can be seen that for a typical value of
the interstellar magnetic field, the maximum energy of electrons in HESS
J1912+101 should be about 10 TeV. The $\gamma$-ray produced via IC by such
electrons should cut off below 10 TeV. The left panel of Figure 2 shows the
evolution of electron distribution under the influence of energy loss due to
synchrotron and IC processes for an injected broken power-law spectrum with a
maximum energy of 1 PeV at the beginning. Compared to the age of this SNR, TeV
and higher energy particles are assumed to be accelerated instantaneously in
the early stage of the SNR evolution (Bell et al., 2013).

Figure 2: Left: evolution of electron distribution due to synchrotron and IC
losses. The initial distribution is a broken power-law with a high-energy
cutoff with parameters indicated in the figure. Right: evolution of the
corresponding emission spectra compared with the spectral energy distribution
(SED) of HESS J1912+101. The radio data are taken from Reich & Sun (2019). The
lower limit is for the polarized emission. The X-ray upper limit is from Chang
et al. (2008), and the blue data points are from the H.E.S.S. collaboration
(H. E. S. S. Collaboration et al., 2018b).

Figure 3: Left: $4^{\circ}\times 4^{\circ}$ TS map of photons above 1 GeV
around HESS J1912+101 after subtracting $\gamma$-ray emission from 4FGL
J1913.3+1019. The green crosses and circle represent the 4FGL sources and the
best-fit radius of the uniform disk is marked by the white circle. The cyan
contours show the TeV $\gamma$-ray emission of HESS J1912+101 (H. E. S. S.
Collaboration et al., 2018b). Right: The $\gamma$-ray SED of HESS J1912+101
(black dots) and 4FGL J1913.3+1019 (green dots). The blue dots are the HESS
data of HESS J1912+101 (H. E. S. S. Collaboration et al., 2018b).The gray
histogram denotes the TS value of HESS J1912+101 for each energy bin and the
arrows indicate the upper limits with 95% significance level.

Recent analyses of Fermi data by Zhang et al. (2020) uncovered a compact GeV
source with a soft spectrum within the shell of HESS J1912+101, and there
appears to be extended diffuse emission. Here we re-analyzed the Fermi-LAT
data around HESS J1912+101. In the region of HESS J1912+101, there are four
4FGL sources (4FGL J1913.3+1019, 4FGL J1914.7+1012, 4FGL J1911.7+1014, 4FGL
J1912.7+0957). Among them, 4FGL J1913.3+1019 is associated with PSR J1913+1011
(Morris et al., 2002; Smith et al., 2019; Zhang et al., 2020) and others are
unidentified sources with soft spectra. We removed these three unidentified
4FGL sources from the model file and treated them as parts of the diffuse
emission around HESS J1912+101. Using the data above 1 GeV with the Fermipy
package version 0.19.0 (Wood et al., 2017), we tested several spatial
templates for HESS J1912+101, including an uniform disk, a 2D Gaussian model,
and the HESS image. The uniform disk model (R.A., decl.=$288.262^{\circ}\pm
0.033^{\circ}$, $10.075^{\circ}\pm 0.037^{\circ}$ with
$\sigma=0.496^{\circ}\pm 0.035^{\circ}$) is favored, while the 2D Gaussian
template (R.A., decl.=$288.277^{\circ}\pm 0.034^{\circ}$, $10.135^{\circ}\pm
0.041^{\circ}$ with $\sigma=0.413^{\circ}\pm 0.036^{\circ}$) gives an equally
good representations of the data. We produced a TS map above 1 GeV after
subtracting the $\gamma$-ray emission from 4FGL J1913.3+1019, which is shown
in the left panel of Figure 3.

In the energy range of 1 GeV - 500 GeV, the TS value of HESS J1912+101 is
fitted to be 469.0, and the spectral index of a power-law model is
2.54$\pm$0.07. The integral photon flux from 1 GeV to 500 GeV is $(7.72\pm
0.44)\times 10^{-9}$ photon cm-2 s-1. For the point source 4FGL J1913.3+1019,
the TS value and the power-law spectral index are fitted to be 24.2 and
1.90$\pm$0.17. The $\gamma$-ray SEDs of HESS J1912+101 and 4FGL J1913.3+1019
shown in the right panel of Figure 3 were produced by dividing all data
between 1 GeV and 500 GeV into 10 bins with identical width on the logarithmic
of energy. And the 95% upper limits are for energy bins with the TS value of
HESS J1912+101 smaller than 4.0. The $\gamma$-ray SED shows a new spectral
component below 10 GeV, which may be attributed to shock interaction with
molecular clouds. The nature of this soft extended component will be explored
in a separate paper. Here we focus on spectral modeling above 10 GeV.

The right-panel of Figure 2 shows that even with a magnetic field of $3\
\mu$G, the age of HESS J1912+101 needs to be less than 10 kyrs to reproduce
the $\gamma$-ray spectrum, which contradicts molecular cloud observations (Su
et al., 2017) and properties of the associated pulsar (Smith et al., 2019).
Here the $\gamma$-ray is produced via IC processes (Jones, 1968), and besides
the cosmic microwave background radiation (CMB), we assume an infrared photon
background with $T=30$ K and an energy density of 1 eV cm-3 (Porter et al.,
2006). To facilitate model comparison, these values for background photons
will be used in this paper unless specified otherwise. Although the
$\gamma$-ray spectrum is flat, it implies an electron distribution with an
index of $\sim$3 in the leptonic scenario. Considering the low radio flux
density, it is evident that energetic electrons need to have a broken power
law distribution with a hard low-energy spectrum. Here the spectral index
$\alpha=1.7$ is indicated in the left panel of Figure 2. Beyond the break
energy of $E_{\rm br}=50$ GeV, the spectral index is 2.7.

Figure 4: Fit to the multi-wavelength SED of HESS J1912+101 in the hadronic
scenario for the $\gamma$-ray emission. The total energy of protons above 1
GeV is $10^{50}$ erg (Left) and $10^{49}$ erg (Right). The model parameters
are indicated on the Figures. The solid line is for $\gamma$-ray emission via
hadronic processes, while the dotted, dotted-dashed, and dashed lines are for
electron bremsstrahlung and IC of Infrared and CMB photons, respectively.

On the other hand, the $\gamma$-ray spectrum can be readily fitted in the
hadronic scenario. For the sake of simplicity, we assume a power-law
distribution with an identical index for electrons and ions.
$N(R_{i})=N_{0,i}R_{i}^{-\alpha}$exp$[-(E_{i}/E^{i}_{\rm cut})^{\delta}]$
where $R_{i}=p_{i}/q_{i}$ is the rigidity of the particle, $E$, $p$ and $q$
are the particle energy, momentum and charge, respectively, and $"i"$
represents different particle species. Considering that the flux ratio of TeV
CR electrons and protons is less than $0.1\%$ and high-energy electrons are
subject to radiative energy loss as they propagate in the Milky Way galaxy, we
fix the density ratio of electrons and protons at 1 GeV $K_{\rm ep}=N_{0,\rm
e}/N_{0,\rm p}=5\times 10^{-3}$. The total energy content of protons above 1
GeV ($W_{\rm p}$) determines the normalization of the particle distributions.
We will usually consider two cases: $W_{\rm p}=10^{49}$ and $W_{\rm
p}=10^{50}$. The mean background density $n_{\rm H}$ and magnetic field $B$,
and the spectral index $\alpha$ and cutoff energy of protons $E_{\rm cut}^{p}$
can be adjusted to fit the multi-wavelength SED. The cutoff energy of
electrons $E_{\rm cut}^{e}$ is obtained by requiring the corresponding
synchrotron energy loss time be equal to the age of the remnant and we assume
a super exponential high-energy cutoff for the electron distribution with
$\delta=2$ unless specified otherwise. The high-energy cutoff of ion
distribution is always exponential with $\delta=1$. When calculating the
$\gamma$-ray emission via the hadronic processes, we only consider protons and
contributions from other ions are approximated by multiplying the proton
produced $\gamma$-ray flux by a factor of 1.84 (Mori, 2009).

Figure 4 shows the results of the spectral fit and the model parameters are
listed in Table 2. When calculating the total energy of the magnetic field
$W_{B}$, we assume a uniform magnetic field with a volume filling factor of 1.
This magnetic field energy therefore should be considered as an upper limit.
Although both set of parameters give good fit to the radio and $\gamma$-ray
spectra, the model with a weaker magnetic field and therefore more energy in
energetic particles (left panel of Figure 4) is favored for the more
reasonable value of the magnetic field energy. Moreover, the case with a
strong magnetic field has a synchrotron spectrum cutting off in the radio band
(the right panel of Figure 4), which appears to be too low. Of course, the
magnetic field is not well constrained. For the given radio emission, one can
always compensate the increase in the magnetic field with a decrease in the
energetic particle energy.

Araya (2020) recently carried out a detailed analyses of G279.0$+$1.1, a huge
SNR with a radius of $\sim 40$ pc. It has a very hard $\gamma$-ray spectrum
with a spectral index of $1.86\pm 0.09$. Although it has not been detected in
the TeV range, the $\gamma$-ray spectrum obtained with the Fermi-LAT extends
to $0.5$ TeV without any indication of a spectral softening at high energies.
The age of this remnant is greater than $100$ kyr. Figure 5 shows our spectral
fits in the hadronic scenario for the $\gamma$-ray emission. It is interesting
to note that there appears to be a low energy spectral component near 1 GeV,
reminiscence of the spectral component below 10 GeV for HESS J1912+101.
Compared with HESS J1912+101, SNR G279.0$+$1.1 has a similar $\gamma$-ray
luminosity, a few times higher radio luminosity and radius, leading to a
higher magnetic field and a large magnetic energy. However, the small volume
filling factor of the radio emission region can reduce the magnetic energy
significantly.

Figure 5: Same as Figure 4 but for G279.0$+$1.1. The total energy of protons
above 1 GeV is $10^{50}$ erg (Left) and $5\times 10^{49}$ erg (Right).

##  3 G296.5+10.0

G296.5$+$10.0 is a bilateral morphology SNR in radio and X-rays with an
angular extension of 90 arcmin $\times$ 65 arcmin, and the distance is
estimated as 2.1 kpc (Giacani et al., 2000). It has relatively bright radio
emission with a typical spectral index of $-0.5$ (Milne & Haynes, 1994). Due
to its large size, Chandra, XMM-Newton, and Suzaku have not mapped the entire
SNR, and only ROSAT PSD data provided the thermal spectrum for the whole SNR
(Kellett et al., 1987). Nevertheless, five XMM-Newton observations taken in
2016 roughly cover the bright limbs of G296.5+10.0 (PI: Brian Whillianms).

Figure 6: The left panel shows X-ray images of the northeast part of SNR
G296.5+10.0, where the red, green, and blue colors represent emission in the
0.2-0.5 keV, 0.5-1.0 keV, and 1.0-2.0 keV bands, respectively. The color bars
in the bottom are in units of photons cm-2 s-1. The blue circle with a radius
of 5.15′′ is the region for spectral extraction, while the other circles and
the square are for the ‘local background’ regions. The green box is for the PN
data, and the white and the yellow circles are for the MOS1 and MOS2 data,
respectively. A few point sources are marked out with red. The middle panel
shows the source-region spectrum and the scaled ‘local background’ spectrum
(according to the effective area) of the PN data, with the model of unresolved
AGN emission in the cosmic X-ray background over-plotted. The right panel
shows the XMM-Newton PN and MOS spectra with their ‘local background’
subtracted, and the best-fit models are shown in the red color. The black
curve, as modeled by an absorbed power law, represents the non-thermal
emission, which is likely caused by the cosmic X-ray background.

Figure 7: TS maps of a $3^{\circ}\times 3^{\circ}$ region around G296.5+10.0.
The top two panels is for photons in the range of 1 GeV - 1 TeV, and the
bottom two panels is for photons from 5 GeV to 1 TeV. The magenta dashed
circle and the white solid circle represent the size of spatial template of
G296.5+10.0 in 4FGL and this work by Fermi-py, respectively. The white cross
show the position of the radio-quiet X-ray-emitting neutron star (1E
1207.4-5209) in it (Helfand & Becker, 1984; Zavlin et al., 2000). The green
(left) and cyan (right) contours describe the radio and X-ray emission of
G296.5+10.0, respectively.

Here we analyze the XMM-Newton data to estimate the flux of non-thermal
emission. This emission could be fairly dim as it was not reported in the
ROSAT data, and even the residual contamination from the soft proton flares
may affect the detection. As a result, only the ID:0781720101 observation is
used, which didn’t suffer apparent soft proton flares during its 28 ks
exposure. The data reduction employs the standard procedure in the XMM-Newton
Science Analysis System (SAS; version: 15.0.0) and the Extended Source
Analysis Software (ESAS) package. The diagnostic file created by the ‘pn-
filter’ shows a steady light curve with PN count rate of $\sim$3.2 count/s in
this observation, suggesting the absence of soft proton flares. We use the
‘eimageget’ script to generate the PN and MOS images, and the ‘eimagecombine’
to combine them and produce a background-subtracted, vignetting-corrected, and
smoothed image, where the background is based on the filter wheel closed
images. This image (the left panel of Figure 6) covers the northeast part of
the SNR.

The PN and MOS spectra are extracted from a circular region with a radius of
5.15′, well confined within the central CCD of the MOS1 and MOS2. A few bright
point sources are manually removed when making the mask. The ‘mos/pn-spectra”
scripts are used to create spectra, and the ‘mos/pn-back’ scripts to generate
model quiescent particle background (QPB) spectra. Using the same method we
also produce spectra from ‘local background’ regions at the northeast corner
of the field of view, where the SNR emission seems negligible. These regions
are slightly different for the MOS1, MOS2, and PN as shown in the left panel
of Figure 6, due to differences in the effective boundaries of different
detectors. We present the PN source spectrum and the ‘local background’
spectrum in the middle panel of Figure 6 with their QPB spectra subtracted,
and compared with the model of the unresolved AGN emission in the cosmic X-ray
background (De Luca & Molendi, 2004; Kuntz & Snowden, 2008). The model is an
absorbed power law component with a spectral index of 1.46 and a normalization
of $\sim 11.6\,\rm photons$ $\rm cm^{-2}s^{-1}sr^{-1}keV^{-1}$ at 1 keV, and
the absorption is determined by the Galactic HI column density towards this
SNR of about $1.1\times 10^{21}$ cm-2 (HI4PI Collaboration et al., 2016). In
the hard X-ray band (2-7 keV), the PN source spectrum is slightly under the
cosmic X-ray background curve, suggesting absence of prominent non-thermal
emission.

Nevertheless, to estimate an upper limit of the non-thermal emission, we fit
the PN and MOS spectra jointly, but using the ‘local background’ spectra as
the background spectra. A single temperature APEC model represents the SNR
thermal emission, and a power law model is for the non-thermal emission. The
redshift is set to zero, and the Galactic HI column density is assigned to be
the maximum column density for the foreground absorption. The best-fit spectra
are shown in the right panel of Figure 6. The best-fit temperature of the
thermal emission is 0.15 keV, and the metal abundances are slightly sub-solar
(Table 1). The non-thermal flux in the range of 0.5–10 keV is $6.0(\pm
1.5)\times 10^{-13}$ erg cm-2 s-1 , with a photon index of $1.4\pm 0.4$. Since
that the ‘local background’ spectra do not have a good statistic in the hard
X-ray band, this component is likely from the unresolved AGN emission. Note
that the photon index of 1.4 for SNR is too hard if the hard tail X-ray
component is synchrotron X-rays. We will treat this fitted flux as an upper
limit to nonthermal emission from the SNR, which in turn suggests an upper
limit of the non-thermal flux of $\sim 1.4\times 10^{-11}$ erg cm-2 s-1 for
the entire SNR region (Kellett et al., 1987).

Table 1: The best-fit results for the XMM-Newton data of G296.5+10.0

Model: TBabs*(APEC+Powerlaw)
---
TBabs | $n_{\rm H}$ |  $(1.0-1.1)\times 10^{21}$ cm-2
APEC |  $\eta_{\mathrm{apec}}$**Normalization parameter of the APEC model has the physical meaning of $\frac{10^{-14}}{4\pi[D_{\rm A}(1+z)^{2}]}\int n_{\rm e}n_{\rm H}dV$, where $D_{\rm A}$ is the angular diameter distance to the source (cm) and $z$ is the redshift, $n_{\rm e}$ and $n_{\rm H}$ are the electron and the hydrogen densities (cm-3), and $V$ is the volume.  | kT [keV] | C | N | O | Ne, Mg, Si, S | Fe, Ni
| $0.04\pm 0.01$ | $0.150\pm 0.001$ | $<0.6$ | $0.2\pm 0.1$ | $0.3\pm 0.1$ | $0.7\pm 0.2$ | $1.1\pm 0.3$
Powerlaw | $flux_{(0.5-10\,\mathrm{keV})}$ |  $6.0(\pm 1.5)\times 10^{-13}$ erg cm-2 s-1
| Photon Index | $1.4\pm 0.4$
$\chi^{2}$/$d.o.f.$ | 967/480 $\sim$ 2

Note. — Compared to the MOS spectra, the PN spectrum may have some systematic
deficit in flux. The total flux of the PN is allowed to vary during the joint
fitting, and the best-fit value here is $0.84\pm 0.01$.

Using 52 months of Pass 7 data recorded by Fermi-LAT, Araya (2013) detected
extended $\gamma$-ray emission toward the region of SNR G296.5+10.0 with a
$\sim$ 5$\sigma$ confidence level. The $\gamma$-ray spectrum of it can be
fitted by a power-law with an index of 1.85$\pm$0.13. In the fourth Fermi-LAT
source catalog (4FGL; Abdollahi et al., 2020), G296.5+10.0 corresponds to the
extended source 4FGL J1208.5-5243e, and the spatial template of it is adopted
to be an uniform disk centered at (R.A.$=182.13^{\circ}$,
Dec.$=-52.73^{\circ}$) with an radius of 0.76∘.

Figure 8: The $\gamma$-ray SED of G296.5+10.0. The red solid line shows the
best-fitting power-law spectrum in the energy range of 100 MeV - 1 TeV, and
the red dashed lines show its 1$\sigma$ statistic error. The gray histogram
denotes the TS value for each energy bin and the arrows indicate the upper
limits with 95% significance level.

Figure 9: Same as Figure 4 but for G296.5+10.0.

With the latest Fermi-LAT Pass 8 data collected from 2008 August 4 to 2018
August 4, we updated the $\gamma$-ray morphology and spectrum of G296.5+10.0.
Figure 7 shows the TS maps of a $3^{\circ}\times 3^{\circ}$ region around
G296.5+10.0. In each panel, the magenta dashed circle is the size of spatial
template of G296.5+10.0 in 4FGL (Abdollahi et al., 2020). As can be seen, the
$\gamma$-ray emission of G296.5+10.0 has a significant deviation with this
spatial template. Therefore, we use the Fermi-py tool to refit the data and
get the updated centered position and radius of the uniform disk, which is
shown as the white solid circle in Figure 7. The best-fit position and radius
of the uniform disk are (R.A.$=182.11^{\circ}$, Dec.$=-52.38^{\circ}$) and
0.79∘, respectively. With this spatial template, the TS value of G296.5+10.0
is fitted to be 203.68, and the spectrum can be well fitted by a power-law
with an index of 1.90$\pm$0.06. In the energy range from 1 GeV to 1 TeV, the
integral photon flux is $(1.19\pm 0.11)\times 10^{-9}$ photon cm-2 s-1 with
statistical error only. Moreover, to derive the $\gamma$-ray SED of
G296.5+10.0, all of the data from 100 MeV to 1 TeV are divided to be 14 equal
logarithmic bins. And for any energy bin with TS value of G296.5+10.0 smaller
than 5.0, an upper limit at 95% confidence level is calculated. The SED and
fit results are shown in Figure 8.

G296.5$+$10.0 has an age of $\sim 10$ kyrs (Vasisht et al., 1997) and a GeV
$\gamma$-ray spectrum similar to that of HESS J1912+101. Adopting the hadronic
model for HESS J1912+101, one can fit the radio and $\gamma$-ray spectra of
G296.5$+$10.0. Due to the absence of TeV observations, the high-energy cutoff
of the proton distribution is not well constrained and we adopt a value of 70
TeV derived from the spectral fit of G296.5$+$10.0. Figure 9 shows the
spectral fit, and the model parameters are listed in Table 2, which are very
similar to those for HESS J1912+101 except for a slightly stronger magnetic
field. As we will see below, this is due to the higher radio to $\gamma$-ray
flux ratio of G296.5$+$10.0 than that of HESS J1912+101. Similarly, we favor
the model with a weaker magnetic field (left panel of Figure 9).

##  4 Hadronic and/or Leptonic models for hard GeV spectra of other 10 SNRs

Figure 10: The multi-wavelength spectral data of 13 SNRs with hard
$\gamma$-ray spectra. The $\gamma$-ray spectra are fitted with a hadronic
model with the normalization of individual spectrum as free parameters. The
model assumes that protons have a single power-law energy distribution with an
exponential high-energy cutoff. Note that the TeV spectra of G78.2$+$2.1
(HAWC) and N132D (HESS) cut off at relatively lower energies, and the soft
spectral component of GeV of HESS 1912$+$101 may be from other contributors,
and are not considered in SED fitting. The best-fit model parameters are
indicated on the Figure. References for the observational data are as follows.
RX J0852.0$-$4622: radio (Duncan & Green, 2000), GeV (Tanaka et al., 2011),
X-ray (Aharonian et al., 2007), TeV (H. E. S. S. Collaboration et al., 2018c);
RX J1713.7$-$3946: radio (Lazendic et al., 2004), X-ray (Tanaka et al., 2008),
GeV and TeV (H. E. S. S. Collaboration et al., 2018a); HESS J1731$-$347: radio
(Tian et al., 2008), GeV (Condon et al., 2017; Guo et al., 2018), X-ray
(Doroshenko et al., 2017), TeV (H. E. S. S. Collaboration et al., 2011a); RCW
86: radio (Clark et al., 1975; Lemoine-Goumard et al., 2012), X-ray (Lemoine-
Goumard et al., 2012), GeV (Ajello et al., 2016), TeV (H. E. S. S.
Collaboration et al., 2018d); SN 1006: radio Dyer et al. (2009), X-ray (Bamba
et al., 2008), GeV (Condon et al., 2017), TeV (Acero et al., 2010);
G$150.3+4.5$: radio (Gerbrandt et al., 2014), X-ray and GeV (Devin et al.,
2020); G$296.5+10.0$: radio (Milne & Haynes, 1994), GeV (this work), HESS
J1534$-$571: radio (Maxted et al., 2018), GeV (Araya, 2017), X-ary and TeV (H.
E. S. S. Collaboration et al., 2018b); RCW 103: radio (Dickel et al., 1996),
GeV (Xing et al., 2014); G78.2$+$2.1: radio (Wendker et al., 1991; Zhang et
al., 1997; Kothes et al., 2006; Gao et al., 2011), X-ray (Leahy et al., 2013),
GeV (Abeysekara et al., 2018) and TeV (Fleischhack, 2019); G279.0$+$1.1: radio
(Woermann & Jonas, 1988; Duncan et al., 1995), GeV Araya (2020); N132D: radio
(Dickel & Milne, 1995), X-ray (Hughes et al., 1998; Bamba et al., 2018), GeV
(Xin et al. 2020, in preparation) and TeV (H. E. S. S. Collaboration et al.,
2015).

We now extend the above study to all SNRs with hard $\gamma$-ray spectra.
Figure 10 shows the multi-wavelength non-thermal spectra of 12 SNRs with hard
GeV spectra. The spectra have been normalized at $\gamma$-ray band by fitting
the $\gamma$-ray spectra with a hadronic model. The model assumes that protons
have a single power law with an exponential high energy cutoff. The
normalization of the SED of each source is adjusted to minimize the $\chi^{2}$
of the $\gamma$-ray data of all sources. We then have an index of 1.74 and
$E_{\rm cut}^{p}=47$ TeV. Then we may classify these sources based on their
radio and/or X-ray spectra.

Based on the radio flux, these SNRs may be divided into three categories. SN
1006, RCW 86, G296.5+10.0, G78.2+2.1 and N132D have strong radio emission,
while the radio emission from G150.3+4.5 and HESS J1534-571 are very weak.
Their normalized radio flux densities can differ by about two orders of
magnitudes. The famous TeV bright SNRs RX J1713-3946, RX J0852.0-4622, HESS
J1731-347, and HESS J1912+101, RCW 103 have normalized radio flux densities
between the above two categories. In Table 2, we use double horizontal lines
to separate these categories and give the estimated age, distance, radius of
these SNRs and the related references. Non-thermal X-ray emission is detected
from five relatively young SNRs: SN 1006, RCW 86, RX J1713-3946, RX
J0852.0-4622, and HESS J1731-347. As will be shown below, a broken power-law
spectrum is needed to explain their multi-wavelength non-thermal emission
spectra, while for other sources, the above hadronic model for HESS J1912+101
and G296.5+10.0 with a single power-law high-energy particle spectrum is
sufficient (See Table 2).

Table 2: Physical and fitting parameters for our sample.

Source | Age | Distance | Radius | $W_{\rm p}$ | $n_{\rm H}$ | $\alpha$ | $E^{\rm p}_{\rm cut}$ | $\Delta\alpha$ | $\delta$ | $E^{\rm e}_{\rm br}$ | $E^{\rm e}_{\rm cut}$ | $B$ | $W_{B}$ | $W_{B}/W_{\rm e}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
| kyr | kpc | pc |  $10^{49}$ erg | cm-3 |  | TeV |  |  | GeV | TeV |  $\mu$G | 1049 erg | 102
G150.3$+$04.5 | 6.0 | 1.0 | 24 | 1.0 | 4.2 | 1.75 | 70 | $-$ | 2.0 | $-$ | 2.0a | 32 | $6.8$ | 34.6
|  |  |  | 10.0 | 0.48 | 1.75 | 70 | $-$ | 2.0 | $-$ | 0.2 | 6.0 | $0.24$ | 0.26
HESS J1534$-$571 | 10.0 | 3.5 | 22.4 | 1.0 | 12 | 1.45 | 20 | $-$ | 2.0 | $-$ | 0.14a | 300 | 496 | $7.2\times 10^{4}$
G323.7$-$01.0 |  |  |  | 10.0 | 1.2 | 1.45 | 20 | $-$ | 2.0 | $-$ | 0.59a | 46 | 11.6 | ${\bf 16}$
HESS J1912$+$101 | 170 | 4.1 | 15 | 1.0 | 27 | 1.9 | 70 | $-$ | 2.0 | $-$ | 0.0051a | 120 | 23.9 | $6.0\times 10^{2}$
G044.5$-$00.2 |  |  |  | 10.0 | 2.7 | 1.9 | 70 | $-$ | 2.0 | $-$ | 0.1225a | 24.5 | 1.0 | 0.63
G279.0$+$1.1 | $100$ | 3.0 | 41.5 | 5.0 | 5.2 | 1.9 | 40 | $-$ | 2.0 | $-$ | 0.006a | 145 | $7.3\times 10^{2}$ | $3.0\times 10^{3}$
|  |  |  | 10.0 | 2.6 | 1.9 | 40 | $-$ | 2.0 | $-$ | 0.0186a | 82 | $2.3\times 10^{2}$ | ${\bf 256}$
RCW 103 | 4.4 | 3.3 | 4.8 | 1.0 | 65 | 2.05 | 70 | $-$ | 2.0 | $-$ | 0.56a | 226 | 2.8 | 12
G332.4$-$00.4 |  |  |  | 10.0 | 6.5 | 2.05 | 70 | $-$ | 2.0 | $-$ | 1.14a | 50 | 0.14 | 0.034
RX J1713.7$-$3946 | 1.6 | 1.0 | 8.7 | 0.6 | 11.7 | 1.6 | 40 | 1.0 | 2.0 | 25.8a | 13.5 | 550 | 97.4 | $3.9\times 10^{3}$
|  |  |  | 5.0 | 1.2 | 1.55 | 40 | 1.75 | 2.0 | 390a | 40 | 142 | 6.5 | ${14}$
G347.3$-$00.5 |  |  |  | 10.0 | 0.6 | 1.6 | 40 | 1.0 | 2.0 | 82 | 38 | 62 | 1.24 | ${\bf 1.73}$
(Leptonic) |  |  |  | 7.7 | 0.1 | 1.89 | 70 | 1.0 | 2.0 | 1800 | 70 | 22 | 0.16 | 0.055
RX J0852.0$-$4622 | 2.7 | 0.75 | 13 | 1.0 | 5.5 | 1.75 | 70 | 1.0 | 2.0 | 28 | 20 | 165 | 29.9 | $3.9\times 10^{2}$
|  |  |  | 1.0 | 5.5 | 1.75 | 70 | 1.4 | 2.0 | 170a | 20 | 165 | 29.9 | $2.4\times 10^{2}$
|  |  |  | 10.0 | 0.45 | 1.5 | 35 | 2.32 | 2.0 | 460a | 30 | 100 | 11 | ${14}$
G266.2$-$01.2 |  |  |  | 10.0 | 0.55 | 1.75 | 70 | 1.0 | 2.0 | 55 | 50 | 31 | 1.06 | ${\bf 1.08}$
(Leptonic) |  |  |  | 350 | 0.001 | 1.5 | 70 | $1.0$ | 2.0 | $20$ | 70 | 8 | 0.7 | 0.012
(Leptonic) |  |  |  | 6.0 | 0.2 | 2.2 | 70 | $-$ | 1.0 | $-$ | 22 | 11 | 0.13 | 0.037
HESS J1731$-$347 | 2.5 | 3.2 | 14 | 1.0 | 7.5 | 1.5 | 30 | 1.0 | 2.0 | 38 | 6.5 | 550 | 403 | $1.1\times 10^{4}$
|  |  |  | 1.0 | 7.5 | 1.5 | 30 | 0.85 | 2.0 | 16.5a | 6.4 | 550 | 403 | $1.5\times 10^{4}$
|  |  |  | 10.0 | 0.6 | 1.5 | 30 | 1.82 | 2.0 | 692a | 25 | 85 | 9.6 | ${8.3}$
G353.6$-$00.7 |  |  |  | 10.0 | 0.75 | 1.5 | 30 | 1.0 | 2.0 | 100 | 17 | 85 | 9.6 | ${\bf 16.2}$
(Leptonic) |  |  |  | 55 | 0.05 | 1.5 | 70 | 1.0 | 2.0 | 170 | 25 | 32 | 1.05 | 0.29
(Leptonic) |  |  |  | 5.5 | 0.05 | 1.95 | 70 | $-$ | 1.0 | $-$ | 13 | 23 | 0.71 | 0.30
G296.5$+$10.0 | 10.0 | 2.1 | 24.1 | 1.0 | 4.2 | 1.9 | 70 | $-$ | 2.0 | $-$ | 0.0097a | 360 | $8.9\times 10^{2}$ | $1.5\times 10^{4}$
|  |  |  | 10.0 | 0.35 | 1.9 | 70 | $-$ | 2.0 | $-$ | 0.26a | 70 | 33.7 | ${\bf 18}$
SN 1006 | 1.0 | 2.2 | 9.6 | 1.0 | 1.0 | 1.90 | 70 | 1.0 | 1.0 | 386a | 8.1 | 180 | 14 | $52$
G327.6$+$14.6(Lep) |  |  |  | 1.5 | 0.1 | 2.1 | 70 | $-$ | 1.0 | $-$ | 7.3 | 65 | 1.8 | ${\bf 2.4}$
RCW 86 | 1.8 | 2.5 | 15.3 | 1.0 | 6.5 | 1.45 | 20 | 0.85 | 1.0 | 0.048a | 0.85 | $1.2\times 10^{4}$ | $2.5\times 10^{5}$ | $9.0\times 10^{8}$
G315.4$-$02.3(Lep) |  |  |  | 10 | 0.05 | 2.25 | 70 | 1.0 | 1.0 | 7200a | 30 | 31 | 1.68 | 0.27
Gamma Cygni | 8.25 | 2.0 | 17 | 1.0 | 29 | 2.00 | 10 | $-$ | 2.0 | $-$ | 0.0036a | 650 | $1.0\times 10^{3}$ | $1.6\times 10^{4}$
G78.2$+$01.2 |  |  |  | 10 | 2.9 | 2.00 | 10 | $-$ | 2.0 | $-$ | 0.077a | 140 | 47.2 | ${\bf 18}$
N132D | 2.5 | 50 | 11.4 | 10.0 | 32 | 2.10 | 70 | $-$ | 2.0 | $-$ | 0.028a | 423 | $1.3\times 10^{2}$ | 56
N132D | 2.5 | 50 | 11.4 | 50.0 | 6.4 | 2.10 | 70 | $-$ | 2.0 | $-$ | 0.23a | 148 | 15.9 | 0.89

$a:$ Determined by requiring the synchrotron energy loss time being equal to
the SNR age;

Note. — References of physical parameters $-$ HESS J1912$+$101 (Su et al.,
2017; Zhang et al., 2020); RCW 103 (Reynoso et al., 2004; Braun et al., 2019);
RX J1713.7$-$3946 (Fukui et al., 2003; Tsuji & Uchiyama, 2016); RX
J0852.0$-$4622 (Katsuda et al., 2008); HESS J1731$-$347 (Tian et al., 2008; H.
E. S. S. Collaboration et al., 2011b); G150.3$+$4.5 (Cohen, 2016; Devin et
al., 2020); HESS J1534$-$571 (Maxted et al., 2018); G279.0$+$1.1 (Duncan et
al., 1995; Araya, 2020); G296.5$+$10 (Vasisht et al., 1997; Giacani et al.,
2000); SNR 1006 (Winkler et al., 2003; Katsuda et al., 2009); RCW 86 (Bocchino
et al., 2000; Helder et al., 2013); $\gamma$ Cygni (Higgs, 1977; Leahy et al.,
2013); N132D (Dickel & Milne, 1995; Vogt & Dopita, 2011).

Figure 11: Same as Figure 4 but for G150.3$+$4.5 (upper) and HESS J1534$-$571
(lower) with very weak radio emission.

For the two SNRs with very weak radio emission, G150.3$+$4.5 is similar to
G296.5+10.0 in the sense that there is no TeV data. We therefore set the high-
energy cutoff of the proton distribution at 70 TeV. The spectral fits are
shown in the upper panels of Figure 11. We favor the model with a low value of
$10^{49}$ erg for the total proton energy $W_{\rm p}$ (left panel) with a
magnetic field of $32\ \mu$G. An even lower value of $W_{\rm p}$ is disfavored
since it implies an even stronger magnetic field and therefore an energy ratio
of the magnetic field to the high-energy electrons greater than 3500. For
higher values of the total proton energy, the electron cutoff energy needs to
be lower than that determined by requiring the synchrotron energy loss time
being equal to the age of the SNR. Otherwise the IC emission will dominate the
$\gamma$-ray emission with a harder spectrum than the observed one. The right
panel shows such a case with the electron cutoff energy of 200 GeV to suppress
contributions to the $\gamma$-ray via the IC processes. Equating the
synchrotron energy loss time to the SNR age will lead to a cutoff energy of
electron $E^{\rm e}_{\rm cut}$ greater than 50 TeV for the relatively weaker
magnetic field of $6\,\mu$G.

We notice that there are some degeneracies among $W_{\rm p}$, $n_{\rm H}$, $B$
and $K_{\rm ep}$. The product of $W_{\rm p}$ and $n_{\rm H}$ is determined by
the $\gamma$-ray spectrum. The total proton energy discussed above should be
re-scaled by the actual value of the mean number density of the background
$n_{\rm H}$. We fix $K_{\rm ep}$ to $5\times 10^{-3}$ in the spectral fit
above. For a given radio flux and $W_{\rm p}$, an increase of $K_{\rm ep}$
will lead to more energetic electrons and a weaker magnetic field.

The lower panels of Figure 11 show the spectral fits for HESS J1534$-$571\.
The cutoff energy of the proton distribution for HESS J1534$-$571 is well
constrained by TeV observations. However, Table 2 shows that the cutoff energy
of $20$ TeV is much lower than the value of 70 TeV for the three sources
studied above. Similar to G296.5+10.0 and HESS J1912+101, we favor the model
with a higher value of $W_{\rm p}$ (left panel) for the weaker magnetic field
and higher synchrotron cutoff energy. The ratio of $W_{B}/W_{\rm e}=1600$ is
also more reasonable. However, the magnetic field is not as well constrained
as for G150.3$+$4.5. We also notice that with a spectral index of 1.45, HESS
J1534$-$571 has a much harder energetic particle spectrum than other sources.
Such a hard spectrum is needed to fit the $\gamma$-ray spectrum. More radio
flux density measurements are needed to test this model.

Figure 12: Same as Figure 4 but for RCW 103.

Among the 5 SNRs with intermediate radio emission, RCW 103 does not have non-
thermal X-ray emission, similar to HESS J1912$+$101\. There is also no TeV
data, we then fix the proton cutoff energy at $70$ TeV. Although RCW 103 is
relatively young with an age of 4.4 kyrs, radio and GeV observations lead to a
high-energy particle spectrum slightly softer than that for HESS J1912$+$101\.
The GeV flux of RCW 103 is actually about 3 times higher than that of HESS
J1912$+$101\. Since the distance to the two source is comparable, the
$\gamma$-ray luminosity of RCW 103 of about 2 times higher. For the same
proton energy $W_{\rm p}$ and comparable radio flux densities, RCW 103
therefore has a higher background density $n_{\rm H}$ and stronger magnetic
field. However, the radius of RCW 103 is more than 3 times smaller than HESS
J1912$+$101, leading to a factor of 30 difference in the volume. We therefore
have a much lower value of $W_{B}/W_{\rm e}$. Even the case with a $W_{\rm p}$
of $10^{50}$ erg (left panel of Figure 12) has a very low value of 0.034 for
$W_{B}/W_{\rm e}$, we still favor this model for its relatively higher cutoff
energy of the synchrotron spectrum. To increase the value of $W_{B}/W_{\rm
e}$, one needs to consider lower values for $W_{\rm p}$ as discussed above
(See right panel of Figure 12).

Figure 13: Several spectral fits to the SED of RX J17137.7$-$3946\. The upper
left is our favored model. The upper right corresponds to the leptonic
scenario for the $\gamma$-ray emission. Although the model shown in the lower
left has fewer parameters, the magnetic field appears to be too strong. The
model in the lower right panel has a relatively weaker magnetic field but
gives a poorer fit to the SED.

The other 3 SNRs with intermediate radio emission RX J1713.7$-$3946, RX
J0852$-$4622, and HESS J1731$-$347 have been studied extensively for their
prominent non-thermal X-rays and TeV emission. There are still debates on the
nature of the $\gamma$-ray emission. Detailed TeV observations of RX
J17137.7$-$3946 have shown that a broken power-law distribution is needed to
fit the $\gamma$-ray spectrum in both the leptonic and hadronic scenarios for
the $\gamma$-ray emission (H. E. S. S. Collaboration et al., 2018a). In
general, for sources with strong non-thermal X-ray emission, a single power
law particle distribution will lead to a poor fit to the multi-wavelength SED.
But for the sake of simplicity and considering the overall quality of
$\gamma$-ray data, we will still adopt a single power-law distribution with an
exponential high-energy cutoff for the rigidity of ions. The electron
distribution, however, can be a broken power law with an exponential
($\delta=1$) or super exponential ($\delta=2.0$) cutoff. Since the
acceleration of low energy particles are not affected by radiative energy
loss, we assume that the electrons and ions have the same spectral index at
low energies.

The upper-right panel of Figure 13 shows a leptonic scenario for the
$\gamma$-ray emission. The corresponding model parameters are given in the
fourth row for RX J1713.7$-$3946 in Table 2. Since we adopt a much higher
energy density of 1 eV cm-3 for the infrared background photons and ion
processes also have significantly contribution to the $\gamma$-ray emission
(solid line), a high value of 22 $\mu$G is inferred for the magnetic field.
However, the total energy of the magnetic field is on the order of $10^{48}$
ergs, which is still much lower than that for ions. Since the ion spectral
cutoff is not well-constrained by the data, we set it at $70$ TeV as we did
above for other sources. An increase of the magnetic field will lead to a
shift to the dominance of $\gamma$-ray fluxes by the hadronic processes. The
upper left panel of Figure 13 shows our favored hadronic model for the
$\gamma$-ray emission. The corresponding model parameters are given in the
third row for RX J1713.7$-$3946 in Table 2. The cutoff energy of protons is 40
TeV, which is lower than that inferred from $\gamma$-ray observations (H. E.
S. S. Collaboration et al., 2018a) for the adoption of a single power-law ion
distribution here. The spectral index of 1.6 is comparable to the low energy
spectral index inferred from $\gamma$-ray observations and the product of
$W_{\rm p}$ and $n_{\rm H}$ is also compatible with these observations. For
such a hard spectrum, a broken power law electron distribution is needed to
fit the radio to X-ray spectrum via the synchrotron process. For the magnetic
field of 62 $\mu$G, we find a break energy of $82$ GeV and a cutoff energy of
$38$ TeV. Although the electron cutoff energy is very close to the proton
cutoff energy, their distributions are quite different at high energies for
the differences in the spectral index and shape of the cutoff $\delta$. The
total energy of the magnetic field is comparable to that of ions and is more
than 2 orders of magnitude higher than that of electrons. Note that although
$K_{\rm ep}=0.005$, the total energy of protons is more than 3 orders of
magnitude higher than that of electron for the hard spectrum and much lower
break energy of the electron distribution than the cutoff energy of protons
(Butt et al., 2008). Further increase of the magnetic field will lead to a
higher value of $W_{B}/W_{\rm e}$ and a slight decrease in $E_{\rm br}^{\rm
e}$, giving rise to a slightly higher ratio of $W_{\rm p}/W_{\rm e}$.

Figure 14: The upper two rows are the same as Figure 13 but for RX
J0852$-$4622\. The bottom left panel is similar to the top left panel except
for a lower $W_{\rm p}$ and a stronger magnetic field. The bottom right panel
is for a relatively simple single power-law leptonic model. However, the
cutoff of the electron distribution is exponential in this case.

To reduce the number of model parameters, one may increase the magnetic field
and set the electron radiative energy loss timescale at the break energy to be
equal to the age of the SNR. This leads to the fit in the lower left panel of
Figure 13. The corresponding model parameters are given in the first row for
RX J1713.7$-$3946 in Table 2. The $\gamma$-ray emission is completely
dominated by the hadronic processes. However, the total energy of the magnetic
field $W_{B}$ is more than 5 orders of magnitude higher than that of electrons
and more than 2 orders of magnitude higher than the energy of ions $W_{\rm
p}$. One may adjust the change of the spectral index $\Delta\alpha$ from low
to high energies to reduce the magnetic field. This leads to the fit in the
lower right panel of Figure 13, which is not as good as the others, especially
in the X-ray band. The corresponding model parameters are given in the second
row for RX J1713.7$-$3946 in Table 2. The slight decrease in $\alpha$ is due
to contribution to GeV $\gamma$-ray via the leptonic processes.

Figure 15: Same as Figure 14 but for HESS J1731$-$347\.

The spectral fits to the SED of RX J0852$-$4622 shown in Figure 14 are very
similar to those for RX J1713.7$-$3946\. The upper left panel shows our
favored model with the model parameters given in the fourth row for RX
J0852$-$4622 in Table 2. The spectral index of 1.75 for the ions is slightly
larger than 1.6 for the favored model of RX J1713.7$-$3946, and the cutoff
energy of the proton distribution is also slightly higher. With a magnetic
field of $31\ \mu$G, the magnetic field energy is about one tenth of the
proton energy and about 100 times higher than the energy of electrons. The
upper right panel shows a leptonic model with the model parameters given in
the fifth row for RX J0852$-$4622 in Table 2. The model slight overproduces
$\gamma$-rays at tens of TeV. Compared with the leptonic model for RX
J1713.7$-$3946, its break energy is about 100 times smaller, leading to a 20
times higher total energy of electrons. The total energy of protons is about
$3.5\times 10^{51}$ erg for the very hard spectrum and very high cutoff
energy, which is a bit too high to be reasonable. The total energy of the
magnetic field however is about $10^{49}$ erg and comparable to that of the
electrons.

The model of the middle right panel of Figure 14 is similar to that of the
lower right panel of Figure 13 for RX J1713.7$-$3946\. The model parameters
are given in the third row for RX J0852$-$4622 in Table 2. The model of the
middle left panel of Figure 14 is similar to that of the lower right panel of
Figure 13 for RX J1713.7$-$3946 with the model parameters given in the second
row for RX J0852$-$4622 in Table 2. To reduce the ratio of $W_{B}/W_{\rm e}$,
$\Delta\alpha=1.4$ instead of 1 is adopted for RX J1713.7$-$3946.

The model of the lower left panel of Figure 14 is similar to our favored model
in the upper left with the model parameter given in the first row for RX
J0852$-$4622 in Table 2. Here the total energy of proton $W_{\rm p}$ is 10
times smaller, leading to a very strong magnetic field of 165 $\mu$G and a
very high value of $3.9\times 10^{4}$ for $W_{B}/W_{\rm e}$. The model of the
lower right panel of Figure 14 has a single power-law electron distribution
with the model parameter given in the sixth row for RX J0852$-$4622 in Table
2. However to fit the SED, the shape of the cutoff needs to be exponential
instead of super exponential for other models. The energy of the magnetic
field is comparable to that of electrons, reminiscence of the leptonic model
with a broken power law electron distribution (upper right panel of Figure
14).

Figure 16: The spectral fits for SN 1006 and RCW 86. The left panels
correspond to the favored leptonic scenarios for the $\gamma$-ray emission.

The six spectral fits for HESS J1731$-$347 shown in Figure 15 are very similar
to those in Figure 14. The model parameters are given in Table 2. Due to its
relatively higher X-ray to $\gamma$-ray flux ratio, the magnetic fields in the
leptonic scenarios are higher than those for RX J0852$-$4622\. For the
leptonic models with a broken power-law electron distribution, the break
energy for HESS J1731$-$347 is between those for the other two sources, so is
the total energy of protons. For the hadronic models, we first notice that the
cutoff energy of protons is the lowest among these three sources and its
particle distributions are also harder than the other two sources. The
$\gamma$-ray luminosity of these three sources are comparable as can be seen
from the product of $W_{\rm p}$ and $n_{\rm H}$.

Compared with G296.5$+$10.0, SN 1006 and RCW 86 have prominent non-thermal
X-ray emission and are relatively younger. They all have relatively strong
radio emission. Figure 16 shows the spectral fits to these two SNRs with the
model parameters given in Table 2. Although both the leptonic (left panels)
and hadronic (right panels) models give reasonable fits to the SEDs, the
hadronic models are disfavored for their relatively stronger magnetic fields,
especially for RCW 86, whose hadronic model requires a magnetic field of more
than 10 mG. This is because the hard GeV spectrum is not compatible to the
soft radio spectrum if we assume electrons and ions have the same spectral
index. To fit the spectrum, the synchrotron spectrum needs to have a spectrum
break below the radio band, implying very strong magnetic field if this break
is associated with radiative energy loss processes. However, we notice that
multi-wavelength images of RCW 86 reveal complicated structure (Ajello et al.,
2016). Multi-zone hadronic models may still work. More detailed studies are
warranted. For SN 1006, $W_{p}$ is on the order of $10^{49}$ erg. A much
higher value can be ruled out without reducing $K_{\rm ep}$ since electrons
already have significant contributions to the $\gamma$-rays via the IC
emission process in both scenarios. We also notice that both leptonic models
for these two SNRs require an exponential cutoff instead of the super
exponential one and SNR 1006 has a single power law electron distribution
while RCW 86 has a broken one. The leptonic models for RWC 86 explored by
Ajello et al. (2016) only consider the CMB for the IC process and have a
single power law electron distribution, which is different from our favored
model.

Figure 17: Same as Figure 4 but for G78.2+01.2 (upper) and N132D (lower).

Besides the three SNRs studied above, there are two more radio bright SNRs:
G78.2+01.2 and N132D. Both of them have very soft spectra in the TeV band,
indicating a spectral cutoff. The latter is an SNR in the Large Magellanic
Cloud and is the most powerful SNR in our sample (Bamba et al., 2018).
Although thermal X-ray emission has been detected from both sources, there is
no evidence for non-thermal X-rays. A simple single power-law model can
readily fit their SEDs. The model parameters are given in Table 2. For
G78.2+01.2, we favor the model with $W_{\rm p}=10^{50}$ erg for the more
reasonable value of magnetic field and $W_{B}/W_{\rm e}$. Stronger magnetic
fields are needed for lower values of $W_{\rm p}$. Considering the fact that
N132D is the most powerful SNR, a value of $5\times 10^{50}$ erg appears to be
reasonable for $W_{\rm p}$. To reduce $W_{\rm p}$, the magnetic field needs to
be increased to reproduce the observed radio flux. The strong magnetic field
of $423\mu$ G appears to be reasonable to this SNR. So both models for N132D
shown in Figure 17 are favored. N132D has a relatively soft $\gamma$-ray
spectrum. The proton distribution cuts off at 70 TeV. The cutoff energy of the
proton distribution for G78.2+01.2, however, is only about 10 TeV, the lowest
in our sample. It is likely that protons above 10 TeV have already escaped
from the SNR and may illuminate surrounding molecular clouds. More
observations in the TeV band are warranted.

##  5 Discussions

In general, all the models presented above give good fits to the SEDs. We
picked up the favored ones mostly based on the model parameters. Firstly, we
favor models with a total energy of protons below $10^{50}$ erg and a magnetic
field between 10 to 100 $\mu$G or as close as to this range as possible. Since
these SNRs are expected to dominate the flux of TeV cosmic rays, on average
each SNR should inject less than $10^{50}$ erg of energy to the cosmic rays
(Zhang et al., 2017). Secondly, the total energy of the magnetic field should
be comparable to that of the electrons or protons. In Table 2, we highlighted
these models with a bold face for the ratio of $W_{B}/W_{\rm e}$. It can be
seen that most of the favored models satisfies the first criteria except for
SNRs G78.2+01.2 and N132D, where the magnetic field is above 100 $\mu$G. A
magnetic field below $100\ \mu$G will require a total energy of protons
exceeding $10^{50}$ ergs. Most of the favored models have a $W_{\rm p}$ of
$10^{50}$ erg. Only for G150.3+04.5 and SN 1006, $W_{\rm p}$ is one order of
magnitude lower. This is reasonable since these SNRs are likely resulting from
Type Ia SNs for their lack of compact neutron stars within the remnants. N132D
is the most powerful SNR, a value of $W_{\rm p}$ exceeding $10^{50}$ erg is
acceptable.

Figure 18: Gas density $n_{\rm H}$ vs magnetic field $B$ for the favored
models for the 13 SNRs studied in this paper. The open signs below $n_{\rm
H}=0.21$cm-3 are for the leptonic models. The dashed lines satisfy the
inserted equation, which assumes synchrotron process for the radio flux
density $f_{\rm R}$ and hadronic processes for the $\gamma$-ray flux density
$f_{\rm GeV}$ and gives the model-predicted relationship between $n_{\rm H}$
and $B$ in the hadronic scenario for different SNRs. The dotted line indicates
the gas density that can give a $\gamma$-ray flux density at 1 GeV via
hadronic processes $F_{\rm pp}$ equal to the 1 GeV flux density produced via
IC scattering of CMB $F_{\rm IC}$ by energetic electrons with $K_{\rm
ep}=0.005$ and $\alpha=2.0$.

Figure 18 shows the scatter plot between $B$ vs $n_{\rm H}$ for the favored
models. The mean density of the emission region is always less than $10$ cm-3
except for N132D, which is consistent with the X-ray observations. It appears
that $n_{\rm H}=0.21$ cm-3 gives the dividing line between the leptonic and
hadronic models with the former having a lower density. Assuming a single
power-law distribution with $\alpha=2$ and $K_{\rm ep}=5\times 10^{-3}$, the
dotted line shows the density when the $\gamma$-ray flux at 1 GeV produced via
the hadronic processes equals to that produced via IC of the CMB. With the
decrease of $\alpha$, this line will shift toward high densities since the GeV
emission efficiency via the IC process increases. Along the dotted line, the
radio to $\gamma$-ray flux density ratio increases with the increase of $B$,
which is consistent with Figure 10 with radio brighter sources having stronger
magnetic fields. The radio flux density also depends on the spectral index
with weaker radio emission for sources with harder spectra. This explains the
relative strong magnetic field for the two weak radio sources HESS J1534-571
and G150.3+4.5. The exact strength of the magnetic field also depends on
contributions to the $\gamma$-ray via the leptonic processes in the hadronic
scenario. Below the dotted line, the leptonic process dominates the
$\gamma$-ray emission. The dashed lines indicate the correlation between
$n_{\rm H}$ and $B$ if the $\gamma$-ray emission is solely produced via the
hadronic scenario for 13 SNRs studied in the paper, where $f_{\rm R}$ and
$f_{\rm GeV}$ represent the flux density at 1 GHz and 1 GeV, respectively.
Equally good fits to the SEDs can be obtained along these lines. Note that via
the CO and/or HI observations, much higher densities are obtained for some
sources in our sample, e.g., $\sim 130$ cm-3 for RX J1713.7-3946 (Fukui et
al., 2012), $\sim 60$ cm-3 for HESS J1731-347 (Fukuda et al., 2014), $\sim
100$ cm-3 for RX J0852-4622 (Fukui et al., 2017), $\sim 75$ cm-3 for RCW 86
(Sano et al., 2019), $\sim 30-80$ cm-3 for N132D (Bamba et al., 2018; Sano et
al., 2020) and $\sim 45$ cm-3 for Gamma Cygni (MAGIC Collaboration et al.,
2020). A higher average density in the $\gamma$-ray emission region will
result in a larger magnetic field and a lower total energy of relativistic
protons. We obtained a lower limit to the total energy of cosmic rays of
$10^{48}-10^{49}$ ergs by adopting these densities. The number density of some
individual cloud can be as high as $10^{4}-10^{5}$ cm-3 in RX J1713.7-3946
(Sano et al., 2010; Maxted et al., 2012). Most of these high density regions
are associated with local clouds with a small volume filling factor, implying
that the shock acceleration mainly operates in a low density inter-cloud
medium as we assumed (see Table 2). Under the strong magnetic field of these
clouds, the broad-band SED can be well reproduced without considering the
higher value of $n_{\rm H}$ due to the energy depended-penetration of cosmic-
rays into the dense clouds. Considering that the distribution of high energy
particles should be more or less uniform due to diffusion, a good spatial
correspondence between TeV gamma-rays and interstellar neutral gases is
expected in the hadronic scenario for the $\gamma$-ray emission.

To explain the origin of hard $\gamma$-ray spectra from supernova remnants in
the hadronic scenario, Gabici & Aharonian (2014) first showed that the energy
dependent penetration of cosmic rays into dense emission regions can lead to a
very hard $\gamma$-ray spectrum at low energies. Detailed studies of shock
interaction with molecular clouds show that the magnetic field strength will
be enhanced not only on the surface of targeted clouds, but also inside these
clouds (Zirakashvili & Aharonian, 2010; Inoue et al., 2012; Inoue, 2019). The
model was further developed by Celli et al. (2019). Detailed modeling of RX
J1713.7-3946 also favors the hadronic scenario (Zhang & Chen, 2016). Our
results suggest that the hard spectra may be due to very efficient particle
acceleration in a low density environment (Zhang et al., 2017; Zhang & Liu,
2019c). A softer spectrum can be produced when shocks slow down dramatically
due to interaction with molecular clouds (Bell et al., 2013; Tang & Chevalier,
2014), which may explain the low energy spectral component seen in HESS
J1912+101 and G279.0+1.1.

Figure 19: The total energy content of electron above 1 GeV $W_{e}$ vs
magnetic field $B$ for the 13 SNRs studied in this paper.

Figure 19 shows the scatter plot between $B$ and $W_{\rm e}$. The total energy
of electrons $W_{\rm e}$ is about $10^{47}$ erg, with N132D having the highest
value of $1.8\times 10^{48}$ erg and G150.3+04.5 having the lowest value of
$2.0\times 10^{46}$ erg. These results look reasonable. Table 2 shows that the
cutoff energies of these electrons are always around 1 TeV for sources without
nonthermal X-ray emission, which may explain the spectral break near 1 TeV in
the cosmic ray electron spectrum (DAMPE Collaboration et al., 2017). For the
three sources RX J1713.7$-$3946, RX J0852$-$4622, and HESS J1731$-$347, both
the leptonic and hadronic scenarios can lead to reasonable model parameters.
The leptonic models always have weaker magnetic fields and more energy in
electrons than the corresponding hadronic models. However, Table 2 shows that
leptonic models always have less energetic protons even we fix $K_{\rm ep}$
for the difference in their spectral shape. The leptonic model always has a
softer spectrum than the hadronic one. Better spectral measurements in the
radio band may distinguish these models.

Figure 20: The ratio of the magnetic energy $W_{B}$ to the total energy
content of electron above 1 GeV $W_{\rm e}$ for the 13 SNRs studied in this
paper.

Figure 20 shows the dependence of $W_{B}/W_{\rm e}$ on the age of the SNR. In
the leptonic models, the magnetic field energy is comparable to the total
energy of electrons. In the hadronic scenario, the magnetic energy is
comparable to that of protons and $W_{B}/W_{\rm e}$ appears to increase with
the age of the SNR.

##  6 Conclusions

In this paper, we carried out detailed spectral modeling of 13 SNRs with hard
GeV spectra. We re-analyzed the Fermi data of HESS J1912+101, and found its
TeV emission can not be attributed to leptonic process for the old age of the
SNR inferred from molecular cloud observations and the spin-down age of the
associated pulsar. The same is for SNR G279.0+1.1. A detailed analysis of XMM-
Newton observations of G296.5+10.0 failed to uncover nonthermal emission,
which in combination with Fermi data analyses also favors the hadronic
scenario for the $\gamma$-ray.

Of the 13 sources studied here, only SN 1006 and RCW 86 favor the leptonic
scenario for their $\gamma$-ray emission. RX J1713.7$-$3946, RX J0852$-$4622,
and HESS J1731$-$347 can be explained with both leptonic and hadronic models.
In the leptonic models, the total energy of the magnetic field is comparable
to that of the electrons. In the hadronic models, the magnetic fields and
protons are close to energy equipartition. All these sources have prominent
nonthermal X-ray emission. For the other 8 sources without evident nonthermal
X-ray emission, the hadronic models with a single power-law particle
distribution are favored. And in the hadronic scenario, the magnetic field of
older remnants tends to contain more energy than relativistic particles, which
may be attributed to escape of high energy particles from SNRs (Profumo et
al., 2018; MAGIC Collaboration et al., 2020).

Although our results do not completely address the origin of hard $\gamma$-ray
spectra from SNRs, young remnants with prominent nonthermal X-ray emission
favors the leptonic scenario, while absence of nonthermal X-ray emission
strongly favors the hadronic scenario. For RX J1713.7$-$3946, RX J0852$-$4622,
and HESS J1731$-$347, both scenarios can fit the SEDs with reasonable
parameters. The leptonic model always predicts a softer radio spectrum than
the corresponding hadronic one, which may be tested with future observations.

The proton spectrum always cuts off below 70 TeV, implying that SNRs may not
be able to accelerate cosmic ray to PeV energies. Alternative PeV sources are
needed to explain CR observations. SNRs of Type Ia SNs produce much fewer TeV
ions than core collapse SNs. TeV cosmic ray fluxes are therefore likely
dominated by more powerful SNRs with compact neutron stars in the middle. The
total energy of relativistic protons is on the order of $10^{50}$ erg for each
core collapse SNR, which indicates very efficient ion acceleration and is
compatible to the relatively hard spectra inferred from observations (Zhang &
Liu, 2019c). G78.2+01.2 has the lowest cutoff energy indicating escape of
particles beyond 10 TeV from the SNR. Nearby molecular clouds may be
illuminated by these escaping particles and produce $\gamma$-ray in the TeV
band. HESS 1912+101 has the oldest age in our sample, yet the cutoff energy of
protons is relatively high, indicating structure of the magnetic field may
play an important role on the escape process. The strong linear polarization
of the radio emission indeed indicates presence of large scale regular
magnetic field, which may trap high-energy particles in this source
effectively. For sources with prominent nonthermal X-ray emission, the
electron distribution always cuts off above 7 TeV, while those without
nonthermal X-ray emission, their electron distribution always cuts off near 1
TeV, which may explain the cosmic ray electron spectrum in the TeV range
(DAMPE Collaboration et al., 2017). Further exploration of this issue is
warranted.

Figure 21: The correlation between the electron synchrotron cooling time at
the cutoff energy and the age of SNRs.

To study the acceleration of high-energy particles in SNRs, Zeng et al. (2019)
plotted the synchrotron energy loss time at the cutoff energy of electron
distribution vs the ages of a sample of SNRs. This figure is updated here as
shown in Figure 21, confirming the early discovery that high-energy electrons
are mostly accelerated in young SNRs and the radiative energy loss and escape
processes dominate in old SNRs.

We thank the anonymous referee for very helpful suggestions that help to
improve the manuscript significantly. This work is partially supported by
National Key R&D Program of China: 2018YFA0404203, NSFC grants: U1738122,
U1931204, 11761131007, 11573070, the Natural Science Foundation for Young
Scholars of Jiangsu Province, China (No. BK20191109), and by the International
Partnership Program of Chinese Academy of Sciences, grant No.
114332KYSB20170008.

## References

  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, Science, 327, 1103, doi: 10.1126/science.1182787
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33, doi: 10.3847/1538-4365/ab6bcb
  * Abeysekara et al. (2018) Abeysekara, A. U., Archer, A., Aune, T., et al. 2018, ApJ, 861, 134, doi: 10.3847/1538-4357/aac4a2
  * Abeysekara et al. (2020) Abeysekara, A. U., Archer, A., Benbow, W., et al. 2020, ApJ, 894, 51, doi: 10.3847/1538-4357/ab8310
  * Acero et al. (2010) Acero, F., Aharonian, F., Akhperjanian, A. G., et al. 2010, A&A, 516, A62, doi: 10.1051/0004-6361/200913916
  * Aguilar et al. (2015a) Aguilar, M., Aisa, D., Alpat, B., et al. 2015a, Phys. Rev. Lett., 114, 171103, doi: 10.1103/PhysRevLett.114.171103
  * Aguilar et al. (2015b) —. 2015b, Phys. Rev. Lett., 115, 211101, doi: 10.1103/PhysRevLett.115.211101
  * Aharonian et al. (2007) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007, ApJ, 661, 236, doi: 10.1086/512603
  * Aharonian et al. (2008) Aharonian, F., Akhperjanian, A. G., Barres de Almeida, U., et al. 2008, A&A, 484, 435, doi: 10.1051/0004-6361:20078715
  * Ajello et al. (2016) Ajello, M., Baldini, L., Barbiellini, G., et al. 2016, ApJ, 819, 98, doi: 10.3847/0004-637X/819/2/98
  * Araya (2013) Araya, M. 2013, MNRAS, 434, 2202, doi: 10.1093/mnras/stt1162
  * Araya (2017) —. 2017, ApJ, 843, 12, doi: 10.3847/1538-4357/aa7261
  * Araya (2020) —. 2020, MNRAS, 492, 5980, doi: 10.1093/mnras/staa244
  * Bamba et al. (2008) Bamba, A., Fukazawa, Y., Hiraga, J. S., et al. 2008, PASJ, 60, S153, doi: 10.1093/pasj/60.sp1.S153
  * Bamba et al. (2018) Bamba, A., Ohira, Y., Yamazaki, R., et al. 2018, ApJ, 854, 71, doi: 10.3847/1538-4357/aaa5a0
  * Bell et al. (2013) Bell, A. R., Schure, K. M., Reville, B., & Giacinti, G. 2013, MNRAS, 431, 415, doi: 10.1093/mnras/stt179
  * Bocchino et al. (2000) Bocchino, F., Vink, J., Favata, F., Maggio, A., & Sciortino, S. 2000, A&A, 360, 671.  https://arxiv.org/abs/astro-ph/0005369
  * Braun et al. (2019) Braun, C., Safi-Harb, S., & Fryer, C. L. 2019, MNRAS, 489, 4444, doi: 10.1093/mnras/stz2437
  * Butt et al. (2008) Butt, Y. M., Porter, T. A., Katz, B., & Waxman, E. 2008, MNRAS, 386, L20, doi: 10.1111/j.1745-3933.2008.00452.x
  * Celli et al. (2019) Celli, S., Morlino, G., Gabici, S., & Aharonian, F. A. 2019, MNRAS, 487, 3199, doi: 10.1093/mnras/stz1425
  * Chang et al. (2008) Chang, C., Konopelko, A., & Cui, W. 2008, ApJ, 682, 1177, doi: 10.1086/589225
  * Clark et al. (1975) Clark, D. H., Caswell, J. L., & Green, A. J. 1975, Australian Journal of Physics Astrophysical Supplement, 37, 1
  * Cohen (2016) Cohen, J. M. 2016, PhD thesis, University of Maryland, College Park
  * Condon et al. (2017) Condon, B., Lemoine-Goumard, M., Acero, F., & Katagiri, H. 2017, ApJ, 851, 100, doi: 10.3847/1538-4357/aa9be8
  * DAMPE Collaboration et al. (2017) DAMPE Collaboration, Ambrosi, G., An, Q., et al. 2017, Nature, 552, 63, doi: 10.1038/nature24475
  * DAMPE Collaboration et al. (2019) DAMPE Collaboration, An, Q., Asfandiyarov, R., et al. 2019, Science Advances, 5, doi: 10.1126/sciadv.aax3793
  * De Luca & Molendi (2004) De Luca, A., & Molendi, S. 2004, A&A, 419, 837, doi: 10.1051/0004-6361:20034421
  * Devin et al. (2020) Devin, J., Lemoine-Goumard, M., Grondin, M.-H., et al. 2020, arXiv e-prints, arXiv:2009.08397.  https://arxiv.org/abs/2009.08397
  * Dickel et al. (1996) Dickel, J. R., Green, A., Ye, T., & Milne, D. K. 1996, AJ, 111, 340, doi: 10.1086/117786
  * Dickel & Milne (1995) Dickel, J. R., & Milne, D. K. 1995, AJ, 109, 200, doi: 10.1086/117266
  * Doroshenko et al. (2017) Doroshenko, V., Pühlhofer, G., Bamba, A., et al. 2017, A&A, 608, A23, doi: 10.1051/0004-6361/201730983
  * Duncan & Green (2000) Duncan, A. R., & Green, D. A. 2000, A&A, 364, 732.  https://arxiv.org/abs/astro-ph/0009289
  * Duncan et al. (1995) Duncan, A. R., Haynes, R. F., Stewart, R. T., & Jones, K. L. 1995, MNRAS, 277, 319, doi: 10.1093/mnras/277.1.319
  * Dyer et al. (2009) Dyer, K. K., Cornwell, T. J., & Maddalena, R. J. 2009, AJ, 137, 2956, doi: 10.1088/0004-6256/137/2/2956
  * Fleischhack (2019) Fleischhack, H. 2019, in International Cosmic Ray Conference, Vol. 36, 36th International Cosmic Ray Conference (ICRC2019), 675.  https://arxiv.org/abs/1907.08572
  * Fukuda et al. (2014) Fukuda, T., Yoshiike, S., Sano, H., et al. 2014, ApJ, 788, 94, doi: 10.1088/0004-637X/788/1/94
  * Fukui et al. (2003) Fukui, Y., Moriguchi, Y., Tamura, K., et al. 2003, PASJ, 55, L61, doi: 10.1093/pasj/55.5.L61
  * Fukui et al. (2012) Fukui, Y., Sano, H., Sato, J., et al. 2012, ApJ, 746, 82, doi: 10.1088/0004-637X/746/1/82
  * Fukui et al. (2017) —. 2017, ApJ, 850, 71, doi: 10.3847/1538-4357/aa9219
  * Gabici & Aharonian (2014) Gabici, S., & Aharonian, F. A. 2014, MNRAS, 445, L70, doi: 10.1093/mnrasl/slu132
  * Gao et al. (2011) Gao, X. Y., Han, J. L., Reich, W., et al. 2011, A&A, 529, A159, doi: 10.1051/0004-6361/201016311
  * Gerbrandt et al. (2014) Gerbrandt, S., Foster, T. J., Kothes, R., Geisbüsch, J., & Tung, A. 2014, A&A, 566, A76, doi: 10.1051/0004-6361/201423679
  * Giacani et al. (2000) Giacani, E. B., Dubner, G. M., Green, A. J., Goss, W. M., & Gaensler, B. M. 2000, AJ, 119, 281, doi: 10.1086/301173
  * Giuliani et al. (2011) Giuliani, A., Cardillo, M., Tavani, M., et al. 2011, The Astrophysical Journal, 742, L30, doi: 10.1088/2041-8205/742/2/l30
  * Guo et al. (2018) Guo, X.-L., Xin, Y.-L., Liao, N.-H., et al. 2018, ApJ, 853, 2, doi: 10.3847/1538-4357/aaa3f8
  * H. E. S. S. Collaboration et al. (2011a) H. E. S. S. Collaboration, Abramowski, A., Acero, F., et al. 2011a, A&A, 531, A81, doi: 10.1051/0004-6361/201016425
  * H. E. S. S. Collaboration et al. (2011b) —. 2011b, A&A, 531, A81, doi: 10.1051/0004-6361/201016425
  * H. E. S. S. Collaboration et al. (2015) H. E. S. S. Collaboration, Abramowski, A., Aharonian, F., et al. 2015, Science, 347, 406, doi: 10.1126/science.1261313
  * H. E. S. S. Collaboration et al. (2018a) H. E. S. S. Collaboration, Abdalla, H., Abramowski, A., et al. 2018a, A&A, 612, A6, doi: 10.1051/0004-6361/201629790
  * H. E. S. S. Collaboration et al. (2018b) —. 2018b, A&A, 612, A8, doi: 10.1051/0004-6361/201730737
  * H. E. S. S. Collaboration et al. (2018c) —. 2018c, A&A, 612, A7, doi: 10.1051/0004-6361/201630002
  * H. E. S. S. Collaboration et al. (2018d) H. E. S. S. Collaboration, Abramowski, A., Aharonian, F., et al. 2018d, A&A, 612, A4, doi: 10.1051/0004-6361/201526545
  * Helder et al. (2013) Helder, E. A., Vink, J., Bamba, A., et al. 2013, MNRAS, 435, 910, doi: 10.1093/mnras/stt993
  * Helfand & Becker (1984) Helfand, D. J., & Becker, R. H. 1984, Nature, 307, 215, doi: 10.1038/307215a0
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116, doi: 10.1051/0004-6361/201629178
  * Higgs (1977) Higgs, L. A. 1977, AJ, 82, 329, doi: 10.1086/112054
  * Hughes et al. (1998) Hughes, J. P., Hayashi, I., & Koyama, K. 1998, ApJ, 505, 732, doi: 10.1086/306202
  * Inoue (2019) Inoue, T. 2019, ApJ, 872, 46, doi: 10.3847/1538-4357/aafb70
  * Inoue et al. (2012) Inoue, T., Yamazaki, R., Inutsuka, S.-i., & Fukui, Y. 2012, ApJ, 744, 71, doi: 10.1088/0004-637X/744/1/71
  * Jones (1968) Jones, F. C. 1968, Physical Review, 167, 1159, doi: 10.1103/PhysRev.167.1159
  * Katsuda et al. (2009) Katsuda, S., Petre, R., Long, K. S., et al. 2009, ApJ, 692, L105, doi: 10.1088/0004-637X/692/2/L105
  * Katsuda et al. (2008) Katsuda, S., Tsunemi, H., & Mori, K. 2008, ApJ, 678, L35, doi: 10.1086/588499
  * Kellett et al. (1987) Kellett, B. J., Branduardi-Raymont, G., Culhane, J. L., et al. 1987, MNRAS, 225, 199, doi: 10.1093/mnras/225.2.199
  * Kothes et al. (2006) Kothes, R., Fedotov, K., Foster, T. J., & Uyanıker, B. 2006, A&A, 457, 1081, doi: 10.1051/0004-6361:20065062
  * Kuntz & Snowden (2008) Kuntz, K. D., & Snowden, S. L. 2008, A&A, 478, 575, doi: 10.1051/0004-6361:20077912
  * Lagage & Cesarsky (1983) Lagage, P. O., & Cesarsky, C. J. 1983, A&A, 125, 249
  * Lazendic et al. (2004) Lazendic, J. S., Slane, P. O., Gaensler, B. M., et al. 2004, ApJ, 602, 271, doi: 10.1086/380956
  * Leahy et al. (2013) Leahy, D. A., Green, K., & Ranasinghe, S. 2013, MNRAS, 436, 968, doi: 10.1093/mnras/stt1596
  * Lemoine-Goumard et al. (2012) Lemoine-Goumard, M., Renaud, M., Vink, J., et al. 2012, A&A, 545, A28, doi: 10.1051/0004-6361/201219896
  * MAGIC Collaboration et al. (2020) MAGIC Collaboration, Acciari, V. A., Ansoldi, S., et al. 2020, arXiv e-prints, arXiv:2010.15854.  https://arxiv.org/abs/2010.15854
  * Maxted et al. (2012) Maxted, N. I., Rowell, G. P., Dawson, B. R., et al. 2012, MNRAS, 422, 2230, doi: 10.1111/j.1365-2966.2012.20766.x
  * Maxted et al. (2018) Maxted, N. I., Braiding, C., Wong, G. F., et al. 2018, MNRAS, 480, 134, doi: 10.1093/mnras/sty1797
  * Milne & Haynes (1994) Milne, D. K., & Haynes, R. F. 1994, MNRAS, 270, 106, doi: 10.1093/mnras/270.1.106
  * Mori (2009) Mori, M. 2009, Astroparticle Physics, 31, 341, doi: 10.1016/j.astropartphys.2009.03.004
  * Morris et al. (2002) Morris, D. J., Hobbs, G., Lyne, A. G., et al. 2002, MNRAS, 335, 275, doi: 10.1046/j.1365-8711.2002.05551.x
  * Porter et al. (2006) Porter, T. A., Moskalenko, I. V., & Strong, A. W. 2006, ApJ, 648, L29, doi: 10.1086/507770
  * Profumo et al. (2018) Profumo, S., Reynoso-Cordova, J., Kaaz, N., & Silverman, M. 2018, Phys. Rev. D, 97, 123008, doi: 10.1103/PhysRevD.97.123008
  * Qiao et al. (2019) Qiao, B.-Q., Liu, W., Guo, Y.-Q., & Yuan, Q. 2019, J. Cosmology Astropart. Phys, 2019, 007, doi: 10.1088/1475-7516/2019/12/007
  * Reich & Sun (2019) Reich, W., & Sun, X.-H. 2019, Research in Astronomy and Astrophysics, 19, 045, doi: 10.1088/1674-4527/19/3/45
  * Reynoso et al. (2004) Reynoso, E. M., Green, A. J., Johnston, S., et al. 2004, PASA, 21, 82, doi: 10.1071/AS03053
  * Sano et al. (2010) Sano, H., Sato, J., Horachi, H., et al. 2010, ApJ, 724, 59, doi: 10.1088/0004-637X/724/1/59
  * Sano et al. (2019) Sano, H., Rowell, G., Reynoso, E. M., et al. 2019, ApJ, 876, 37, doi: 10.3847/1538-4357/ab108f
  * Sano et al. (2020) Sano, H., Plucinsky, P. P., Bamba, A., et al. 2020, ApJ, 902, 53, doi: 10.3847/1538-4357/abb469
  * Smith et al. (2019) Smith, D. A., Bruel, P., Cognard, I., et al. 2019, ApJ, 871, 78, doi: 10.3847/1538-4357/aaf57d
  * Su et al. (2017) Su, Y., Zhou, X., Yang, J., et al. 2017, The Astrophysical Journal, 845, 48, doi: 10.3847/1538-4357/aa7f2a
  * Tanaka et al. (2008) Tanaka, T., Uchiyama, Y., Aharonian, F. A., et al. 2008, ApJ, 685, 988, doi: 10.1086/591020
  * Tanaka et al. (2011) Tanaka, T., Allafort, A., Ballet, J., et al. 2011, ApJ, 740, L51, doi: 10.1088/2041-8205/740/2/L51
  * Tang & Chevalier (2014) Tang, X., & Chevalier, R. A. 2014, ApJ, 784, L35, doi: 10.1088/2041-8205/784/2/L35
  * Tian et al. (2008) Tian, W. W., Leahy, D. A., Haverkorn, M., & Jiang, B. 2008, ApJ, 679, L85, doi: 10.1086/589506
  * Tsuji & Uchiyama (2016) Tsuji, N., & Uchiyama, Y. 2016, PASJ, 68, 108, doi: 10.1093/pasj/psw102
  * Vasisht et al. (1997) Vasisht, G., Kulkarni, S. R., Anderson, S. B., Hamilton, T. T., & Kawai, N. 1997, ApJ, 476, L43, doi: 10.1086/310493
  * Vogt & Dopita (2011) Vogt, F., & Dopita, M. A. 2011, Ap&SS, 331, 521, doi: 10.1007/s10509-010-0479-7
  * Wendker et al. (1991) Wendker, H. J., Higgs, L. A., & Landecker, T. L. 1991, A&A, 241, 551
  * Winkler et al. (2003) Winkler, P. F., Gupta, G., & Long, K. S. 2003, ApJ, 585, 324, doi: 10.1086/345985
  * Woermann & Jonas (1988) Woermann, B., & Jonas, J. L. 1988, MNRAS, 234, 971, doi: 10.1093/mnras/234.4.971
  * Wood et al. (2017) Wood, M., Caputo, R., Charles, E., et al. 2017, in International Cosmic Ray Conference, Vol. 301, 35th International Cosmic Ray Conference (ICRC2017), 824.  https://arxiv.org/abs/1707.09551
  * Xing et al. (2014) Xing, Y., Wang, Z., Zhang, X., & Chen, Y. 2014, ApJ, 781, 64, doi: 10.1088/0004-637X/781/2/64
  * Yuan et al. (2012) Yuan, Q., Liu, S., & Bi, X. 2012, ApJ, 761, 133, doi: 10.1088/0004-637X/761/2/133
  * Zavlin et al. (2000) Zavlin, V. E., Pavlov, G. G., Sanwal, D., & Trümper, J. 2000, ApJ, 540, L25, doi: 10.1086/312866
  * Zeng et al. (2019) Zeng, H., Xin, Y., & Liu, S. 2019, ApJ, 874, 50, doi: 10.3847/1538-4357/aaf392
  * Zhang et al. (2020) Zhang, H.-M., Xi, S.-Q., Liu, R.-Y., et al. 2020, ApJ, 889, 12, doi: 10.3847/1538-4357/ab5af6
  * Zhang & Chen (2016) Zhang, X., & Chen, Y. 2016, ApJ, 821, 43, doi: 10.3847/0004-637X/821/1/43
  * Zhang & Liu (2019a) Zhang, X., & Liu, S. 2019a, ApJ, 876, 24, doi: 10.3847/1538-4357/ab14df
  * Zhang & Liu (2019b) —. 2019b, ApJ, 874, 98, doi: 10.3847/1538-4357/ab09fe
  * Zhang et al. (1997) Zhang, X., Zheng, Y., Landecker, T. L., & Higgs, L. A. 1997, A&A, 324, 641
  * Zhang & Liu (2019c) Zhang, Y., & Liu, S. 2019c, MNRAS, 482, 5268, doi: 10.1093/mnras/sty3136
  * Zhang et al. (2017) Zhang, Y., Liu, S., & Yuan, Q. 2017, ApJ, 844, L3, doi: 10.3847/2041-8213/aa7de1
  * Zhang & Liu (2020) Zhang, Y.-R., & Liu, S.-M. 2020, Chinese Astron. Astrophys., 44, 1, doi: 10.1016/j.chinastron.2020.04.001
  * Zirakashvili & Aharonian (2010) Zirakashvili, V. N., & Aharonian, F. A. 2010, ApJ, 708, 965, doi: 10.1088/0004-637X/708/2/965