A revisited study of the unidentified -ray emission towards the SNR Kes 41
Key Words.:ISM: supernova remnants — ISM: individual objects: Kes 41 — gamma rays: ISM
Kes 41 is among the Galactic supernova remnants (SNRs) proposed to be physically linked to -ray emission at GeV energies. Although not conclusively, the nature of the -ray photons has been explained by means of hadronic collisions of particles accelerated at the SNR blast wave with target protons in an adjacent molecular clump. We performed an analysis of Fermi-Large Area Telescope (LAT) data of about 9 years to assess the origin of the -ray emission. To investigate this matter we also used spectral modelling constraints from the physical properties of the interstellar medium towards the -ray emitting region along with a revised radio continuum spectrum of Kes 41 (, ). We demonstrated that the -ray fluxes in the GeV range can be explained through bremsstrahlung emission from electrons interacting with the surrounding medium. We also consider a model in which the emission is produced by pion-decay after hadronic collisions, and confirmed that this mechanism cannot be excluded.
Molecular clouds (MCs) house stellar objects at different stages of their evolution, from star-forming regions (SFRs) to the remains of supernovae (SNe). The detection of rays at GeV and/or TeV energies from the molecular gas can serve to illuminate the high-energy particle production taking place in these objects located within the cloud and even in its vicinity. Indeed, at a first glance the spatial correspondence of a cloud with the -ray emission is susceptible to explain the latter via collisions of accelerated protons, for instance at supernova remnant (SNR) shock fronts, with target protons (and maybe heavier nuclei, Banik & Bhadra, 2017) in the ambient matter. However, the recognition of this mechanism is not always straightforward since leptonic processes can also result in -ray emission via either bremsstrahlung or inverse Compton scattering produced by relativistic electrons.
Certainly, the number of -ray emitting sources detected at GeV energies by the Fermi Gamma-ray Space Telescope is growing rapidly. More than 3000 sources were reported in the last catalogues presented by Acero et al. (2015) and Ackermann et al. (2016). However, only for 35 of these sources a reliable counterpart originated in the radio emission from SNRs was found. Remarkably, approximately only half of these SNRs are well found in interaction with surrounding molecular gas.111See http://www.physics.umanitoba.ca/snr/SNRcat/. As a counterpart to the GeV emission, several of the Fermi sources were also detected at TeV energies222See http://tevcat.uchicago.edu/., and/or in the radio/X-ray domains (see for instance, Castro et al., 2013; Acero et al., 2016). However, regardless of the existence of a spatially coincident counterpart, the comprehension of the relative contribution of hadronic and leptonic processes responsible for -ray production remains unclear for a large fraction of cases (e.g. Tanaka et al., 2011; Pivato et al., 2013; H. E. S. S. Collaboration et al., 2018a).
Here we deal with the nature of the GeV -ray emission identified towards (, )(,) with the Large Area Telescope (LAT) detector onboard the Fermi satellite. The source known as 3FGL J1838.64654 in the Fermi-LAT source catalogue (3FGL, Acero et al., 2015) lies in the GeV emitting region surveyed in this work.333The source 3FGL J1838.64654 is catalogued as FL8Y J1638.54654 in the preliminary version of the Fermi-LAT 8-year source list, which will be replaced by the forthcoming 4FGL catalogue. Given the spatial coincidence on the plane of sky, Liu et al. (2015) considered natural a link between the GeV emission and the Galactic SNR Kes 41. Moreover, using the properties of the molecular gas emission derived by Zhang et al. (2015) in a restricted 7 4 area along the western rim of Kes 41, Liu et al. (2015) proposed the hadronic interactions as the most feasible process to explain the production of the -ray emission in the GeV domain. The large-scale structure of the ambient matter in which Kes 41 evolves was recently investigated in the companion paper presented by Supan et al. (2018a). As a result, on the basis of molecular and atomic line emission data, respectively from The Mopra Galactic Plane CO Survey (Burton et al., 2013) and the Southern Galactic Plane Survey (SGPS, McClure-Griffiths et al., 2005), along with mid-infrared Spitzer (Churchwell et al., 2009; Carey et al., 2009) information, the authors reported the discovery of the natal cloud (26 in size, mass 10-30 M) of the SNR and several star-forming regions around it. Now, in the light of the newly-determined physical conditions of the -ray emitting cloud, together with an updated set of Fermi-LAT data and on radio observations, we re-analised the feasibility of hadronic and leptonic models to fit the broadband spectral energy distribution (SED) for the SNR/-ray source system.
2 Observations and data analysis
2.1 -ray observations
We analysed the field of Kes 41 at -ray energies by using GeV data acquired with the Fermi-LAT during 9 years of the mission since the beginning of it on August 08, 2008 to July 10, 2017 (that is, mission elapsed time between 239557417 and 521337605, respectively). This time interval greatly increases the observing time by 71% with respect to the analysis previously presented in Liu et al. (2015), by adding 3.3 years of observations with the Fermi-LAT. Our ojective is to obtain an up-to-date morphological and spectral picture of the GeV sky in the direction of the remnant Kes 41 (i.e., towards the source 3FGL J1838.64654), in order to investigate the nature of the observed -ray emission.
Events from a region of interest (ROI) with a radius of 10 centred in the source 3FGL J1838.64654 ( , ) were selected. The energy range 0.5-400 GeV was chosen for the spatial analysis, as a compromise between limited statistics and angular resolution. For the spatial analysis, the ROI was reconstructed with a pixel size of . Scientific products were obtained following a standard reduction chain, employing the Science Tools (ST) package version v10r0p5 and python, with the current version of event reconstruction (Pass 8, Atwood et al. 2013) and the latest instrument response functions (IRFs P8R2_SOURCE_V6).444Details about specifications of the mission, instrumental issues, current event reconstruction framework, background models, etc, can be found at the Fermi-LAT Science Support Center (FSSC) web page: https://fermi.gsfc.nasa.gov/ssc/. In order to consider “good” photons for our analysis, we only kept events filtered according good-time-intervals (GTIs) and source class events (evtype = 3), which are indicated for point-source analysis. Additionally, we discarded photons coming from a zenith angle greater than 90 to minimize contamination due to rays generated by interactions of cosmic rays with the upper atmosphere.
From the above mentioned filtering process we obtained a set of high-quality -ray data, which we analysed by means of the maximum likelihood technique (Mattox et al., 1996) implemented through the ST routine gtlike. We performed a binned likelihood analysis, dividing the 0.5-400 GeV range in 30 logarithmically-spaced energy bins. Background emission was modelled through the user-contributed script make3FGLxml.py including point as well as extended sources in the ROI from the 3FGL catalogue along with the gll_iem_v06 and iso_P8R2_SOURCE_V6_v06 models for the diffuse Galactic and isotropic extragalactic components, respectively. The likelihood optimization procedure was carried out fixing spectral parameters (i.e., spectral index, normalization, etc) of the sources beyond 5 from the ROI centre, in order to ensure convergence. A first approximation of the spectral parameters was obtained by using the Minuit optimizer, which was then refined with the NewMinuit optimizer until convergence. Normalization of the Galactic and extragalactic backgrounds were left free.
2.1.1 Spatial distribution
In order to reveal the appearance of the GeV source in the field, we obtained a test-statistics (TS) map of the region by comparing the emission from the source with that of the constructed background model by using the tool gttsmap. The TS parameter is defined as , where is the likelihood with source included and corresponds to the null hypothesis (i.e., the sky subtracting the GeV excess of the source). It is a useful parameter to evaluate the statistical significance of a -ray excess, as it is related to the detection significance of the source according to . The new TS map for the region in which we are interested is presented in Fig. 1, where the radio emission at 843 MHz obtained from the Sydney University Molonglo Sky Survey (SUMSS, Bock et al., 1999) is depicted using white contours. For ease of reference, the HII regions in the field are labelled in Fig. 1 with numbers from 1 to 8. The TS angular distribution, with an overall significance of 27, is consistent with those presented by Liu et al. (2015).
2.1.2 Spectral analysis
We performed a new binned likelihood analysis of the data in order to study the spectral energy distribution of the source. In this case we restricted the range to 0.3-143 GeV, which was divided into 6 spectral bands, and we carried out a separated likelihood analysis for each of them. The extreme energies and the corresponding (logarithmic) central energy , as well as the GeV fluxes for each interval are reported in Table 1, and represented in Figs. 3 and 4 (Sect. 4). The distribution of the spectral points reveals a tendency for the flux of the source to fall below detection limits at an energy 50 GeV, suggesting a possible low-energy cutoff. This spectral shape may support the idea that events confidently associated with the source for energies above this value are scarce. For energies 50 GeV only upper limits were obtained.
Assuming a power-law spectral shape, the resulting photon index is , which is consistent with the one obtained by Liu et al. (2015). The flux in the 0.3-143 GeV band is determined to be ph s. The detailed modelling of the broadband spectral energy distribution from radio up to GeV energies is presented in Sect. 4. The uncertainties reported in Table 1 include statistical as well as systematic errors. Systematic errors are mainly associated with the propagation of inaccuracies in the effective area of the instrument, the point spread function, and with uncertainties related to the normalization of the Galactic diffuse background. Uncertainties associated with the effective area of the Fermi-LAT vary across the energy range, reaching maximum values of 10% at the highest energy considered in our study (Ackermann et al. 2012, see also the FSSC web site555https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/Aeff_Systematics.html.). On the other hand, systematics associated with the PSF are of the order of 5% for energies below 100 GeV in our range, linearly increasing up to 20% at higher energies.666https://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html. To estimate the systematics associated with the Galactic background, we followed Abdo et al. (2009), that is, after obtaining the best-fit values from the binned likelihood procedure, we varied the normalization of this background in 6% in order to determine the departure of the spectral parameters under this artificially-fixed model from the previously determined best-fit values.
|Energy range||Flux ()||TS|
|(GeV)||(GeV)||( erg s)|
|0.5||0.30 0.84||13.8 2.3 4.4||133|
|1.4||0.84 2.34||11.5 1.1 2.8||163|
|3.9||2.34 6.55||8.4 0.8 2.1||150|
|10.9||6.55 18.3||2.9 0.7 1.3||28|
|30.6||18.3 51.2||$a$$a$footnotetext: 8|
|85.6||51.2 143||$a$$a$footnotetext: 3|
3 The supernova remnant Kes 41
Kes 41 (G337.80.1) is a Galactic SNR classified, on the basis of the thermal X-ray emission detected towards its interior, as a member of the thermal-composite SNR class (Zhang et al., 2015). Although it was proposed that the remnant could be the result of the collapse of a star no later than a B0-type, with a mass 18 (Zhang et al., 2015), the lack of a central compact object reported for this remnant brings into question the proposed nature of the progenitor star.
Here, we present a thorough look at Kes 41 which as displayed in Fig. 1 is seen superimposed in projection on the GeV emission detected by Fermi-LAT. A zoomed-in view of the remnant at 843 MHz from SUMSS is shown in Fig. 2a. The remnant exhibits a slightly elongated radio shell of about 6 in size (equivalent to an average size of 24 pc at the SNR distance of 12 kpc888The detection of the OH maser spot evidenced the encounter of the remnant with dense interstellar gas and served to place the remnant at a distance of 12 kpc (Koralesky et al., 1998).), brighter towards its Galactic western portion. In Supan et al. (2018a) it has been demonstrated that this bright emission is morphologically correlated with an enhancement in the molecular gas emission observed between 63 and 48 km s (see Fig. 3 in that paper).
To calculate the global radio spectral index of the remnant we measured the integrated flux density of the source on images taken from public surveys at 843 and 5000 MHz, and constructed the spectrum shown in Fig. 2b by adding these estimates to previously published flux density measurements for Kes 41. The list of the integrated flux density estimates is presented in Table 2. The reported values, provided that the information about the primary calibrators is available, were brought onto the flux density scale of Perley & Butler (2017). From the spectrum it is evident that the flux density value at 80 MHz lies below the general trend of the data, which in principle may indicate the presence of thermal absorption along the line of sight. As only one point is a scarce evidence to define a spectral turnover at low frequencies, we first excluded the flux density value at 80 MHz and calculated the integrated spectral index in the radio band by a single weighted power-law least-squares fit (where ) to the remaining data points in Fig. 2b. The derived value of the integrated spectral index is = 0.10, which is fairly consistent with the estimation presented by Whiteoak & Green (1996) based only on integrated fluxes measured at 408 and 843 MHz. The integrated spectral index is also similar to that measured in other SNRs without a compact remnant in their interior. On the other hand, if the low frequency turnover is considered valid, a weighted least-squares fit to all the integrated flux densities shown in Fig. 2b using a power law plus an exponential turnover at a fiducial frequency of 408 MHz (), produces a flux = 23.32.6 Jy and an optical depth =0.0370.012 indicative of a non significant absorption level at this frequency. The free-free optical depth at 80 MHz is =1.130.37 (). The global spectral index has the same value derived when absorption from thermal ionised gas along the line is not present. Certainly, more flux density measurements at frequencies below 100 MHz are needed to clarify the presence of the spectral turnover.
|80||Dulk & Slee (1972)|
|408||Shaver & Goss (1970)|
|408||Dulk & Slee (1972)|
|843||This work$a$$a$footnotetext: 843||18.0 ? 6.0||Whiteoak & Green (1996)|
|4850||Wright et al. (1994)|
|5000||Dulk & Slee (1972)|
4 Physical scenarios for the GeV rays
It is believed that the -ray radiation produced in supernova remnants originates by the interaction of accelerated hadrons and/or leptons with surrounding ambient matter and radiation. This high energy -ray flux is enhanced in presence of dense MCs (see Slane et al., 2015, for a review).
In Liu et al. (2015) the -ray flux of Kes 41 is modelled considering leptonic and hadronic scenarios. Assuming that the rays originate by inverse Compton scattering of accelerated electrons with the low energy photons of the Cosmic Microwave Background, they found that the required energy in accelerated electrons to reproduce the observed flux is of the order of erg. This result is about one order of magnitude larger than the canonical value of erg, which corresponds to a fraction of 10% of the typical energy released in a supernova explosion. They also developed hadronic models compatible with the observed flux and also with the canonical values of the energy that goes to the accelerated particles in a supernova explosion.
In the current work, we take the properties derived for the large cloud into account (Supan et al. 2018a) and use, for the first time, a broadband spectrum including observations at radio energies to review the plausibility of a hadronic and leptonic origin of the -ray flux in the region of the SNR Kes 41.
The energy distribution of the accelerated particles in both, the hadronic and leptonic, scenarios considered here is assumed to be a power law with an exponential cutoff,
where indicates the particle species (electrons or protons), is a normalization constant, is the spectral index, and is the cutoff energy.
Let us first consider a model in which leptons are the component responsible for the -ray flux. In this type of model the radio flux is due to synchrotron radiation emitted by the propagation of the accelerated electrons in the ambient magnetic field. On the other side, the -ray radiation is generated mainly by the interaction of the accelerated electrons with the ambient matter, via the nonthermal bremsstrahlung process, and with the ambient radiation field, via inverse Compton. We refer to Aharonian et al. (2010) and references therein for calculation details of the synchrotron emissivity. Regarding the inverse Compton and the nonthermal bremsstrahlung emission, we followed the method given in Jones (1968) for the former, while the method presented in Baring et al. (1999) was used to model the electron-electron bremsstrahlung interaction and in Koch & Motz (1959) and Sturner et al. (1997) for the electron-ion bremsstrahlung process.
In modelling the spectrum from radio to -ray energies, we consider the radio observations reported in Table 2 and the updated Fermi-LAT data presented in Sect. 2.1. The chi-square used to fit the data is given by,
where corresponds to the observed flux at energy with an uncertainty , and is the model under consideration with parameters . Here, is the ambient magnetic field intensity and the proton density. We notice that the last term in includes the uncertainty on the determination of the proton density. The values of and correspond to the E1 region enclosing the SNR Kes 41, which is drawn in Fig. 1. A comprehensive analysis of the physical conditions in this region is presented in Supan et al. (2018a). The inclusion of the mean proton density as one of the fitting parameters makes the fit to take into account its correlation with the other fitting parameters. Moreover, the covariance matrix of the fit depends on the estimated mean proton density. Therefore, the mean proton density uncertainty is properly propagated (including correlations) in subsequent calculations that are based on the fitting parameters like, for instance, the total energy in accelerated electrons corresponding to the model.
The spectral index of the accelerated electrons energy distribution is fixed during the minimization procedure. Figure 3 shows the fit of the experimental data for , this value is motivated by the prediction of the first order Fermi mechanism and it is contained in the one-sigma region of the radio data fit performed in Sect. 3. As Fig. 3 shows, the nonthermal bremsstrahlung process is the dominant contribution to the -ray part of the flux. The inverse Compton component is several orders of magnitude smaller than the corresponding one to nonthermal bremsstrahlung. The fitted parameters are , , and cm (the hat is included to emphasize the maximum likelihood estimates character of the parameters). The obtained cutoff energy of the accelerated electrons is quite small. This is due to the suppression present in the -ray data. It is worth mentioning that, the high value of the proton density requires a smaller number of accelerated electrons to reproduce the -ray data and then, a high value of the magnetic field is required to reproduce the radio data.
The energy in accelerated electrons is calculated by using Eq. (1). For the spectral index the following value is obtained,
The total energy in accelerated electrons is smaller than the canonical value of erg. The maximum energy in accelerated electrons is obtained by fixing the spectral index in , which is the upper limit of the one-sigma region obtained from the fit of the radio data (see Sect. 3). In this case, a good fit is still obtained with the following value for the total energy,
This value is still smaller than erg. Therefore, a leptonic origin of the -ray flux is compatible with the present data.
Let us now consider the hadronic model. In this type of models the rays originate in the decay of particles, mainly neutral pions, generated in the interactions of accelerated protons with ambient protons. As mentioned before, the energy spectrum of the accelerated protons is assumed to be a power law with an exponential cutoff (see Eq. (1)). The -ray flux at Earth can be written as
where is the distance to the source, is the differential cross-section in proton-proton collisions resulting in the -ray emission, and is the speed of light.
The differential cross section of rays originated in proton-proton interactions has been extensively studied in the literature. A comprehensive study have been performed in Kafexhiu et al. (2014), in which a parametrization of the -ray energy spectrum, originated in proton-proton collisions, has been obtained. The projectile proton energy range of the parametrization starts at the kinematic threshold reaching a maximum value of 1 PeV. The differential cross section at low proton energies ( GeV) is obtained from a compilation of experimental data, while at high proton energies it is derived from Monte Carlo simulations. It is very well known that the hadronic interactions at the highest energies are unknown. However, there are models that extrapolate low energy accelerator data to the highest energies. The parametrization in Kafexhiu et al. (2014) includes several options for the high energy part, in fact it has been performed for the hadronic interaction models implemented in Geant 4.10.0 (Agostinelli et al., 2003), PYTHIA 8.18 (Sjöstrand et al., 2008), Sibyll 2.1 (Fletcher et al., 1994), and QGSJET01 (Kalmykov et al., 1997). In this work, we used this parametrization of the proton-proton collisions resulting in the -ray emission with the PYTHIA 8.18 option for the high-energy part. It is worth mentioning that one of the most important aspects of the new parametrization is the detailed description of the -ray spectrum at low proton energies, which represents an important improvement over earlier approaches.
Figure 4 shows the fit of the -ray part of the SNR Kes 41 energy spectrum observed at Earth. In this case we considered the value of the proton density, cm, in the SNR surroundings, derived inside the E1 region (see also analysis in Supan et al. 2018a). In addition, for this part of our analysis the spectral index is fixed during the fitting procedure. As in the case of the leptonic model we consider the canonical value and the values in the one-sigma region corresponding to the fit of the radio data. The cutoff energy obtained for is .
The energy in accelerated protons calculated for the ambient proton density measured in the E1 region and for is given by,
which is smaller than the canonical value of erg. Moreover, considering the spectral index , the upper limit of the one sigma region for the fit of the radio data, the maximum value obtained for the energy in accelerated protons is
which, also in this case, is smaller than the canonical value.
It is noteworthy that from the leptonic and hadronic models developed in this section in which , a ratio between the differential number of accelerated electrons and the differential number of accelerated protons is obtained for eV. This ratio, defined as (Merten et al., 2017)
is crucial in modelling the nonthermal emission produced in cosmic ray sources. In our leptonic model it is assumed that the -ray part of the flux coming from the hadronic component is much smaller than the one corresponding to bremsstrahlung. Therefore, a value of is required in order to obtain a contribution from the proton component one order of magnitude smaller than the one corresponding to the bremsstrahlung process. This value is more than one order of magnitude larger than the one inferred from cosmic rays observations (Merten et al., 2017), which is generally assumed to be of the order of . However, we recall that such a number represents an average value and it is not possible to know the differential particle number ratio in each individual source that contributes to the observed average. The canonical value of can be obtained theoretically by assuming that the number of protons and electrons that are accelerated are the same and that the spectral indexes of both components after acceleration are also the same, taking a value of . As discussed in Merten et al. (2017), is affected by several factors which can strongly modify its canonical value. In particular, there are evidences, both from a theoretical and observational point of view, which show that the spectral indexes of both components can be different. In that work, it is also shown that even for smaller values of the difference between the two spectral indexes a large variation of is obtained. For instance, for the differential particle ratio varies in such a way that . Therefore, models which requires large values of even of order one cannot be discarded (see for instance H. E. S. S. Collaboration et al. 2018b). We cautiously note that the justification of a negligible contribution of the bremsstrahlung process in the leptonic model of Liu et al. (2015) is based on the assumption of the canonical value of .
The analysis using the Fermi-LAT data presented in this paper shows that the -ray part of the spectrum of Kes 41 can be dominated by either accelerated leptons producing bremsstrahlung emission or accelerated hadrons that generate -ray emission by interacting with the surrounding ambient matter. The -ray flux can be properly fitted by these two types of models and also the total energy in accelerated particles is, in all cases, smaller than the canonical 10% of the typical energy released in a supernova explosion.
5 Concluding remarks
This paper is the second in a series focused on the multiwavelength analysis of the counterparts to the unidentified -ray emission detected at GeV energies towards the SNR Kes 41. Motivated by the new results recently presented in Supan et al. (2018a) revealing the natal molecular cloud of Kes 41, we performed, for the first time, a global assessment of all available data from radio to rays in order to determine the relative contribution of the hadronic and leptonic processes to the overall spectrum. To do this, we re-analised the Fermi-LAT data to study the -ray radiation spatially coincident (in projection) with the recently unveiled large-scale molecular material in the region of the SNR Kes 41. Using the leptonic models analysed in this work, we demonstrated that if the main contribution to the -ray emission comes from the SNR then, both radio and -ray data can be successfully modelled by synchrotron radiation and bremsstrahlung mechanism, respectively. A leptonic contribution from the inverse Compton emission to explain the production of the emission at GeV energies can be easily excluded. We also modelled the Fermi-LAT spectral data by a hadronic scenario considering the molecular, atomic, and ionised interstellar gases and found that hadronic interactions in a region relatively close to Kes 41 provide a viable mechanism for explaining the observed -ray emission, which qualitatively agree with the results of Liu et al. (2015).
- Abdo et al. (2009) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, ApJ, 706, L1
- Acero et al. (2015) Acero, F., Ackermann, M., Ajello, M., et al. 2015, ApJS, 218, 23
- Acero et al. (2016) Acero, F., Ackermann, M., Ajello, M., et al. 2016, ApJS, 224, 8
- Ackermann et al. (2012) Ackermann, M., Ajello, M., Albert, A., et al. 2012, ApJS, 203, 4
- Ackermann et al. (2016) Ackermann, M., Ajello, M., Atwood, W. B., et al. 2016, ApJS, 222, 5
- Agostinelli et al. (2003) Agostinelli, S., Allison, J., Amako, K., et al. 2003, Nuclear Instruments and Methods in Physics Research A, 506, 250
- Aharonian et al. (2010) Aharonian, F. A., Kelner, S. R., & Prosekin, A. Y. 2010, Phys. Rev. D, 82, 043002
- Atwood et al. (2013) Atwood, W. B., Baldini, L., Bregeon, J., et al. 2013, ApJ, 774, 76
- Banik & Bhadra (2017) Banik, P. & Bhadra, A. 2017, Phys. Rev. D, 95, 123014
- Baring et al. (1999) Baring, M. G., Ellison, D. C., Reynolds, S. P., Grenier, I. A., & Goret, P. 1999, ApJ, 513, 311
- Bock et al. (1999) Bock, D. C.-J., Large, M. I., & Sadler, E. M. 1999, AJ, 117, 1578
- Burton et al. (2013) Burton, M. G., Braiding, C., Glueck, C., et al. 2013, PASA, 30, e044
- Carey et al. (2009) Carey, S. J., Noriega-Crespo, A., Mizuno, D. R., et al. 2009, PASP, 121, 76
- Castro et al. (2013) Castro, D., Slane, P., Carlton, A., & Figueroa-Feliciano, E. 2013, ApJ, 774, 36
- Churchwell et al. (2009) Churchwell, E., Babler, B. L., Meade, M. R., et al. 2009, PASP, 121, 213
- Dulk & Slee (1972) Dulk, P. & Slee, O. 1972, Australian Journal of Physics, 25
- Fletcher et al. (1994) Fletcher, R. S., Gaisser, T. K., Lipari, P., & Stanev, T. 1994, Phys. Rev. D, 50, 5710
- H. E. S. S. Collaboration et al. (2018a) H. E. S. S. Collaboration, Abdalla, H., Abramowski, A., et al. 2018a, A&A, 612, A6
- H. E. S. S. Collaboration et al. (2018b) H. E. S. S. Collaboration, Abdalla, H., Abramowski, A., et al. 2018b, A&A, 612, A5
- Haynes et al. (1978) Haynes, R. F., Caswell, J. L., & Simons, L. W. J. 1978, Australian Journal of Physics Astrophysical Supplement, 45, 1
- Jones & Dickey (2012) Jones, C. & Dickey, J. M. 2012, ApJ, 753, 62
- Jones (1968) Jones, F. C. 1968, Physical Review, 167, 1159
- Kafexhiu et al. (2014) Kafexhiu, E., Aharonian, F., Taylor, A. M., & Vila, G. S. 2014, Phys. Rev. D, 90, 123014
- Kalmykov et al. (1997) Kalmykov, N. N., Ostapchenko, S. S., & Pavlov, A. I. 1997, Nuclear Physics B Proceedings Supplements, 52, 17
- Koch & Motz (1959) Koch, H. W. & Motz, J. W. 1959, Reviews of Modern Physics, 31, 920
- Koralesky et al. (1998) Koralesky, B., Frail, D. A., Goss, W. M., Claussen, M. J., & Green, A. J. 1998, AJ, 116, 1323
- Liu et al. (2015) Liu, B., Chen, Y., Zhang, X., et al. 2015, ApJ, 809, 102
- Mattox et al. (1996) Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396
- McClure-Griffiths et al. (2005) McClure-Griffiths, N. M., Dickey, J. M., Gaensler, B. M., et al. 2005, ApJS, 158, 178
- Merten et al. (2017) Merten, L., Becker Tjus, J., Eichmann, B., & Dettmar, R.-J. 2017, Astroparticle Physics, 90, 75
- Perley & Butler (2017) Perley, R. A. & Butler, B. J. 2017, ApJS, 230, 7
- Pivato et al. (2013) Pivato, G., Hewitt, J. W., Tibaldo, L., et al. 2013, ApJ, 779, 179
- Shaver & Goss (1970) Shaver, P. A. & Goss, W. M. 1970, Australian Journal of Physics Astrophysical Supplement, 14, 77
- Sjöstrand et al. (2008) Sjöstrand, T., Mrenna, S., & Skands, P. 2008, Computer Physics Communications, 178, 852
- Slane et al. (2015) Slane, P., Bykov, A., Ellison, D. C., Dubner, G., & Castro, D. 2015, Space Sci. Rev., 188, 187
- Sturner et al. (1997) Sturner, S. J., Skibo, J. G., Dermer, C. D., & Mattox, J. R. 1997, ApJ, 490, 619
- Tanaka et al. (2011) Tanaka, T., Allafort, A., Ballet, J., et al. 2011, ApJ, 740, L51
- Whiteoak & Green (1996) Whiteoak, J. B. Z. & Green, A. J. 1996, A&AS, 118, 329
- Wright et al. (1994) Wright, A. E., Griffith, M. R., Burke, B. F., & Ekers, R. D. 1994, ApJS, 91, 111
- Zhang et al. (2015) Zhang, G.-Y., Chen, Y., Su, Y., et al. 2015, ApJ, 799, 103