Characterization of the GeV emission from the Kepler supernova remnant

The Kepler supernova remnant (SNR) is the only historic supernova remnant lacking a detection at GeV and TeV energies which probe particle acceleration. A recent analysis of Fermi-LAT data reported a likely GeV gamma-ray candidate in the direction of the SNR. Using approximately the same dataset but with an optimized analysis configuration, we confirm the gamma-ray candidate to a solid $>6\sigma$ detection and report a spectral index of $2.14 \pm 0.12_{\rm stat} \pm 0.15_{\rm syst}$ for an energy flux above 100 MeV of $(3.1 \pm 0.6_{\rm stat} \pm 0.3_{\rm syst}) \times 10^{-12}$ erg~cm$^{-2}$~s$^{-1}$. The gamma-ray excess is not significantly extended and is fully compatible with the radio, infrared or X-ray spatial distribution of the SNR. We successfully characterized this multi-wavelength emission with a model in which accelerated particles interact with the dense circumstellar material in the North-West portion of the SNR and radiate GeV gamma-rays through $\pi^{o}$ decay. The X-ray synchrotron and inverse-Compton (IC) emission mostly stem from the fast shocks in the southern regions with a magnetic field B$\sim$100 $\mu$G or higher. Depending on the exact magnetic field amplitude, the TeV emission could arise from either the South region (IC dominated) or the interaction region ($\pi^{o}$ decay dominated).


Introduction
The last Galactic supernova to be observed from Earth occurred on October 9, 1604 and a detailed report was produced by Johannes Kepler whose name is now attached to the supernova and its remnant. The Kepler SNR is most certainly the remnant of a Type Ia explosion but the large scale asymmetry with brighter emission towards the North from radio to X-rays (DeLaney et al. 2002;Cassam-Chenaï et al. 2004;Blair et al. 2007;Reynolds et al. 2007) has caused some confusion with a core collapse origin (see Vink 2017, for a review). This asymmetry is now thought to be associated with circumstellar medium (CSM) from a runaway supernova progenitor system with significant mass loss prior to the explosion in a single degenerate scenario (e.g. Bandiera 1987;Burkey et al. 2013;Katsuda et al. 2015).
Estimates for the distance to the SNR range widely, from 3 to 7 kpc in the literature (e.g. Reynoso & Goss 1999;Sankrit et al. 2005;Katsuda et al. 2008). The measurement of the proper motion of Balmer-dominated filaments using the Hubble space telescope at a 10-year interval combined with the independently derived shock velocity from spectroscopy (Hα line width) provides the most robust estimation at d = 5.1 +0.8 −0.7 kpc (Sankrit et al. 2016). Throughout the paper we will use a distance of 5 kpc and rescale the values from the literature (e.g. shock speed) to match this distance whenever possible.
In the X-ray band, the emission is dominated by thermal emission with strong lines and in particular Fe lines supporting a Type Ia origin (e.g. Cassam-Chenaï et al. 2004;Reynolds et al. 2007). Non-thermal emission from thin synchrotron-dominated filaments was later revealed by Chandra observations (Bamba et al. 2005;Reynolds et al. 2007). Proper motion studies of these synchrotron rims (Vink 2008;Katsuda et al. 2008) show fast shocks with velocities 1 ranging from ∼2000 km s −1 in the northern region to ∼5000 km s −1 in the South.
The slower velocities in the North are related to the higher CSM density in this direction. Measurement of the thickness of these filaments suggests a high magnetic field of 150-300 µG if the width is energy loss limited (Bamba et al. 2005;Parizot et al. 2006).
Despite being one of the youngest SNRs in our Galaxy with high velocity shocks and signs of dense material and interaction, Kepler was the only historic SNR without detected γ-ray emission until now. This has changed with the recent report by Xiang & Jiang (2021) of a ∼ 3.8 σ detection 2 with Fermi-LAT in the direction of the Kepler SNR. In addition, the recent detection of TeV γ-rays after a deep exposure with the H.E.S.S. telescopes (152 hours, Prokhorov et al. 2021, H.E.S.S. collaboration submitted) opens the window for a detailed γ-ray study of this young and historic SNR.
In this work we aim to transform the status of the Fermi-LAT discovery from likely candidate to solid detection by using a more sophisticated analysis with approximately the same dataset (see Sect. 2). In addition to the modest significance, Xiang & Jiang (2021) find a slightly offset best-fit position from the SNR. We thus analyze in detail if this offset is statistically compatible with the SNR morphology, as realized by multi-wavelength spatial templates. We conclude in Sect. 3 by modeling Kepler's multi-wavelength emission under the assumption that γ-rays are 1 Velocities were rescaled to a 5 kpc distance. 2 Significance associated to a Test Statistic of 22.94 with 4 degrees of freedom.
Article number, page 1 of 7 arXiv:2201.05567v1 [astro-ph.HE] 14 Jan 2022 A&A proofs: manuscript no. kepler-Fermi-AandA-accepted Fig. 1. Left panel: Zoom on the Fermi-LAT TS map at the Kepler SNR position above 1 GeV. The green contours are from the infrared 24 µm Spitzer map. The plus symbol and circles illustrate the best-fit position and the 68%/95% confidence contours. Middle panel: Fermi-LAT TS map of the 15 • × 15 • region of interest around the Kepler SNR above 100 MeV. Magenta plus symbols represent the sources from the 4FGL-DR2 catalog and the green symbols the added point sources. For both TS maps, the Kepler SNR was not included in the model. Right panel: Spectral energy distribution of the Kepler SNR obtained using the infrared spatial template. Blue error bars represent the statistical uncertainties, while the red ones correspond to the statistical and systematic uncertainties added in quadrature. emitted from the northern interacting region while synchrotron and inverse Compton emission arises mostly from the fast shocks in the southern region.

