Gamma-ray emission from middle-aged supernova remnants interacting with molecular clouds: the challenge for current models Free

ABSTRACT

1 INTRODUCTION

In the past few years, both space-based GeV (10⁹ eV) observatories (Fermi and AGILE) and ground-based TeV (10¹² eV) observatories (H.E.S.S., MAGIC, and VERITAS) detect γ-ray emission from several middle-aged supernova remnants (SNRs), e.g. W44 (Abdo et al. 2010a; Giuliani et al. 2011), IC443 (Albert et al. 2007; Acciari et al. 2009; Abdo et al. 2010b), W28 (Aharonian et al. 2008a; Abdo et al. 2010c; Hanabata et al. 2014), and W51C (Abdo et al. 2009; Aleksić et al. 2012; Jogler & Funk 2016). Multiwavelength observations further reveal the spatial correlation between the γ-ray emission region and the molecular cloud (MC) interaction region with robust tracers like OH maser and/or powerful diagnostic like molecular line broadening (see e.g. Jiang et al. 2010 and Slane et al. 2015). The γ-ray emission is produced either by energetic electrons with Bremsstrahlung and Inverse Compton (IC) emission mechanism or by accelerated protons with π⁰-decay emission mechanism. Since dense MCs are ideal sites for proton–proton interaction, the observed spatial association with MCs implies that the γ-ray emission is likely to have a hadronic origin. The characteristic π⁰-decay signature around 67.5 MeVis considered to be a unique feature to distinguish hadronic emission from leptonic emission unambiguously. Recently, the π⁰-decay signature is proposed in SNRs W44, IC443 (Ackermann et al. 2013), and W51C (Jogler & Funk 2016), which is believed to be the first direct evidence for CR proton acceleration in SNRs, making old SNRs interacting with MCs (hereafter SNR/MC) an important class of objects in γ-ray sky.

Despite above exciting progress in observation, our theoretical understanding about CR acceleration and emission in old SNRs is still very limited. It is partly because the evolution of old SNRs is more complicated and is strongly affected by the surrounding interstellar medium. To date, two scenarios are developed to explain the MC association and the γ-ray emission with hadronic origin. One is the direct interaction scenario (Bykov et al. 2000; Uchiyama et al. 2010; Tang & Chevalier 2014, 2015; Lee et al. 2015; Cardillo, Amato & Blasi 2016), in which the remnant directly interacts with the MCs. The collision creates a cooling shock region with enhanced density and magnetic fields, where the accelerated protons and electrons are able to produce enhanced π⁰-decay emission in γ-ray and synchrotron emission in radio, respectively. The other one is the escaping scenario (Aharonian & Atoyan 1996; Fujita et al. 2009; Gabici, Aharonian & Casanova 2009; Li & Chen 2010; Ohira, Murase & Yamazaki 2011), in which the MCs interact with the CR particles escaping from an adjacent SNR passively. Due to the high-density and magnetic fields in the MCs, runaway CR protons and electrons are able to illuminate the clouds in γ-ray with π⁰-decay emission and in radio with synchrotron emission, respectively. In the following discussion, the non-thermal particles in the vicinity of an SNR are referred to as CR particles while the pre-existing CR background is instead referred to as ambient CRs.

The growing number of SNR/MC detected in γ-ray enables us to investigate their physical properties in general and put better constraint on the theoretical models. Previous studies about γ-ray emission in SNR/MC focus on the individual source. In this work, we compare the γ-ray spectra of 10 SNR/MC in the First Fermi SNR Catalog (Acero et al. 2016) to obtain deeper insight into the physical origin of the γ-ray emission. As a first attempt, in this study, we mainly focus on the shape of γ-ray spectra without providing detailed modelling for each individual SNR/MC.

In Section 2, we present the γ-ray spectra from a sample of 10 SNR/MC and then discuss the interesting features shown in the spectra. In Sections 3 and 4, we discuss the escaping scenario and the direct interaction scenario in detail, respectively. We start with a brief introduction to the scenario and then compare the model spectrum with observation to gain more insight into the origin of γ-ray emission. Section 5 is the discussion section.

2 COMPARISON OF γ-RAY SPECTRUM

Based on the spatial overlap of Fermi detection and the radio extension of known SNRs, Fermi collaboration classified 30 sources as likely GeV SNRs in their First Supernova Remnant Catalog (Acero et al. 2016). According to multiwavelength observation, 11 of the 30 sources are further identified as SNR/MC. In this paper, we focus on 10 of the 11 SNR/MC, as HB21 in the 1st Fermi SNR catalogue was not detected in the recent work by Ackermann et al. (2017).

2.1 γ-ray spectrum from 10 SNR/MC

In Fig. 1, we present the γ-ray flux (upper panel) and luminosity (middle panel) of 10 SNR/MC with data available in the literature. The distances of all SNR/MC are taken from table 6 in Acero et al. (2016). For W28, only the emission from the northern part of the remnant, which is spatially overlapped with HESS J1801–233, is plotted. In the lower panel, we scale the γ-ray flux around 1GeV to be |${\sim } 10^{-10}\rm \ erg \ s^{-1}\, cm^{-2}$| to compare the shape of all 10 SNR/MC spectra. Statistical error bars are also included in the lower panel to demonstrate the comparison more clearly.

Figure 1.

Upper panel: γ-ray flux E²dF/dE as a function of photon energy. Middle panel: γ-ray luminosity as a function of photon energy. Distance adopted for the luminosity calculation is indicated in the brackets in kpc. Lower panel: scaled γ-ray flux as a function of photon energy, where we normalize the flux around 1 GeV to be about 10⁻¹⁰ erg cm⁻² s⁻¹. Reference: W44 and IC443 (Ackermann et al. 2013), W51C (Jogler & Funk 2016), W28N (Abdo et al. 2010c), W49B (H.E.S.S. Collaboration et al. ), W30 (Ajello et al. 2012), CTB 37A (Aharonian et al. 2008b; Brandt & Fermi-LAT Collaboration 2013), W41 (Castro et al. 2013; H.E.S.S. Collaboration et al. 2015b), G349.7+0.2 (H.E.S.S. Collaboration et al. 2015a), and G357.7-0.1 (Castro et al. 2013).