LAT data reduction and preparation
The Fermi-LAT is a γ-ray telescope which detects photons by conversion into electron-positron pairs in the range from 20 MeV to higher than 500 GeV (Atwood et al. 2009). The following analysis was performed using 12 years of Fermi-LAT data (2008 August 04 -2020 August 03). A maximum zenith angle of 90 • below 1 GeV and 105 • above 1 GeV was applied to reduce the contamination of the Earth limb, and the time intervals during which the satellite passed through the South Atlantic Anomaly were excluded. Our data were also filtered removing time intervals around solar flares and bright GRBs. The data reduction and exposure calculations were performed using the LAT f ermitools version 1.2.23 and f ermipy (Wood et al. 2017) version 0.19.0. We performed a binned likelihood analysis and accounted for the effect of energy dispersion (when the reconstructed energy differs from the true energy) by using edisp_bins = −3. This means that the energy dispersion correction operates on the spectra with three extra bins below and above the threshold of the analysis 3 . Our binned analysis is performed with 10 energy bins per decade, spatial bins of 0.02 • for the morphological analysis and 0.05 • for the spectral analysis over a region of 15 • × 15 • . We included all sources from the LAT 10-year Source Catalog (4FGL-DR2 4 ) up to a distance of 15 • from Kepler. Sources with a predicted number of counts below 1 and a significance to zero were removed from the model.
The summed likelihood method was used to simultaneously fit events with different angular reconstruction quality (PSF0 to PSF3 event types 5 ). The Galactic diffuse emission was modeled by the standard file gll_iem_v07.fits and the residual background and extragalactic radiation were described by a single isotropic component with the spectral shape in the tabulated model iso_P8R3_SOURCE_V3_v1.txt. The models are available from the Fermi Science Support Center (FSSC) 6 .
Since the point spread function (PSF) of the Fermi-LAT is energy dependent and broad at low energy, we started the morphological analysis at 1 GeV while the spectral analysis was made from 100 MeV up to 1 TeV.