Open in new tabDownload slide

According to Fig. 1, most of the γ-ray spectra peak at a few GeV. Below ∼1 GeV, the spectra are characterized by a rising feature, which is consistent with the π⁰-decay signature. Above ∼1 GeV, the spectra start to steepen and follow roughly a power-law profile with no clear sign for an exponential like cut-off. The GeV emission data from Fermi and AGILE appear to be smoothly connected with the TeV emission data from H.E.S.S, MAGIC, and VERITAS except W30 and W41, which implies that the GeV and TeV emissions are probably generated at the same region with the same emission mechanism.

In W30¹ (G8.7-0.1) and W41² (G23.3-0.3), the TeV emission is much harder than the GeV emission and is spatially overlapped with a PWN in the vicinity of SNR/MC (Ajello et al. 2012; H.E.S.S. Collaboration et al. 2015b). It is possible that the excess of TeV emission is due to a PWN in the line of sight. Multiwavelength observation of the two remnants is needed in future to distinguish the contribution from the PWN. In addition, the low-energy spectrum of W41 doesn’t show the rising feature, which might be because the low energy emission from W41 is not dominated by π⁰-decay.

G357.7-0.1 and W44 are not detected in the TeV band. The GeV flux of G357.7-0.1 is lower than the majority of SNR/MC. Hence, the TeV emission of G357.7-0.1 may have similar slope as the other SNR/MC, but is still below the detection limit. W44 is a well-known SNR/MC and is very bright in the GeV band. The non-detection of TeV emission in W44 indicates that its TeV spectrum is much softer than the other SNR/MC. In Section 4.6, we provide a possible explanation to the lack of TeV detection in W44.

The γ-ray luminosity of all 10 SNR/MC at 1 GeV varies from |$10^{34}\rm \ erg \ s^{-1}$| to |$10^{36}\rm \ erg \ s^{-1}$|⁠. The diversity in γ-ray luminosity is probably the result of different MC environments or/and different physical properties of SNRs.

It is found that five SNR/MC prefer a BPL proton spectrum, which are presented in Table 1. Moreover, all five SNR/MC have similar α₁ ∼ 2.3, α₂∼ 3.0, and p_br ∼ 150 GeV/c except W44, which exhibits a much steeper spectrum with no TeV detection. Recently, Neronov et al. (2017) found that γ-ray emission from nearby giant MCs in the Gould Belt also prefers a BPL proton spectrum with |$\alpha _1= 2.33^{+0.06}_{ -0.08}, \, \alpha _2=2.92^{+0.07}_{ -0.04}$|⁠, and |$p_{br}=18.35^{+6.48}_{ -3.57}$|⁠. It is interesting that the BPL proton spectrum derived from SNR/MC here has similar α₁ and α₂ but larger p_br compared to that obtained in the isolated giant MCs. Since the γ-ray emission from isolated giant MCs traces the ambient CRs, the trend discussed above is likely an indication for re-acceleration of pre-existing ambient CRs in SNRs and will be discussed in Section 4.4 with details.

The fitting results shown in Table 1 are derived with the naimapython package (Zabalza 2015) and only statistical errors are taken into account in the calculation. We got slightly different results for W44 and IC 443 comparing with Ackermann et al. (2013). It is possibly because the naima package applies the new parametrization of the π⁰-decay cross-sections in Kafexhiu et al. (2014), while Ackermann et al. (2013) adopt the parametrized cross-section in Kamae et al. (2006). Kafexhiu et al. (2014) pointed out that the parametrization provided in Kamae et al. (2006) show some unphysical features around the threshold energy. Besides, the formulae developed in Kamae et al. (2006) do not agree with the recent measurements of pp (proton-proton) inelastic cross-section by the TOTEM collaboration at the Large Hadron Collider (e.g. Antchev et al. 2013).

The rest five SNR/MC need updated data in future to test whether a BPL proton spectrum is preferred to explain the γ-ray emission. As G357.7-0.1 has only four data points in the GeV band while the TeV emission of W30 and W41 is likely contaminated by a PWN in the line of sight. CTB 37A and W28N lack data points below 1 GeV that are important to constrain the low-energy part of CR proton spectrum.

2.2 π⁰-decay signature

In the literature, the rising feature shown in the log(E²dF/dE) − log(E) plot around a few hundreds MeV is often referred to as the π⁰-decay signature and considered to be the unique feature to identify hadronic emission in SNRs. However, as we will show later, the physical interpretation of this rising feature is not trivial, which require further clarification as pointed out by Strong (2016).

In the proton–proton interaction, the dominant channel for γ-ray production is through the decay of secondary π⁰, i.e. proton–proton collision creates π⁰, which then quickly decays into two photons. If we assume isotropic decay of π⁰ in its rest frame, then the observed γ-ray emission is found to be symmetric about half of the π⁰ mass (67.5 MeV) in the log(dF/dE) − log(E) plot (Stecker 1971). More importantly, above symmetry in the π⁰-decay emission is independent of the primary proton spectrum, which becomes a unique feature to identify hadronic emission and is usually referred to as the π⁰-decay signature.

In the context of SNR/MC, the π⁰-decay emission is produced by the interaction between primary CR protons and thermal nuclei in the vicinity of SNRs. In the upper panel of Fig. 2, we show the π⁰-decay emission for a primary proton spectrum n(p) ∝ p^−2.4 (blue solid line) in log(dF/dE) − log(E) plot, where p is the proton momentum. The π⁰-decay signature is clearly illustrated in the plot as a symmetric bump-like feature around 67.5 MeV(dashed line). γ-ray data below 67.5 MeV are essential to reveal the symmetric π⁰-decay signature unambiguously.

Figure 2.

Upper panel: γ-ray flux dF/dE as a function of photon energy E_γ in arbitrary units. Lower panel: γ-ray flux E²dF/dE as a function of photon energy E_γ in arbitrary units. The colour points are the scaled γ-ray spectra taken from the lower panel of Fig. 1. The blue solid line indicates π⁰-decay emission for a primary proton spectrum n(p) ∝ p^−2.4, where p is the proton momentum. The black dashed line shows the photon energy of 67.5MeV.

Open in new tabDownload slide

In the widely used log(E²dF/dE) − log(E) plot, due to the multiplication of E² factor, the symmetric bump like feature changes into a rising feature with a flattening of spectral slope around 100 MeV (see the lower panel of Fig. 2). Since the unique symmetric feature of π⁰-decay is now hidden in the rising feature and reflected as a changing of spectral slope, the identification of π⁰-decay signature becomes much more difficult in the log(E²dF/dE) − log(E) plot. Here, we want to emphasize that it is more straightforward to search for π⁰-decay signature in the log(dF/dE) − log(E) plot than the commonly used log(E²dF/dE) − log(E) plot.

In Fig. 2, we also present the scaled γ-ray spectra from all 10 SNR/MC. In the upper panel, it is found that below ∼100MeV the scaled γ-ray spectra bend towards a harder spectrum with decreasing energy, which is consistent with the symmetric π⁰-decay signature. However, the lack of data points below 67.5 MeV prevents us from identifying the π⁰-decay signature unambiguously, as the above bending feature can also be explained by Bremsstrahlung emission with a BPL electron spectrum. In the lower panel, the scaled γ-ray spectra are consistent with a rising feature. But it is less clear why the γ-ray data below 67.5 MeV is crucial for identifying the π⁰-decay signature. This again demonstrates that π⁰-decay signature is better revealed in the log(dF/dE) − log(E) plot instead of log(E²dF/dE) − log(E) plot.

2.3 Hadronic origin of γ-ray emission

Next, we want to discuss how to identify hadronic emission in SNR/MC without data below ∼100MeV, as both Fermi and AGILE don’t have sensitivity below ∼100MeV. Recently, hadronic emission is identified in W44, IC443, and W51C (Giuliani et al. 2011; Ackermann et al. 2013; Jogler & Funk 2016) based on a combination of several factors. The first one is the detection of a rising feature in (E²dF/dE) − log(E) plot which is consistent with π⁰-decay. Secondly, the spatial correlation between the γ-ray emission region and the MC interaction region also favors hadronic origin of γ-ray emission (e.g. Jiang et al. 2010; Slane et al. 2015). In the end, detailed calculations show that both IC and Bremsstrahlung emission mechanism cannot reproduce the γ-ray emission naturally. With typical background photon field, the simulated IC emission is too low to explain the γ-ray data. Bremsstrahlung emission is able to reproduce the rising feature in observation but requires an internal break in the electron spectrum. However, there is no physical reason for a break around a few hundreds MeV in the electron spectrum. Besides, Bremsstrahlung emission cannot explain the TeV emission from SNR/MC. It is the combination of all these factors together which makes us believe that the γ-ray emission from W44, IC443, and W51C has a hadronic origin.

Based on the above discussion, we constrain our study to hadronic models in the rest of this paper. IC and Bremsstrahlung emission is assumed to be negligible for our discussion, which should be a good approximation for energy ≳ 1 GeV. With hadronic origin, the rising feature below ∼1 GeV in the observed γ-ray spectra of SNR/MC can be naturally explained by the π⁰-decay signature. Although the log(E²dF/dE) − log(E) plot cannot reveal the π⁰-decay signature very clearly, it provides more details about the spectral shape at high energy part. In the following sections, we mainly focus on the features in the high energy part of spectra above ∼1 GeV. As a result, all the figures will be in the log(dF/dE) − log(E) format.

In the next two sections, we describe the escaping scenario and direct interaction scenario in detail and then compare the model spectra with observation. The main difference between the two scenarios is the source of primary CR protons. The escaping scenario focuses on the CR particles that escaped from the remnant, while the direct interaction scenario instead investigates energetic particles confined within the remnant. In both scenarios, the CR particles are believed to be accelerated at the remnant shock through the diffusive shock acceleration (DSA) process (e.g. Bell 1978; Blandford & Eichler 1987).

In current non-linear theory of DSA, there are still two open questions, one is how do energetic particles manage to escape the shock region and the other is how are seed particle injected into the DSA process. Both problems are not fully understood at this point (e.g. Malkov & Drury 2001) and require a special prescription in the treatment of DSA. In young SNRs, free escape boundary is widely used to describe the particle escaping and thermal injection of seed particles is often assumed for particle injection. In middle-aged SNRs, both prescriptions, however, confront some challenges which will be discussed in the following sections.

In this paper, π⁰-decay emission from the proton–proton interaction is calculated with the parametrized γ-ray production cross-sections derived in Kafexhiu et al. (2014). The formula is found to be accurate within |$20\hbox{ per cent}$| accuracy from the kinematic threshold (280MeV) up to PeV energies. At low energy, the model is fitted with experimental data while at high energy it is tested with public available code results. Please see Appendix  A for more details.

3 ESCAPING SCENARIO

3.1 Basic idea