Morphological analysis
The spectral parameters of the sources in the model were first fit simultaneously with the Galactic and isotropic diffuse emissions from 1 GeV to 1 TeV. During this procedure, a point source fixed at the position (RA J2000 , Dec J2000 = 262.62 • , −21.49 • ) reported by Xiang & Jiang (2021) was used to reproduce the γ-ray emission of the Kepler SNR. To search for additional sources in the region of interest (ROI), we computed a test statistic (TS) map that tests at each pixel the significance of a source with a generic E −2 spectrum against the null hypothesis: TS = 2(ln L 1 − ln L 0 ), where L 0 and L 1 are the likelihoods of the background (null hypothesis) and the hypothesis being tested (source plus background). We iteratively added four point sources in the model where the TS exceeded 25. We localized the four additional sources (RA J2000 , Dec J2000 = 255.56 • , −22.47 • ; 260.13 • , −26.95 • ; 263.04 • , −27.97 • , and 264.78 • , −16.30 • ) and we fit their power-law spectral parameters. We then localized the source associated with Kepler and we obtained our best Point Source model (PS) at RA J2000 = 262.618 • ± 0.023 • , Dec J2000 = -21.488 • ± 0.026 • . The radii at 68% confidence and 95% confidence provided by f ermipy are: 0.036 • and 0.058 • . During this fit, we left free the normalization of sources located closer than 4 • from the ROI center as well as the Galactic and isotropic diffuse emissions. Figure 1 (left panel) presents our best localization and confidence radii superimposed on the TS map above 1 GeV. We tested its extension by localizing a 2D symmetric Gaussian. The significance of the extension is calculated through TS ext = 2(ln L ext − ln L PS ) where L ext and L PS are the likelihood obtained with the extended and point-source model, respectively. For Kepler, we found TS ext = 0.3 indicating that the emission is not significantly extended compatible with the result from Xi-F. Acero et al.: Characterization of the GeV emission from the Kepler supernova remnant Table 1. Results of the fit of the LAT data between 1 GeV and 1 TeV using different spatial models. The second column reports the TS values obtained for each spatial model, while column 3 indicates the number of degrees of freedom k adjusted in the model.  Table 2. Impact on the source significance of different analysis setups above a 700 MeV energy threshold. For comparison our setup corresponds to configuration 1 and that of Xiang & Jiang (2021)  Finally, we also examined the correlation of the γ-ray emission from Kepler with multi-wavelength templates derived from the VLA at 1.4 GHz (DeLaney et al. 2002), Spitzer in infrared at 24 µm (Blair et al. 2007), and Chandra in the 0.5-4 keV energy band (Reynolds et al. 2007). We cannot use the likelihood ratio test to compare the hypothesis of a point source model to that of a multi-wavelength template, because the two models are not nested. However, we can use the Akaike information criterion (AIC; Akaike 1974) which takes into account the number k of independently adjusted parameters in a given model. The standard AIC formula is AIC = 2k − 2ln L. Here we estimate a ∆AIC = 2k − TS comparing the AIC of the null source hypothesis and the AIC of the source being tested.
The lowest ∆AIC value, reported in Table 1, is obtained for the infrared template from Spitzer. While a lower ∆AIC indicates statistical preference for a model, the similar ∆AIC values indicate that each of these models provides an equally good representation of the γ-ray signal detected by the LAT. This result is consistent with our expectations from the Fermi-LAT PSF (∼ 1 • at 1 GeV 7 ) and the 0.03 • SNR radius: all templates have a similar structure when convolved with the relatively broad PSF. Hence, in our analysis below, we adopt the infrared template but emphasise that any of these templates would yield similar residual TS maps or inferred spectral properties.