CR particles accelerated at SNR shock can escape from the remnant and then propagate into the surrounding interstellar medium (ISM) after gaining enough energy. When these escaping CR particles encounter a dense MC, they interact with thermal nuclei in the MCs and illuminate the clouds in the γ-ray sky through π⁰-decay emission (Aharonian & Atoyan 1996). Assuming that the remnant is a point source, Gabici et al. (2009) modelled the multiwavelength emission of an illuminated MC in detail. Later, Li & Chen (2010) and Ohira et al. (2011) extended the model to account for the finite size of an SNR.

3.2 Runaway CR spectrum and the π⁰-decay emission from illuminated MCs

In this section, we compare the γ-ray emission predicted from escaping model with that measured in observation. In the upper panel of Fig. 3, we present the runaway CR proton spectrum in an illuminated MC for different t_age and L₁, where t_age is the age of SNR and L₁ is the distance between the SNR and the nearby MCs. The runaway CR spectrum is calculated with the model in Ohira et al. (2011), which is described with details in Appendix  B. Different combination of t_age and L₁ corresponds to a different spatial configuration of SNR/MC, which are used to demonstrate the diversity of γ-ray spectra expected from the escaping scenario. In the lower panel of Fig. 3, we compare the model spectra with γ-ray data. Given the complexity and diversity of all SNR/MC, we focus on only the comparison of spectral shape, where the normalization factor is left as a free parameter. Both the model spectra and γ-ray data are scaled to have the same peak value.

Figure 3.

From upper to lower panel, it is the runaway CR proton spectrum E²dN/dE as a function of proton kinetic energy, the corresponding π⁰-decay emission E²dF/dE as a function of photon energy, and the scaled γ-ray data and model spectra as a function of photon energy. t_{age, 4} is the remnant age in 10⁴ yr and L₁ is the distance between the inner radius of MCs and the centre of the SNR in pc. See Appendix  B for more details about the calculation.

Open in new tabDownload slide

According to Fig. 3, the runaway CR proton spectrum is characterized by a sharp low-energy cut-off E_low, which is a natural result of the free escape boundary adopted in the escaping scenario. Under the assumption of free escape boundary, CR particles with higher energy escape from the remnant earlier during the SNR evolution and there is a threshold energy E_low for escaping CRs. Particles with energy E < E_low either haven’t escaped from the remnant yet or do not have enough time to diffuse into the nearby MCs, thus leaving a low-energy cut-off E_low in the runaway CR spectrum. The sharp cut-off feature is not shown in the corresponding π⁰-decay emission. It is mainly because photons emitted in π⁰-decay spread in a large energy range, which smooth the sharp feature. However, the peak of γ-ray flux in escaping model does depend on E_low and is shifted to higher energy with larger E_low. The trend is valid as long as the cut-off energy E_low is larger than the threshold energy 280MeV. In order to reproduce the observed γ-ray spectra, the low energy cut-off of runaway CR proton spectrum has to satisfy E_low ≲ 1 GeV as shown in the lower panel of Fig. 3.

3.3 Escape of GeV cr particles

E_low ≲ 1 GeV brings a fundamental problem to the escaping scenario, which is can/how does the low-energy CR particle (∼1 GeV) escape from middle-aged SNRs. Early work by Bell (1978) proposed that damping mechanism such as neutral-ion damping can suppress CR-induced turbulence and facilitate the escaping of CR particles. Recently, Ptuskin & Zirakashvili (2003) further investigate how the maximum momentum p_max of CR particles evolves in the Sedov–Taylor (ST) phase of an SNR. It is argued that both non-linear wave damping and linear neutral-ion damping are crucial for the escape of low-energy CR particles.

Before we continue our discussion, we want to clarify that the maximum momentum at the shock front p_max is equivalent to the escaping momentum of the remnant p_esc under the free escape boundary condition, i.e. particles with momentum p > p_esc = p_max are able to escape the remnant through the free escape boundary. The maximum momentum of particle within the SNRs p_SNR, however, can be larger than the maximum momentum at the shock front p_max, if there are energetic CR particles trapped in the interior of a remnant. In this paper, we assume p_SNR = p_max = p_esc for simplification.

In equations (3) and (4), from the first step to the second step we apply the same κ, a, and ξ_cr as that in Ptuskin & Zirakashvili (2003), see the detailed definition of all the parameters there. Smaller κ, a, and ξ_cr certainly can decrease p_max and make it easier for GeV particles to escape. But smaller κ, a, and ξ_cr also decrease p_max at early time evolution of an SNR and make it even more difficult for the remnant to become a PeV accelerator. Hence, if we assume SNRs are CR accelerators up to the knee energy ∼10¹⁵eV, then κ, a, and ξ_cr can’t be arbitrarily small and deviate too much from those adopted in Ptuskin & Zirakashvili (2003).

According to equations (3) and (4), the escape of GeV particles put strong constraint on the physical properties of ISM and SNR, e.g. the shock velocity u_sh has to satisfy u_sh ≲ 100 km s⁻¹. However, when the shock velocity slows down to u_sh ≲ 200 km s⁻¹ in typical ISM environment, radiative cooling becomes important. Because of cooling loss, the remnant gradually leaves the adiabatic ST phase and eventually enters the radiative phase, which further complicates the discussion. As a result, the picture proposed in Ptuskin & Zirakashvili (2003) for ST phase of SNR evolution may no longer be valid for those middle-aged SNRs. Besides, it is also unclear from observation whether strong wave damping is indeed happening in these middle-aged SNRs, which deserve further investigation.

3.4 Challenge for the model

The main challenge for escaping scenario is how to obtain a low-energy cut-off E_low ≲ 1 GeV. Ohira et al. (2011) propose that the observed SNRs and MCs are very close or even in physical contact. In this special configuration, CR particles with all energy are able to escape the remnant and then diffuse into the adjacent MCs. With above assumption, Ohira et al. (2011) are able to reproduce the γ-ray emission in several SNR/MC. However, when SNR and MCs are directly interacting with each other, the validity of escaping scenario becomes an open question. In the unshocked part of MCs, the γ-ray emission is possibly interpreted by the illuminated cloud model, while in the interaction region the γ-ray emission is most likely described by the direct interaction scenario.