Spectral analysis
Using the best-fit infrared spatial template, we performed the spectral analysis from 100 MeV to 1 TeV. We first verified whether any additional sources were needed in the model by examining the TS maps above 100 MeV. Two additional sources were detected at RA J2000 , Dec J2000 = 264.19 • , −20.01 • ; 263.72 • , −21.36 • (not detected in the 1 GeV -1 TeV range used in Sect. 2.2). The TS map obtained with all the sources consid-ered in the model (see Fig. 1, middle panel) shows no significant residual emission, indicating that the ROI is adequately modeled. We then tested different spectral shapes for Kepler. During this procedure, the spectral parameters of sources located up to 4 • from the ROI center were left free during the fit, like those of the Galactic and isotropic diffuse emissions. We tested a simple power-law model, a logarithmic parabola, and a smooth broken power-law model. Again, the improvement between the power-law model and the two other models is tested using the likelihood ratio test. In our case, ∆T S is 2.4 for the logarithmic parabola model and 3.1 for the smooth broken power-law representation, indicating that no significant curvature is detected. Assuming a power-law representation, the best-fit model for the photon distribution yields a TS of 38.3 above 100 MeV, a spectral index of 2.14 ± 0.12 stat ± 0.15 syst , and a normalization of (2.71±0.57 stat ±0.26 syst ) ×10 −14 MeV −1 cm −2 s −1 at the pivot energy of 2947 MeV. This implies an energy flux above 100 MeV of (3.1 ± 0.6 stat ± 0.3 syst ) × 10 −12 erg cm −2 s −1 and a corresponding γ-ray luminosity of (0.93 ± 0.18 ± 0.09) × 10 34 erg s −1 at a distance of 5 kpc. The systematic errors on the spectral analysis depend on our uncertainties on the Galactic diffuse emission model, on the effective area, and on the spatial shape of the source. The first is calculated using eight alternative diffuse emission models following the same procedure as in the first Fermi-LAT supernova remnant catalog (Acero et al. 2016) and the second is obtained by applying two scaling functions on the effective area. We also considered the impact on the spectral parameters when changing the spatial model from the infrared template to the best point source hypothesis. These three sources of systematic uncertainties were added in quadrature. The Fermi-LAT spectral points shown in Fig. 1 (right panel) were obtained by dividing the 100 MeV -1 TeV energy range into 5 logarithmically-spaced energy bins and performing a maximum likelihood spectral analysis to estimate the photon flux in each interval, assuming a power-law shape with fixed photon index Γ=2 for the source. The normalizations of the diffuse Galactic and isotropic emission were left free in each energy bin as well as those of the sources within 4 • . A 95% C.L. upper limit is computed when the TS value is lower than 1. Spectral data points are given in Table 3. We examined the reason for the significant improvement of the derived TS value (of 38.3) with respect to the TS value of 22.9 above 700 MeV reported by Xiang & Jiang (2021). To do so, we re-analyzed the source above the same threshold 700 MeV using a point source localized at the position reported by Xiang & Jiang (2021). Their setup corresponds to configuration 4 in Table 2 and we find a TS value very similar to theirs (22.9). We tested several analyses with different spatial bin sizes, region sizes and with/without summed likelihood and the improvement step by step is shown in Table 2. The higher TS value that we find in our analysis is most likely due to the summed likelihood analysis and the finer spatial binning of 0.05 • (configuration 1). These improvements together with our lower energy threshold of 100 MeV boost our detection to the TS value of 38.3.

Discussion
We now model the γ-ray emission in a multi-wavelength context. As Kepler is a well studied SNR, we aim to build a coherent model by fixing as many parameters as possible from observations and theoretical grounds.
Article number, page 3 of 7 A&A proofs: manuscript no. kepler-Fermi-AandA-accepted Table 3. Fermi-LAT flux data points using the infrared template. The flux parameter with an asterisk denotes an upper-limit. The first (second) flux errors represent statistical (systematic) errors respectively.

Model motivation
Our assumption is that on the one hand, the observed GeV γ-ray emission is mostly of hadronic nature (π 0 decay) being radiated from the North-West hemisphere where the shock is in interaction with the dense CSM as traced by infrared and optical maps.
On the other hand, the leptonic components (synchrotron and inverse Compton) arise from high velocity regions mostly observed in the South with shock speed 8 ranging from 4000-7000 km s −1 (regions 4-12, Katsuda et al. 2008). X-ray synchrotron emission requires high-speed shock regions, but radio emission can be produced by slower shocks. However for simplicity we model the electron population with a single radio to X-ray population. Consequently we expect our model will underpredict the radio data points. We use the measured and inferred properties of the SNR to fix various parameters of our model. First we fix the fast shocks (leptonic components) at 5000 km s −1 and the slower shocks in the interacting region at 1700 km s −1 for a target density of 8 cm −3 as reported in Sankrit et al. (2016). Secondly we derive the electron spectral distribution assuming that the electron maximal energy is limited by synchrotron losses and that proton maximal energy is limited by the age of the remnant.