In escaping scenario, the shape of model spectrum strongly depends on the spatial configuration of SNR/MC. Since the observed 10 SNR/MC are likely in different spatial configuration, hence we expect to detect a variety of γ-ray spectra with different shape, which however is not seen in current data as indicated by the lower panel of Fig. 3. The diversity of γ-ray spectra in SNR/MC expected from escaping scenario is also discussed in Gabici et al. (2009). Someone may argue that the 10 SNR/MC discussed here is biased to middle-aged SNRs in physical contact with MCs as proposed by Ohira et al. (2011). However, the illuminated MCs external to W28 and W44 also exhibit similar spectrum with strong GeV emission and peak at a few GeV. In addition, γ-ray spectrum from middle-aged SNRs with no signature of MC interaction, e.g. SNRs S147 (Katsuta et al. 2012) and Cygnus Loop (Katagiri et al. 2011), is also similar as those in SNR/MC, which cannot be explained by the escaping scenario.

If we attribute the γ-ray emission from all these objects to illuminated clouds, then it is quite puzzling why we didn’t see any objects with γ-ray spectrum like the blue and black lines in the lower panel of Fig. 3. γ-ray emission from many young SNRs does show strong TeV emission (see e.g. fig. 7 in Funk 2015). But the observed emission is likely dominated by leptonic emission instead of hadronic emission. In summary, the γ-ray spectra from 10 SNR/MC are inconsistent with the prediction from escaping scenario statistically, which need to be confirmed by future observation. This inconsistency not only challenges the escaping scenario but also challenges our understanding of CR escaping in SNRs. It might suggest that the free escape boundary widely adopted in the modelling of CR escaping is not a good assumption (Drury 2011). We have to keep in mind that the introduction of free escape boundary in non-linear DSA was originally intended to achieve a self-consistent treatment of particle acceleration at SNR shock front. It is not designed to investigate the distribution and propagation of escaping CR in the ISM surrounding an SNR.

3.5 Illuminated clouds in w28 and w44

The best candidates for illuminated MCs are γ-ray bright MCs which are adjacent to an SNR but are not in physical contact with the remnant. The γ-ray source HESS J1800-240 external to W28 (hereafter W28 240) and the γ-ray bright regions surrounding W44 (hereafter W44 MCs) are considered to be two good candidates in this category (Uchiyama et al. 2012; Hanabata et al. 2014).

The environment in W28 240 is very complicated and needs further explanation before we continue our discussion (Aharonian et al. 2008a; Hanabata et al. 2014). In the TeV band, W28 240 is resolved into three individual sources, 240A, 240B, and 240C. 240A is spatially coincident with two H ii regions, G6.1-0.6 and G6.225-0.569, while 240B is spatially associated with the ultra-compact H ii region W28A2. 240C is spatially overlapped with SNR G5.71-0.08, which is likely interacting with MCs according to the OH maser detection. Please see Hanabata et al. (2014) and references therein for more details.

In this section, we compare the γ-ray emission from MCs external to W28 and W44 with emission from isolated MCs in the Gould Belt (Neronov et al. 2017). In the upper panel of Fig. 4, we provide the γ-ray luminosity of W28 240, W44 MCs and several isolated MCs in the Gould Belt. The luminosity is normalized for a cloud mass of 10⁵ M_⊙ and is calculated with the distance and mass shown in Table 2. When both CO and dust estimated mass is available, the dust estimated cloud mass is applied in the calculation. It is found that the normalized γ-ray luminosity of 240B and W44 MCs is comparable with that from isolated MCs while the luminosities of 240A and 240C are higher than the isolated MCs. Since there are CO dark clouds (e.g. Grenier, Casandjian & Terrier 2005), CO emission as a tracer of MCs is likely to underestimate the cloud mass. If the cloud mass in W28 240 and W44 MCs is larger than the CO estimated value, then the corresponding normalized luminosity becomes even smaller and thus more consistent with the isolated MCs.

Figure 4.

Upper panel: γ-ray luminosity from W44 MCs, W28 240, and 7 isolated MCs in Gould Belt. The luminosity is normalized to a cloud mass of 10⁵ M_⊙ based on the mass data in Table 2. Lower Panel: scaled γ-ray luminosity to compare the spectral shape of different MCs. Reference: isolated MCs in Gould Belt (Neronov et al. 2017), W44 MCs (Uchiyama et al. 2012), and W28 240 (Hanabata et al. 2014). The distance applied in the luminosity calculation is also listed in Table 2.

Open in new tabDownload slide

In the lower panel of Fig. 4, we present the scaled γ-ray luminosity to compare the shape of different spectra. The spectral shape of W44 MCs and W28 240B qualitatively agree with that of isolated MCs. 240C and 240A show harder spectrum, which is attributed to the escaping CRs in the literature. However, there are other possible explanations which require further investigation. 240C is spatially overlapped with SNR G5.71-0.08, which is interacting with MCs. Hence, the γ-ray emission of 240C can be explained by direct interaction scenario without involving escaping CRs. The GeV data of 240A are harder than the TeV data from H.E.S.S with a break around 100 GeV. The excess of GeV emission in 240A is likely affected by the spatially coincident H ii regions. Based on the above discussion, the γ-ray emission in W44 MCs and W28 240 is possibly attributed to or at least partly attributed to the ambient CRs and other associated sources, which need to be tested by future observation.

4 DIRECT INTERACTION SCENARIO

4.1 Basic idea

When an SNR collides with nearby MCs, the shock front driven by the remnant is slowed down by dense clouds. As a result, the postshock temperature drops and triggers efficient radiative cooling. Because of cooling-induced compression, a thin shell with high-density and magnetic field eventually forms behind the shock front. The thin radiative shell is an ideal site for both π⁰-decay emission in γ-ray and synchrotron emission in radio (Blandford & Cowie 1982). The direct interaction scenario focuses on the interaction between SNRs and MCs, and is motivated by the slow MC shock detected in observation. Slow MC shocks are good indicators for direct MC interaction and have been identified in all 10 SNR/MC discussed here with robust tracers such as OH maser.

Uchiyama et al. (2010) investigate the interaction between a young (non-radiative) SNR and MCs, and then adopt the model to explain the γ-ray and radio emission in W44. Tang & Chevalier (2014, 2015) instead study the collision between an old (radiative) SNR and MCs. In the former case, the fast (non-radiative) shock in a young SNR is slowed down by dense MCs and then drives a radiative shell with enhanced density and magnetic field into the MCs. In this case, both the radio and γ-ray emissions are expected to follow the MC interaction region in morphology. In the latter case, the old (radiative) SNR has a radiative shell itself, which can produce a significant amount of radio and γ-ray emission. The collision between SNR and MCs then creates a region with even higher density and magnetic field, which is also able to generate enhanced radio and γ-ray emission. IC 443 is considered to be a good example for this case. Tang & Chevalier (2014, 2015) propose that such two components model can explain the discrepancy between radio and γ-ray morphology in SNR/MC like IC 443.

If the picture proposed by Tang & Chevalier (2014) is correct, then old (radiative) SNRs without MC interaction are also able to produce a significant amount of γ-ray emission in its cooling shell. The γ-ray luminosity of radiative SNRs without MC interaction is expected to be smaller than those SNR/MC, making it more difficult to detect such objects. S147 is probably a good example of this category with γ-ray peak luminosity of only |$10^{33}\rm \ erg \ s^{-1}$|⁠. Multiwavelength observations of S147 show that the γ-ray, radio, and H α emission regions are spatially correlated with each other very well (Xiao et al. 2008; Katsuta et al. 2012). Since γ-ray, radio, and H α emission are expected to trace the accelerated protons, electrons, and the cooling shell, respectively, the spatial correlation among them is consistent with emission from a radiative SNR. In Katsuta et al. (2012), the multiwavelength emission from S147 is studied with the escaping scenario, but require an extremely low-filling factor for the interaction region. Besides, there is no signature for MC interaction in S147. So it is more natural to attribute the γ-ray emission to the radiative shell of an old SNR without MC interaction. With accumulated Fermi data and CTA in future, we expect to detect more remnants like S147.

Next, we discuss the primary CR spectrum within the framework of direct interaction scenario. Several different ideas have been proposed to interpret the observed γ-ray emission with hadronic origin, including pure adiabatic compression of pre-existing ambient CRs, DSA of thermal injected seed particles, and re-acceleration of pre-existing ambient CRs.

4.2 Pure adiabatic compression

Assuming pure adiabatic compression, the diversity of SNR properties and ISM environment in SNR/MC is simply reflected in the compression ratio s. In Fig. 5, we plot the primary CR proton spectrum based on equation (6) (upper panel) and the corresponding π⁰-decay emission (middle panel) for different s. When s increases from 10 to 100, the compressed CR spectrum is boosted significantly in the vertical y-axis, while the shift of spectrum in the horizontal x-axis is relatively small. In the lower panel of Fig. 5, we present the scaled γ-ray data and model spectra which are normalized to have the same peak value. It is shown that the spectral shape of π⁰-decay emission based on pure adiabatic compression remains almost unchanged with different s. The compressed CR spectrum is capable to explain the spectral shape of SNR/MC like W44 and IC 443, which are consistent with the conclusion in Tang & Chevalier (2014) and Cardillo et al. (2016). But it has difficulty in reproducing the emission of SNR/MC with harder TeV spectrum. It might imply that the ambient CR spectrum in the vicinity of SNRs is harder than the local values measured in our solar neighbourhood, which will be discussed in Section 4.5.

Figure 5.

Same as Fig. 3, but the primary CR proton spectrum is taken from equation (6) due to pure adiabatic compression of ambient CRs. s is the adiabatic compression ratio.

Open in new tabDownload slide

The pure adiabatic compression case might correspond to the situation with a quasi-perpendicular shock, where DSA is inefficient. Hybrid simulations show that at very oblique shocks ions only gain a factor of a few in momentum and energy through shock drift acceleration (Caprioli & Spitkovsky 2014). Hence, in direct interaction scenario the γ-ray emission from SNR/MC like W44 is possibly explained by shock drift acceleration in quasi-perpendicular shock plus adiabatic compression. The main challenge for such an explanation is the origin of quasi-perpendicular shock. If the dense filaments in MCs have a preferred magnetic field direction along the filaments, then quasi-perpendicular shock might form during the interaction between SNR and the dense filaments.

4.3 DSA of thermal injected particles

In this section, we discuss about the situation involving DSA of thermal injected particles. In young SNRs, it is widely accepted that the seed particles injected into DSA process are those energetic particles in the high-energy tail of thermal population (e.g. Blasi 2013). In middle-aged SNRs, as a first guess it is natural to adopt the same assumption. However, as we will show later, thermal injected particles as seed particles have difficulties in reproducing the observed TeV emission.

Figure 6.

Same as Fig. 3, but the primary CR proton spectrum is described by equation (13) based on DSA of thermal injected seed particles. p_max and p_br are the cut-off momentum and break momentum in GeV/c, respectively. The adiabatic compression ratio s is assumed to be 50.

Open in new tabDownload slide

According to Fig. 6, the spectral break p_br induced by neutral-ion damping is crucial in interpreting the steepening of γ-ray spectrum above a few GeV. The main concern about the modified version of neutral-ion damping is the time-scale taken for it to rebuild the CR spectrum, which is not discussed in Malkov et al. (2011). As a result, whether such modified version of neutral-ion damping is efficient enough in middle-aged SNRs is still an open question (Drury 2011).