Theoretical context
Following the prescription of Parizot et al. (2006), the acceleration timescale for the particles to reach an energy E at the forward shock can be written as: where r is the compression ratio assumed at the shock, B 100 the downstream magnetic field in units of 100 µG, and V sh,3 the shock speed in units of 1000 km s −1 . The deviation to Bohm diffusion is parametrized by where D Bohm is the Bohm diffusion coefficient. At a value of one, the acceleration at the shock is the most efficient and the maximal reachable energy decreases for higher values of k 0 . In the loss limited regime, the electron maximal energy can be obtained by equating the acceleration timescale τ acc to the synchrotron loss time at the shock giving: wheref (r) ≡ f (r)/ f (4), with f (r) = √ r − 1/r.
The maximal proton energy is obtained by equating τ acc from Equation 1 with the age of the remnant of 400 years giving: where T 400 is the age of the remnant in units of 400 yrs. For the magnetic fields considered in this modeling (B ∼ 100 µG), the synchrotron cooling is non negligible and is modeled with a broken power-law with E break obtained by equating the age of the remnant and the synchrotron loss time downstream of the shock giving: Because of the cooling and the cut at the maximum energy, the electron population is modeled as an Exponentially Cutoff BrokenPowerLaw with a change of slope after E break to Γ 2 = Γ 1 + 1.
Assuming that the synchrotron emission is limited by cooling, the acceleration efficiency parameter k 0 can be indirectly estimated by comparing the shock speed and the cutoff energy of the X-ray synchrotron spectrum using Equation 34 of Zirakashvili & Aharonian (2007). Such a study was carried out by Tsuji et al. (2021) on a population of SNRs including Kepler. For Kepler's south-eastern regions where fast shocks are observed (Tsuji et al. 2021, reg 4-8 in Table 2), k 0 (their η) ranges from 2.0 to 3.2 which is equivalent to 3.1 to 5.0 9 when rescaled to a distance of 5 kpc instead of 4 kpc. Given the measured values mentioned above, we decided to fix k 0 = 3.4 (the median value of the distribution) for our modeling for both the electron and proton populations (radiative model shown in Fig. 2).
The radiative models from the naima packages (Zabalza 2015) have been used with the Pythia8 parametrization of Kafexhiu et al. (2014) for the π o decay. For the IC, a far infrared field (T=30 K, U ph = 1 eV cm −3 ) was used in addition to the CMB (Porter et al. 2006).

Multi-wavelength data and spectral energy distribution
For the multi-wavelength data presented in Fig. 2, we used the updated compilation of radio fluxes from Castelletti et al. (2021), the X-ray data from the Suzaku XIS + HXD instruments (covering the 3-10 keV and 15-30 keV band, Nagayoshi et al. 2021), and the H.E.S.S. flux points from Prokhorov et al. (2021, H.E.S.S. collaboration submitted). The newly derived Fermi-LAT flux points using the infrared spatial template are presented (same as Fig. 1).
The resulting adjusted models, shown in Fig. 2, are obtained with only four free parameters being the downstream magnetic field, a unique electron/proton spectral index, and the associated energy budgets (see Table 4). The electron population spectral index and the amplitude of the magnetic field are correlated in the spectral energy distribution (SED) fitting. Motivated by theoretical expectations from recent kinetic hybrid simulations (Diesing & Caprioli 2021) predicting spectral indices steeper than 2, we present two different scenarios for spectral indices of 2.3 and 2.2 corresponding to an intermediate (90 µG) and high magnetic field value (170 µG) , respectively. A viable scenario, not shown in Fig. 2, can also be obtained for B=250 µG if the spectral index is changed to 2.15. Such magnetic field value is compatible with the estimated values at the shock given the thin Table 4. List of parameters obtained from the modeling of the spectral energy distribution in different scenarios. Parameters in brackets are fixed from observables or theory (see Sect. 3.1 and Sect. 3.2) while other parameters are adjusted to the data. The energy budget values are integrated above 1 GeV.   Table 4. X-ray synchrotron filaments size (Bamba et al. 2005;Rettig & Pohl 2012). When changing k 0 to 1 for the protons (i.e. maximal acceleration efficiency), E p,max increases and slightly improves the model agreement with the H.E.S.S. flux points in the hadronic dominated case.
Assuming that the hadronic emission arises from a small angular region in the SNR could explain the modest energy budget required (5.6×10 48 erg). On the Spitzer infrared 24 µm and the Hubble Hα images from Sankrit et al. (2016), the North-West interacting region has an opening angle of ∼45 • . Assuming a similar angle in the third dimension, this spherical cap represents ∼15% of the SNR surface. The local proton energy budget is therefore equivalent to about 4% of the local kinetic energy assuming an energy explosion of 10 51 erg.
The intermediate and high magnetic field scenarios reproduce equally well the GeV to TeV flux points and cannot be disentangled from the SED analysis alone. However, we note that the inverse Compton emission dominates above 300 GeV in an intermediate magnetic field case while the hadronic emission dominates the entire γ-ray band for a high magnetic field scenario. Therefore if the IC emission arises from the fast moving shocks in the southern regions, the precise location of the TeV γ-ray emission might be able to constrain the hadronic or leptonic nature of the emission and indirectly the average magnetic field in the SNR. The distance between the dense interacting region in the North-West and the southern rim is of the order of 0.05 • . While this is at the limit of the H.E.S.S. telescopes source localization precision for a faint source, a comparison of the GeV and TeV best-fit positions could shed light on the nature of TeV γ-ray emission. With an increased sensitivity and spatial resolution, the next generation Cherenkov Telescope Array (Cherenkov Telescope Array Consortium et al. 2019) will locate with great accuracy the Kepler SNR γ-ray emission.