The primary CR spectrum derived in equation (13) is able to reproduce the γ-ray emission in W44 and G357.7-0.1, which do not have TeV detection (Uchiyama et al. 2010). For the rest of SNR/MC with TeV detection, the model discussed in this section however fails in the TeV band due to the exponential cut-off (Tang & Chevalier 2014). It is mainly because DSA is inefficient in middle-aged SNRs with slow shocks and is not able to accelerate thermal injected particles to above TeV energy, which suggests that thermal injected particles might not be the dominant component of seed particles in SNR/MC. In Section 4.2, we already show that pure adiabatic compression of ambient CRs is capable to produce a significant amount of TeV emission, which motivates us to consider ambient CRs as seed particles in DSA, i.e. re-acceleration of pre-existing ambient CRs.

4.4 Re-acceleration of pre-existing ambient CRs

In this section, we discuss the model involving re-acceleration of pre-existing ambient CRs (hereafter RPCR). From energetic point of view, in middle-aged SNRs with slow shocks RPCR is more efficient in accelerating particles. Recently, Lee et al. (2015) performed a hydro-simulation for W44 with a self-consistent treatment of non-linear DSA. The authors found that, in order to reproduce the observed γ-ray and radio emission, about |${\sim } 33\hbox{ per cent}$| of shock energy has to contribute to the DSA process for thermal injection, while only |${\lesssim } 1\hbox{ per cent}$| of acceleration efficiency is needed for RPCR.

If RPCR is the correct picture for DSA in middle-aged SNRs, then it indicates a transition of seed particles during the SNR evolution, which is from thermal injected seed particles in young SNRs to pre-existing CRs in middle-aged SNRs. The transition is a natural consequence of the declining shock velocity in the SNR evolution. For slow shock, the postshock temperature is low, thus the thermal particles behind the shock front are shifted to lower energy, which makes the thermal injection less efficient.

RPCR is potentially very interesting as it may help explain the anomaly detected by PAMELA and AMS in the ambient CR spectrum (e.g. Adriani et al. 2011). For example, Thoudam & Hörandel (2014) propose that strong re-acceleration of ambient CRs at slow SNR shocks is able to explain the observed spectral break in CR proton and Helium spectra at a few hundred GV rigidity. The study here provides possibly the first evidence for efficient RPCR in middle-aged SNRs.

Figure 7.

Same as Fig. 3, but the primary CR proton spectrum is replaced with equation (21) involving both RPCR and adiabatic compression. σ is the power-law index for energy-dependent diffusion. t/τ is a dimensionless time ratio, where t is roughly the remnant age and τ is a characteristic time for DSA, see equation (22). The adiabatic compression ratio s is assumed to be 50.

Open in new tabDownload slide

According to Fig. 7, as the ratio t/τ increases the critical momentum p_crit in the proton spectrum is shifted to higher energy, while the power-law index below and above the break remains almost the same; see more discussion in Tang & Chevalier (2015). This feature in RPCR can naturally explain why the BPL spectrum in SNR/MC has similar α₁ and α₂ but larger p_br comparing with that in isolated giant MCs as discussed in Section 2.1.

In the case of RPCR, the modified version of neutral-ion damping (Malkov et al. 2011) is not needed to explain the steepening of spectrum above a few GeV (Tang & Chevalier 2015; Cardillo et al. 2016). According to Fig. 7, the model involving both RPCR and adiabatic compression is able to reproduce the overall profile of γ-ray spectra in all SNR/MC except W30 and W41. The dispersion shown in the TeV spectra is likely partly explained by the different acceleration time, i.e. t/τ, and partly induced by the intrinsic dispersion in the ambient CR spectrum, which will be discussed in the following section. W30 and W41 show much harder TeV emission than the others, which is possibly caused by the spatially overlapped PWN in the line of sight.

4.5 Pre-existing ambient CR spectrum

For SNR/MC with hard TeV spectrum, the ambient CR spectrum measured in our solar neighbourhood seems to be too steep to explain the observed TeV emission. It might imply that the ambient CR spectrum in the vicinity of an SNR is harder than the local value. There are several possibilities which can induce a harder ambient CR spectrum.

At first, the ambient CRs surrounding an SNR is strongly affected by the escaping CR particles, which can potentially harden the ambient CR spectrum. As particles with higher energy escape from the remnant more easily. A self-consistent model involving both ambient CRs and escaping CR particles is needed in future to investigate this possibility. Next, the CR background in our Galaxy may not be uniform as assumed in the standard model of Galactic CRs. It is possibly linked to the very high energy γ-ray emission discovered by H.E.S.S. in the Galactic Centre ridge (Aharonian et al. 2006) and the excess of GeV emission discovered by Fermi in the inner Galaxy (Ackermann et al. 2012). One possible explanation for the excess of diffuse γ-ray emission is that the CR background in our Galaxy is not uniform and the ambient CR spectrum in the Galactic central region is harder than our local value. If the ambient CR spectrum is not uniform in our Galaxy, then the dispersion shown in TeV emission of all SNR/MC simply reflects the non-homogeneity of CR background in our Galaxy. Ambient CR spectrum with spectral index of 2.5−2.7 at very high energy would be enough to explain the TeV emission from different SNR/MC.

4.6 Challenge for the model

In the direct interaction scenario, RPCR plus adiabatic compression appears to the best model to explain the γ-ray emission in SNR/MC. RPCR suggests a transition of seed particles in SNR evolution, which is from thermal injected seed particles in young SNR to the pre-existing ambient CRs in middle-aged SNR. In the transition, the spectrum of SNR/MC gradually changes from an exponential like cut-off at very high energy to a power-law profile as the remnant age increases.

The steep GeV spectrum and non-detection of TeV emission in W44 imply that it may be a good example of SNRs in such transition. It is likely that the thermal injected particles in W44 still dominate the seed particles, while ambient CRs already become the dominant component of seed particle in the other middle-aged SNRs. Multiwavelength observation of both young and middle-aged SNRs is needed in future to test the idea of transition in seed particle.

5 DISCUSSION

In this work, we study the γ-ray emission from 10 SNR/MC in the First Fermi SNR catalogue with a focus on the spectral shape. We compare the γ-ray data available in literature with the model prediction from widely used escaping scenario and direct interaction to obtain deeper insight into the physical origin of γ-ray emission.

In the escaping scenario, the shape of model spectrum strongly depends on the spatial configuration of SNR/MC. Since the observed 10 SNR/MC are likely in different spatial configuration, a variety of γ-ray spectra with different shape are expected in observation, which however is not seen in current data. Moreover, in order to reproduce the γ-ray data, CR protons with energy ≲ 1 GeV must be able to escape from the remnant and then diffuse into the nearby MCs, which however is still an open question. Two best examples of illuminated MCs discussed in the literature are γ-ray bright MCs external to W28 and W44. We argue that the γ-ray emission from them is possibly attributed to or at least partly attributed to the ambient CRs and other associated sources. In summary, the γ-ray spectra from 10 SNR/MC are inconsistent with the prediction of escaping scenario statistically. The inconsistency not only challenges the escaping scenario but also challenges our understanding of CR escaping in SNRs. It may imply that the free escape boundary widely adopted in the escaping scenario is not a good recipe to describe the spatial distribution of escaping CRs. We have to keep in mind that the free escape boundary was originally introduced in non-linear DSA to achieve a self-consistent treatment of particle acceleration at the shock front. It is not designed to study the propagation of escaping CRs in the ISM surrounding an SNR. If we instead assume a finite escape probability at the forward shock for particle with all energies, then the sharp low-energy cut-off in the escaping CR spectrum will disappear, which may relieve the tension between γ-ray data and the escaping model.

In the direct interaction scenario, the model involving RPCR and adiabatic compression is able to explain the γ-ray emission from most SNR/MC. RPCR suggests a transition in seed particle, which is from thermal injected seed particle in young SNRs to pre-existing ambient CRs in middle-aged SNRs. The transition is likely a natural result of declining shock velocity in SNR evolution. Because in old SNRs with slow shock the thermal particles behind the shock front are shifted to lower energy, which makes the thermal injection less efficient. Re-acceleration of ambient CRs is potentially very interesting as it may help explain the anomaly detected by PAMELA and AMS in the ambient CR spectrum (e.g. Adriani et al. 2011). For example, Thoudam & Hörandel (2014) propose that re-acceleration of ambient CRs at slow SNR shocks is able to explain the observed spectral break in proton and Helium spectra at a few hundred GV rigidity. Multiwavelength observation is needed in future to investigate whether such transition exists in SNR evolution. We also propose that radiative SNRs without MC interaction are able to produce a significant amount of γ-ray emission. A good candidate is S147. With accumulated Fermi data and CTA in future, we expect to detect more remnants like S147. Through the discussion, we assume strong shock condition. If the slow shock in middle-aged SNRs happens to have a low Mach number, then it can reduce the compression ratio and steepen the steady state DSA spectrum, which needs to be investigated in future work.

ACKNOWLEDGEMENTS

XT is very grateful to Roger Chevalier, Patrick Slane, Andrew Strong, Shigehiro Nagataki, Shiu-Hang Lee, and the anonymous referee for constructive comments, which help to improve the manuscript significantly. XT would also like to thank Yang Chen, Xiao Zhang, Yasunobu Uchiyama, and Eugene Churazov for reading the manuscript.

Footnotes

REFERENCES

APPENDIX A: π⁰-decay emisison

The π⁰-decay emission from proton–proton interaction is calculated with the parametrized γ-ray production cross-sections dσ/dE_γ developed in Kafexhiu et al. (2014). At low energy, the model is fitted with experimental data while at high energy it is tested with public available code GEANT4, PYTHIA 8, SIBYLL 2.1, and QGSJET-I. The four public codes predict slightly different results at high energy. In this work, we apply the formula of dσ/dE_γ fitted with GEANT 4 results to do the calculation. The analytical formula is found to be accurate within |$20\hbox{ per cent}$| accuracy from the kinematic threshold (280MeV) to PeV energies. We notice that the γ-ray production cross-sections provided by Kafexhiu et al. (2014) are not smooth at some connecting points. But the resulted π⁰-decay emission seems to be unaffected. Emission from the secondary electrons are neglected for simplification in this paper.

APPENDIX B: RUNAWAY CR SPECTRUM

The runaway CR spectrum f(R, t, p) at a given distance R from the remnant centre and at a given time t since supernova explosion depends on several physical processes including SNR evolution, DSA, and the recipe for CR escaping, which is assumed to be free escape boundary here. The resulted γ-ray emission due to interaction between runaway CRs and particles in MCs further depends on the CR propagation in ISM and spatial distribution of MCs, i.e. shape and density profile of MCs.

In the upper panel of Fig. 3, we present N(E, t, L₁, L₂) as a function of proton kinetic energy E for different remnant age t and different distance L₁ of MCs. Through the paper, we fix the thickness of MCs to be 5 pc, i.e. L₂ − L₁= 5 pc. In the calculation, we assume R_sedov = 2.1pc and t_sedov= 210 yras in Ohira et al. (2011) for simplification. N(E, t, L₁, L₂) and N(p, t, L₁, L₂) are related by N(E, t, L₁, L₂) = N(p, t, L₁, L₂)/βc, where β is proton velocity in c.