Kepler SNR γ-ray emission in context
The detection of the Kepler SNR at GeV and TeV energies, completes our high-energy view of historical SNRs. In this section we compare the global spectral properties of the young and likely hadronic-dominated SNRs Kepler, Tycho, and Cassiopeia A with respect to the older middle-aged SNRs W 44, IC 443, and Cygnus Loop. Note that other young SNRs such as SN 1006, RX J1713.7−3946 or RCW 86, showing a spectral slope Γ ∼ 1.5 at GeV energies (see e.g. Acero et al. 2015), are likely dominated by leptonic emission and are not considered in our sample. The distances to the sources are fixed to 3.33 ± 0.10 kpc for Cassiopeia A (Alarie et al. 2014), 4 ± 1 kpc for Tycho (Hayato et al. 2010), 5.1 +0.8 −0.7 kpc for Kepler (Sankrit et al. 2016), 735 ± 25 pc for Cygnus Loop (Fesen et al. 2018), 3.0 ± 0.3 kpc for W 44 (Ranasinghe & Leahy 2018), and 1.7 ± 0.1 kpc for IC 443 (Yu et al. 2019). Figure 3 compares the SED in terms of luminosity of the aforementioned SNRs. As a proxy to discuss spectral curvature, we estimated the hardness ratio HR=νF 1 TeV /νF 1 GeV . Tycho, Kepler, and Cassiopeia A exhibit a nearly flat spectrum (HR=0.2-0.4) while the curvature is stronger for IC 443 (HR=0.015) and W 44 (< 2 × 10 −3 ). Such a contrast is due to differences in the acceleration and emission mechanisms. In the young SNRs sample, high shock speeds (3000-6000 km s −1 ) are observed producing highly energetic CRs interacting with circumstellar material. The second sample exhibits lower shock speeds (few 100 km s −1 ) with the presence of radiative shocks where compression and re-acceleration of pre-existing CRs takes place producing a spectral break at lower energies than for the young SNR sample. We note that the separation in terms of luminosity in our sample is not as clear. While W 44 is 50-100 times more luminous at the shock in W 44 is interacting with dense molecular environments (few 100 cm −3 , Yoshiike et al. 2013) whereas it is closer to ∼1-10 cm −3 in Cygnus Loop (Fesen et al. 2018).

Conclusion
By using ∼12 years of Fermi-LAT data and a summed likelihood analysis with the PSF event types, we were able to confirm GeV γ-ray emission at a >6σ detection level that is spatially compatible with the Kepler SNR. From the analysis of this γ-ray emission we draw the following conclusions: • above 100 MeV, the source is detected with a TS=38.3 with a power-law index of 2.14 ± 0.12 stat ± 0.15 syst . • the source is not significantly extended with an upper limit on its extension of 0.09 • (SNR radius is 0.03 • ). • the SED is modeled in a scenario with only four free parameters (B, Γ e,p , W e , W p ), the rest being fixed from the literature and theoretical grounds. The GeV γ-ray emission is interpreted as π o decay from the North-West interaction region. The TeV emission could be IC dominated (B < 100 µG; expected peak location in the South) or π o decay dominated (B > 100 µG; expected peak location in the North-West). While this is at the limit of current generation instruments, a comparison of the Fermi-LAT and HESS best-fit positions and errors could help to disentangle the two scenarios. • assuming a particle density of 8 cm −3 , derived from infrared observations, and that the interaction region represents 15% of the SNR surface, the local fraction of kinetic energy transferred to accelerated particles is of the order 4%.