Interaction between Molecular Clouds and MeV-TeV Cosmic-ray Protons Escaped from Supernova Remnants

TL;DR

This paper presents a model explaining how MeV and GeV-TeV cosmic-ray protons interact with molecular clouds around supernova remnants, accounting for observed iron line emissions and gamma-ray spectra, supporting a hadronic origin.

Contribution

The study introduces an analytical model linking MeV and GeV-TeV cosmic-ray behaviors with observed emissions around supernova remnants, highlighting the role of CR escape and interaction.

Findings

Model reproduces iron line intensities and gamma-ray spectra for W28 and W44.

Supports hadronic scenario for gamma-ray emissions from SNRs.

Shows MeV CRs are confined during SNR-cloud collision, then leak into clouds.

Abstract

Recent discovery of the X-ray neutral iron line (Fe I K\alpha at 6.40 keV) around several supernova remnants (SNRs) show that MeV cosmic-ray (CR) protons are distributed around the SNRs and are interacting with neutral gas there. We propose that these MeV CRs are the ones that have been accelerated at the SNRs together with GeV-TeV CRs. In our analytical model, the MeV CRs are still confined in the SNR when the SNR collides with molecular clouds. After the collision, the MeV CRs leak into the clouds and produce the neutral iron line emissions. On the other hand, GeV-TeV CRs had already escaped from the SNRs and emit gamma-rays through interaction with molecular clouds surrounding the SNRs. We apply this model to the SNRs W28 and W44 and show that it can reproduce the observations of the iron line intensities and the gamma-ray spectra. This can be another support of a hadronic scenario for the gamma-ray emissions from these SNRs.

Figures

Equations

$$\frac{\partial f}{\partial t}(t,\mbox{\boldmath$r$},p)-D_{\rm ISM}(p)\Delta f(t,\mbox{\boldmath$r$},p)=q_{s}(t,\mbox{\boldmath$r$},p)\>,$$

$$f_{\mathrm{point}}(t,r,p)=\frac{\exp [-(r/R_{\mathrm{d}}{)}^2]}{{\pi }^{3/2}{R}_{\mathrm{d}}^3}{N}_{\mathrm{esc}}(p),$$

$$R_{\mathrm{d}}(t,p)=\sqrt{4D_{\mathrm{ISM}}(p)[t-t_{\mathrm{esc}}(p)]},$$

$$N_{\rm esc}(p)=\int dt\int d^{3}\mbox{\boldmath$r$}\>q_{\rm s}(t,r,p)\>,$$

$$q_s(p)=\frac{N_{\mathrm{esc}}(p)}{4\pi r^2}\delta [r-R_{\mathrm{esc}}(p)]\delta [t-t_{\mathrm{esc}}(p)].$$

$$f(t,r,p)$$

$$\frac{D_{\mathrm{sh}}(p_{\mathrm{esc}})}{u_{\mathrm{sh}}}\sim l_{\mathrm{esc}},$$

$$p_{\mathrm{esc}}=p_{\max }{\left(\frac{R_{\mathrm{sh}}}{R_{\mathrm{Sedov}}}\right)}^{-\alpha },$$

$$R_{\mathrm{sh}}(t)=R_{\mathrm{Sedov}}{\left(\frac{t}{t_{\mathrm{Sedov}}}\right)}^{2/5},$$

$$R_{\mathrm{esc}}(t)=(1+\kappa )R_{\mathrm{sh}}(t).$$

$$t_{\mathrm{esc}}(p)=t_{\mathrm{Sedov}}{\left(\frac{p}{p_{\max }}\right)}^{-5/(2\alpha )}.$$

$$N_{\mathrm{esc}}(p)\propto p^{-(s+\beta /\alpha )},$$

$$D_{\mathrm{ISM}}(p)=10^{28}\chi {\left(\frac{pc}{10\mathrm{GeV}}\right)}^{\delta }\mathrm{c}{\mathrm{m}}^2{\mathrm{s}}^{-1}$$

$$N_{\mathrm{sh}}(t,p)=N_{\mathrm{esc}}(p_{\mathrm{esc}}(t)){\left(\frac{p}{p_{\mathrm{esc}}(t)}\right)}^{-s},$$

$$n_{\mathrm{CR},\mathrm{sh}}(t,p)=\frac{N_{\mathrm{sh}}(t,p)}{V_{\mathrm{c}}},$$

$$I_{6.4keV}=\frac{1}{4\pi }\int dE{\sigma }_{6.4keV}(E)v_{\mathrm{CR}}(E)n_{\mathrm{H}}\int dx{n}_{\mathrm{CR}}(E,x),$$

$$I_{6.4keV}$$

$$R_{\mathrm{esc}}(t_{\mathrm{obs}})=L_1,$$

Full text

1 Introduction

Supernova remnants (SNRs) have been thought to be the site where cosmic-rays (CRs) with an energy of $E\lesssim 10^{15.5}$ eV (the knee energy) are accelerated. The most plausible process of the CR acceleration is a diffusive shock acceleration (DSA) at their shock front (Bell, 1978; Blandford & Ostriker, 1978; Drury, 1983). An excess in GeV-TeV gamma-rays have been observed for SNRs associated with molecular clouds (e.g. Abdo et al. (2009); Aharonian et al. (2008)), which is believed to be evidence that CRs are actually accelerated at SNRs. However, there has been a debate on whether the origin of the gamma-rays is leptonic or hadronic. If it is leptonic, these signals could be caused by bremsstrahlung or inverse Compton scattering of electrons. However, recent detections of the characteristic pion-decay feature in the gamma-ray spectra strongly suggest that the origin should be hadronic (e.g. Ackermann et al. (2013); Jogler & Funk (2016)). However, since gamma-rays are produced only by CR protons with energies of $E\gtsim$ GeV, the gamma-ray observations cannot probe lower-energy ($\sim$ MeV) protons that are likely to be accelerated at the same time. If the existence of the lower-energy protons are confirmed, it could be a further support of the hadronic scenario.

The lower-energy protons could be studied through ionization signatures in molecular gas, because those protons are very efficient in ionizing molecular gas (Schuppan et al., 2012; Krause et al., 2015). Alternatively, they could be investigated through X-ray neutral iron line emissions (Fe I K$\alpha$ at 6.40 keV). Recently, Nobukawa et al. (2018) actually detected the iron line emissions from five SNRs interacting with molecular clouds from Suzaku archive data. However, the line emissions could be produced not only by low-energy protons but also electrons and X-rays. Nobukawa et al. (2018) concluded that protons are most likely because of the observed large equivalent width of the line and non-existence of nearby X-ray sources. The iron line emissions were also reported by Okon et al. (2018) for the SNR W28 (see also Sato et al. (2014, 2016); Saji et al. (2018)).

In the hadronic model, two scenarios have been considered for the gamma-ray emissions associated with molecular clouds. One is the direct interaction scenario, in which an SNR directly interacts with molecular clouds (e.g. Bykov et al. (2000); Lee et al. (2015)). In particular, reacceleration and/or compression of Galactic background CR protons may boost their energy and create gamma-ray emissions from molecular clouds (Uchiyama et al., 2010; Tang & Chevalier, 2014; Cardillo et al., 2016; Tang, 2019). However, this scenario could face difficulty in explaining the neutral iron line emissions from MeV protons. This is because the ionization cooling time of the low-energy protons is very short in high-density molecular clouds and thus it is unlikely that the clouds contain those protons as background particles. Thus, in this study we focus on another scenario called the escaping scenario, in which the molecular clouds passively interact with the CR particles escaping from an adjacent SNR (e.g. Aharonian & Atoyan (1996); Fujita et al. (2009); Gabici et al. (2009); Li & Chen (2010)). We aim to explain both the gamma-ray and neutral iron line emissions based on the escaping scenario for the first time.

This paper is organized as follows. In section 2, we explain our model about the escape of CR protons from SNRs and interaction between the CRs and molecular clouds. In section 3, we apply our model to two SNRs (W28 and W44) and show that both neutral iron line emissions and gamma-ray spectra can be explained. The results are discussed in section 4. Conclusions are given in section 5. Hereafter, we refer to CR protons as CRs.

2 Models

2.1 Distribution of high-energy CRs escaped from an SNR

In this subsection, we summarize the derivation of the distribution function of high-energy ($\gtrsim$ GeV) CRs escaped from an SNR based on the model given by Ohira et al. (2011).

We solve a diffusion equation

[TABLE]

where $r$ is the position, $p$ is the CR momentum, f(t,\mbox{\boldmathr},p) is the distribution function, $D_{\rm ISM}(p)$ is the diffusion coefficient in the interstellar medium (ISM) around the SNR, and q_{\rm s}(t,\mbox{\boldmathr},p) is the source term of CRs.

In the following, we assume that the SNR is spherically symmetric and $r$ is the distance from the SNR center. Moreover, we assume that CRs with a momentum $p$ escape from the SNR at $t=t_{\rm esc}(p)$ (Ptuskin & Zirakashvili, 2005; Ohira et al., 2010). In the case of a point source, the source term is give by q_{\rm s}=N_{\rm esc}(p)\delta(\mbox{\boldmathr})\delta[t-t_{\rm esc}(p)], and the solution is

[TABLE]

where

[TABLE]

and

[TABLE]

which is the spectrum of the whole escaped CRs.

In realty, CRs escape from the surface of the SNR, $R_{\rm esc}(p)$, and the source term should be,

[TABLE]

For this source term, the solution of equation (1) can be derived using equation (2) as the Green function:

[TABLE]

We need to specify the escape time $t_{\rm esc}(p)$, the radius $R_{\rm esc}(p)$, and the spectrum $N_{\rm esc}(p)$. We assume that the SNR is in the Sedov phase and CRs are accelerated through a DSA. Thus, CRs are scattered back and forth across the shock front by magnetic turbulence during the acceleration. The diffusion coefficient around the shock, $D_{\rm sh}(p)$, is expected to be much smaller than $D_{\rm ISM}(p)$ for a given $p$, and the diffusion length of the CRs is $\sim D_{\rm sh}(p)/u_{\rm sh}$, where $u_{\rm sh}$ is the velocity of the shock front. We assume that if the CRs cross an escape boundary outside the shock front, they escape from the SNR. Thus, the momentum of escaping CR, $p_{\rm esc}$, is given by

[TABLE]

where $l_{\rm esc}$ is the distance of the escape boundary from the shock front. We adopt a geometrical confinement condition $l_{\rm esc}=\kappa R_{\rm sh}$ and assume that $\kappa=0.04$ (Ptuskin & Zirakashvili, 2005; Ohira et al., 2010). The escape momentum $p_{\rm esc}$ is expected to be a decreasing function of the shock radius. Here, we adopt a phenomenological power-law relation:

[TABLE]

where $p_{\rm max}$ and $R_{\rm Sedov}$ are the escape momentum and the shock radius at the beginning of the Sedov phase ($t=t_{\rm Sedov}$), respectively. Following Ohira et al. (2011), we assume that the index is $\alpha=6.5$, which well reproduces gamma-ray spectra of SNRs. While Ohira et al. (2011) assumed that $p_{\rm max}c=10^{15.5}$ eV (the knee energy), we assume that $p_{\rm max}c<10^{15.5}$ eV and treat it as a parameter because there has been no direct evidence that CRs are accelerated up to the knee energy at SNRs (e.g. Gabici (2017)).

Since we assumed that the SNR is in the Sedov phase, the shock radius is represented by

[TABLE]

and the escaping radius is given by

[TABLE]

We assume that $R_{\rm Sedov}=2.1$ pc and $t_{\rm Sedov}=210$ yr following Ohira et al. (2011). Eliminating $R_{\rm sh}$ from equations (8) and (9) and replacing $p_{\rm esc}$ and $t$ with $p$ and $t_{\rm esc}$, respectively, we obtain

[TABLE]

We assume that the CR spectrum at the shock front is always represented by a single power-law $\propto p^{-s}$ and the number of CRs in the momentum range $(m_{\rm p}c,m_{\rm p}c+dp)$ in the SNR is $K(R_{\rm sh})dp\propto R_{\rm sh}^{\beta}$, where $m_{\rm p}$ is the proton mass. The factor $K(R_{\rm sh})$ corresponds to the normalization of the CR spectrum confined in the SNR. If we assume a thermal leakage model for CR injection, the index is $\beta=3(3-s)/2$ (Ohira et al., 2010). Based on these assumptions, the spectrum of the escaped CRs ($p>p_{\rm esc}$) is written as

[TABLE]

(Ohira et al., 2010). Note that the spectrum of the whole escaped CRs ($p>p_{\rm esc}$) is represented by $\propto p^{-(s+\beta/\alpha)}$ regardless of time [see equation (5)]. We determine the normalization of equation (12) from the total energy of the escaped CRs with $pc>1$ GeV ($E_{\rm tot,CR}$), which is treated as a parameter.

For the diffusion coefficient in the ISM, we assume the following form,

[TABLE]

(Ohira et al., 2011). In this study, we assume Kolmogorov-type turbulence ($\delta=1/3$), which is theoretically motivated and close to the values estimated based on recent observations ($\delta\sim 0.4$; Evoli et al. (2015); Genolini et al. (2015)). The constant $\chi(\leq 1)$ is introduced because the coefficient around SNRs can be reduced by waves generated through the stream of escaping CRs (e.g. Fujita et al. (2010, 2011)). In this study, we fix it at $\chi=0.5$.111The diffusion coefficient is related to magnetic fluctuations $\delta B$ as in $D_{\rm ISM}/D_{\rm Bohm}\sim(B/\delta B)^{2}$, where $D_{\rm Bohm}=(1/3)r_{\rm L}v_{\rm CR}$ is the Bohm diffusion coefficient, $r_{\rm L}$ is the gyro-radius, $v_{\rm CR}$ is the velocity of the particle, and $B$ is the background magnetic field (e.g. Roh et al. (2016)). Equation (13) indicates that $\delta B/B\sim 0.002$ for $B=3\rm\>\mu G$ and $pc\sim 1$ GeV.

2.2 Low-energy CRs interacting with molecular clouds and iron line emissions

The 6.4 keV neutral iron line emissions have been observed only in the vicinity of SNRs (Nobukawa et al., 2018). Thus, MeV CRs responsible for the line emissions are distributed there. Some of the SNRs show a sign of interaction with molecular clouds through maser emissions (e.g. Pastchenko & Slysh (1974); Wootten (1981); Claussen et al. (1997)). Equation (11) shows that MeV CRs escape from an SNR after GeV CRs escape. For the SNRs we study in section 3 (W28 and W44), Ohira et al. (2011) indicated that while CRs with $E\gtrsim$ GeV have already escaped from the SNRs at this time, MeV CRs have not. In the following, we assume that MeV CRs are still confined around the SNRs when the SNRs contact with the molecular clouds from which the iron line emissions are detected.

The spectrum of the low-energy CRs confined in an SNR is written as

[TABLE]

which is defined for $p<p_{\rm esc}(t)$. We assume that the CRs are confined in a region around the shock front with a width of $W_{\rm sh}\equiv 2\>l_{\rm esc}$ (figure 1a). The number density of the confined CRs is

[TABLE]

where $V_{\rm c}\approx 4\pi R_{\rm sh}^{2}W_{\rm sh}$ is the volume of the confinement region. Note that we do not consider CRs advected into a far downstream region of the shock front ($r<R_{\rm sh}$) because they probably lose their energy through adiabatic cooling.

For the sake of simplicity, we here assume that the shock front is a plane and the molecular cloud is an uniform cuboid, and that the distribution of CRs in the cloud is one-dimensional (figure 1a). We assume that the confined CRs start seeping into the cloud when the escaping boundary ($r=R_{\rm esc}$) contacts the surface of a molecular cloud ($r=r_{\rm MC}$). This is because the CR diffusion coefficient in the cloud is expected be much larger than that in the confinement region due to the wave dumping through collisions between protons and neutral particles (Kulsrud & Cesarsky, 1971). The CRs are continuously leaked into the cloud at a rate of $n_{\rm CR,sh}u_{\rm sh}$ per unit area of the shock front until the confinement region passes the surface of the cloud ($R_{\rm esc}-W_{\rm sh}=r_{\rm MC}$).

The photon number intensity of the neutral iron line is given by

[TABLE]

where $E$, $\sigma_{\rm 6,4keV}(E),v_{\rm CR}(E),n_{\rm H},n_{\rm CR}(E,x)$ are the kinetic energy of the CRs, the cross section to produce the iron line at 6.4 keV, the CR velocity, the number density of hydrogens in the molecular cloud, and the CR density in the cloud, respectively. The depth of the cloud in the direction of line of sight is represented by $x$. In figure 1a, we assume that the angle between $x$-direction and $r$-direction is zero ($\theta=0$), and $x=0$ corresponds to $r=r_{\rm MC}$. For the the cross-section $\sigma_{{\rm 6,4keV}}(E)$, we use the one for the solar metallicity and $1<E<10^{4}$ MeV calculated by Tatischeff et al. (2012). We assume that $\sigma_{\rm 6.4keV}(E)=0$ for $E>10^{4}$ MeV and $E<1$ MeV, which does not affect the results.

If the injection and the cooling of CRs are balanced in the cloud, the column density of the CRs is written as $n_{\rm CR,sh}u_{\rm sh}t_{\rm cool}$, where $t_{\rm cool}$ is the cooling time of the CRs and is calculated using the ionization loss rate given by Mannheim & Schlickeiser (1994). The column density corresponds to the second integral of equation (16). Thus, the line intensity is represented by

[TABLE]

This equation is correct if the injection of CRs into the molecular cloud is endless. However, the width of the confinement region $W_{\rm sh}$ is finite and the CR column density $n_{\rm CR,sh}u_{\rm sh}t_{\rm cool}$ cannot be larger than $n_{\rm CR,sh}W_{\rm sh}$. Since we do not know how deep the confinement region is immersed in the cloud at present, we simply assume that the region is half immersed ($d_{\rm MC}=W_{\rm sh}/2$ in figure 1a). Thus, equation (17) is modified as

[TABLE]

where $t^{\prime}_{\rm cool}(E)=\min[t_{\rm cool}(E),0.5\>t_{\rm pass}]$ and $t_{\rm pass}\equiv W_{\rm sh}/u_{\rm sh}$ is the time scale in which the confinement region passes the surface of the molecular cloud. This means that the CRs that were originally in the overlapped region (the shaded region in figure 1a) when $d_{\rm MC}=W_{\rm sh}/2$ was satisfied have escaped into the cloud.

However, this correction is not significant. The time scale in which the confinement region passes the surface of the molecular cloud ($t_{\rm pass}$) is $\sim 20$% of the age of the SNR ($\sim 10^{4}$ yr) for the parameters we choose in section 3. While this time scale is too short to change the gamma-ray spectrum produced by $\gtrsim$ GeV CRs, it is much larger than the cooling time of MeV CRs. For example, CRs with $E\sim 10$ MeV are most effective to create the iron line emissions due to a large $\sigma_{\rm 6.4keV}(E)$; they have a cooing time of $t_{\rm cool}\lesssim 100$ yr for $n_{\rm H}>1000\rm\>cm^{-3}$, which means that $t^{\prime}_{\rm cool}=t_{\rm cool}$ at $E\sim 10$ MeV.

For the sake of simplicity, we assume that CRs that have entered clouds do not escape from the clouds before they lose their energy through the rapid cooling. This may be realized if magnetic fields are oriented in the clouds so that they trap the CRs. In other words, the iron line emissions are produced only in the clouds where the fields are properly distributed. Under this assumption, the intensity represented by equation (18) does not depend on the details of the molecular could such as the total length in the direction of $x$ or the CR diffusion coefficient inside it. If the viewing angle $\theta$ is not zero (figure 1b), we expect that $I_{\rm 6.4keV}$ becomes larger and $I_{\rm 6.4keV}(\theta)\sim I_{\rm 6.4keV}(\theta=0)/\cos\theta$. From now on, we assume that $\theta$ is not too close to $90^{\circ}$ unless otherwise mentioned.

2.3 Gamma-ray emissions from molecular clouds

Gamma-rays are produced through $pp$-interaction between CRs and hydrogens in molecular clouds. We assume that molecular clouds with density $n_{\rm H}$ are distributed in a shell region between $r=L_{1}\approx R_{\rm esc}$ and $r=L_{2}$ with a filling factor $f_{\rm gas}$ (figure 1b). The current time $t_{\rm obs}$ is given by

[TABLE]

and equations (9) and (10)222We implicitly assumed that the current time is given by $R_{\rm sh}(t_{\rm obs})=L_{1}$ when we calculate the neutral iron line intensity [equation (18)]. We ignore the difference of the current times when we calculate gamma-ray spectra.. CRs with $p>p_{\rm esc}$ have escaped from the SNR, and their distribution function is given by equation (6). The momentum spectrum for the CRs with $p<p_{\rm esc}$ that have seeped into the cloud is given by $f_{\rm gas}N_{\rm sh}/2$ if cooling is ignored. The factor of two comes from the assumption that the CRs that were originally in the overlapped region (the shaded region in figure 1a) when $d_{\rm MC}=W_{\rm sh}/2$ was satisfied have escaped into the cloud (section 2.2). We calculate the gamma-ray spectra using a model by Kamae et al. (2006) and Karlsson & Kamae (2008).

3 Results

We apply our model to the SNRs W28 and W44. We chose these objects because both neutral iron line emissions and gamma-ray emissions have been detected (Nobukawa et al., 2018). In particular, detailed gamma-ray spectra are available for these objects. Our procedure is as follows. (1) Since the distribution and the mass of molecular clouds are determined mainly through radio observations, we fix the parameters for the clouds at their observed values (table 3.2). (2) Using the observed gamma-ray spectra, we constrain our model parameters for CRs through $\chi^{2}$ fitting (table 3.2). (3) From the fitting results, we estimate the intensities of neutral iron line emissions (table 3.2) and compare them with the observed intensities.

3.1 W28

W28 is a middle-aged SNR from which gamma-rays have been observed in the GeV (Abdo et al., 2010; Hanabata et al., 2014; Cui et al., 2018) and the TeV bands (Aharonian et al., 2008). Since previous studies have shown that the distance to the SNR is $d\sim 2$ kpc (e.g. Goudis (1976); Velázquez et al. (2002)), we assume that $d=2$ kpc in this study (table 3.2).

We focus on the northern gamma-ray component (HESS J1801-233; Aharonian et al. (2008)), which appears to be associated with the neutral iron line emissions (Nobukawa et al., 2018). For W28, we assume that $L_{1}=12$ pc, $L_{2}=15$ pc and $f_{\rm gas}=0.1$ (table 3.2). They are estimated from the distribution of molecular gas and the gamma-ray images (figure 2 in Aharonian et al. (2008)), assuming that the cloud is rather spherical. We fix the mass of the molecular cloud at the observed value ($M_{\rm gas}\sim 5\times 10^{4}\>M_{\odot}$; Aharonian et al. (2008)). For these parameters, the gas number density is $n_{\rm H}=3000\>\rm cm^{-3}$. The current time is $t_{\rm obs}=1.5\times 10^{4}$ yr [equation (19) and table 3.2].

The gamma-ray spectrum is sensitive to the total CR energy ($E_{\rm tot,CR}$), the maximum momentum of CRs ($p_{\rm max}$), and the index of the CR energy spectrum ($s$) at the shock front. Thus, we fit the observed spectrum with our model by varying these three parameters; the other parameters are fixed. The results are shown in table 3.2. The escape momentum $p_{\rm esc}$ and the iron line intensity $I_{\rm 6.4keV}$ can be obtained as a result of the fit (table 3.2). Figure 2 shows that the best-fit model well reproduces the Fermi (Abdo et al., 2010; Cui et al., 2018) and the HESS observations (Aharonian et al., 2008). The iron line intensity is $I_{\rm 6.4keV}=0.07^{+0.01}_{-0.05}\rm\>photons\>s^{-1}\>cm^{-2}\>sr^{-1}$ (table 3.2), which is consistent with $I_{\rm 6.4keV}=0.10\pm 0.05\rm\>photons\>\>s^{-1}\>cm^{-2}\>sr^{-1}$ obtained by Nobukawa et al. (2018)333Nobukawa et al. (2018) did not represent $I_{\rm 6.4keV}$ for individual SNRs in their paper.. On the other hand, Okon et al. (2018) observed more outside regions (closer to the rim) of the SNR. Their obtained values of $I_{\rm 6.4keV}$ are generally larger than that reported by Nobukawa et al. (2018). In particular, for the region where the shock is interacting with clouds or the rim of the SNR (region 1 in their paper), the intensity is $I_{\rm 6.4keV}=0.48\pm 0.28\rm\>photons\>\>s^{-1}\>cm^{-2}\>sr^{-1}$ if the contribution from the Galactic ridge X-ray emission is subtracted. This is probably because the angle between the line of sight and the radial direction of the SNR is close to $\theta=90^{\circ}$ (figure 1b).

3.2 W44

W44 is another middle-aged SNR from which gamma-rays have been observed in the GeV band; the decrement below $\sim 200$ MeV suggests a hadronic origin (Abdo et al., 2010; Ackermann et al., 2013; Cardillo et al., 2014). On the other hand, TeV gamma-rays have not been detected (Buckley et al., 1998; Aharonian et al., 2002). Since the distance has been estimated to be $d\sim 3$ kpc (Caswell et al., 1975; Wolszczan et al., 1991), we assume that $d=3$ kpc. The mass of cold gas directly associated with W44 is rather uncertain. Yoshiike et al. (2013) estimated that the total molecular mass around W44 is $\sim 4\times 10^{5}\>M_{\odot}$. However, the radio and gamma-ray images show that only part of the molecular gas seems to be responsible for the gamma-ray emissions (their figures 3 and 7). From the overlapped area of the radio and gamma-ray emitting regions, we assume that $M_{\rm gas}=4\times 10^{4}\>M_{\odot}$ and $f_{\rm gas}=0.1$ (table 3.2), which result in $n_{\rm H}=2400\rm\>cm^{-3}$. The mass is larger than that of the shocked gas ($\sim 7\times 10^{3}\>M_{\odot}$; Yoshiike et al. (2013)), which should be the minimum mass of the gas directly associated with the SNR. We also assume that $L_{1}=12$ pc and $L_{2}=15$ pc by reference to the gamma-ray image (table 3.2). The uncertainties of the parameters are discussed in the next section. The current time is $t_{\rm obs}=1.5\times 10^{4}$ yr (table 3.2).

In figure 3, we present the gamma-ray spectrum of W44. Our model results are consistent with the Fermi results (Ackermann et al., 2013). The peak of the spectrum is attributed to the break of the CR momentum spectrum at $p=p_{\rm esc}$ [equations (12) and (14); see also Ohira et al. (2011)]. The gamma-ray energy at the peak is larger than that of W28 (figure 2), which reflects the larger value of $p_{\rm esc}$ (table 3.2). The iron line intensity is $I_{\rm 6.4keV}=0.10^{+0.04}_{-0.02}\rm\>\>photons\>\>s^{-1}\>cm^{-2}\>sr^{-1}$ (table 3.2), which is consistent with $I_{\rm 6.4keV}=0.15\pm 0.08\rm\>photons\>\>s^{-1}\>cm^{-2}\>sr^{-1}$ obtained by Nobukawa et al. (2018).

4 Discussion

We have shown that observational results of W28 and W44 can be explained by an CR escaping scenario for SNRs. In our model, the SNRs that show neutral iron line emissions are interacting with surrounding molecular clouds. While CRs with $E\gtrsim$ GeV have already escaped from the SNRs, those with $E\sim$ MeV are now leaking into the clouds and producing the iron line emissions.

Our model constrains parameters for an SNR by the iron line and gamma-ray emissions, respectively. From the iron line intensity $I_{\rm 6.4keV}$, the product of $n_{\rm CR,sh}(E\sim 10~{}{\rm MeV})u_{\rm sh}$ can be obtained from equation (18). Note that since the ionization cooling time satisfies $t_{\rm cool}n_{\rm H}=$ const (e.g. Mannheim & Schlickeiser (1994)), the iron line intensity is not dependent of $n_{\rm H}$ as long as $t_{\rm cool}<0.5\>t_{\rm pass}$, which is fulfilled when $E$ is relatively small ($E\lesssim 100$ MeV). On the other hand, observations of the gamma-ray luminosity determines the product of $n_{\rm CR}(\gtrsim{\rm GeV})M_{\rm gas}$ because the clouds are the targets of CRs in $pp$-interaction.

The parameters regarding molecular clouds, especially for W44 (section 3.2), have some uncertainties, and there is a kind of degeneracy in terms of gamma-ray luminosity. For example, if $M_{\rm gas}$ is doubled by doubling the gas density $n_{\rm H}$, the total CR energy $E_{\rm CR,tot}$ needs to be halved. This is because the CR density $n_{\rm CR}$ needs to be halved to be consistent with the observed gamma-ray luminosity. In this case, $I_{\rm 6.4keV}$ is reduced but it is not exactly halved because of the contribution of CRs with $E\gtrsim 100$ MeV for which the relation of $t_{\rm cool}n_{\rm H}=$ const is not valid. The intensity $I_{\rm 6.4keV}$ is marginally consistent with observations for W28 and W44. Similarly, if $M_{\rm gas}$ is doubled by doubling the filling factor $f_{\rm gas}$, $E_{\rm CR,tot}$ needs to be halved. As a result, the CR density $n_{\rm CR}$ is halved, and $I_{\rm 6.4keV}$ is halved. These could be used to verify our model by future observations of cold gas. For example, given the currently observed gamma-ray luminosity and $I_{\rm 6.4keV}$, it is unlikely that $M_{\rm gas}$ is larger than those we assumed by an order of magnitude, although there is an uncertainty about the dependence of $I_{\rm 6.4keV}$ on the angle $\cos\theta$ (section 2.2). The lower limit of $M_{\rm gas}$ can be constrained by the upper limit of $E_{\rm CR,tot}$ because $E_{\rm CR,tot}$ cannot exceed the explosion energy of an supernova ($\sim 10^{51}$ erg). Thus, it is unlikely that $M_{\rm gas}$ is smaller than those we assumed by an order of magnitude.

For $n_{\rm H}>1000\rm\>cm^{-3}$, the cooling time of CRs is very short ($t_{\rm cool}\lesssim 100$ yr at $E\sim 10$ MeV). This means that MeV CRs that enter the cloud almost immediately lose their energy. Thus, the neutral iron line emissions can be observed in a time scale of $\sim t_{\rm pass}$. Table 3.2 shows that $t_{\rm pass}$ is 20% of the current time $t_{\rm obs}$ and thus the duration is not extremely small compared with the age of the SNRs. In other words, the possibility of observing the iron line emissions is not tiny. Moreover, if multiple clouds are randomly located around the SNR, the iron line emissions could blink on and off as the shock front passes the clouds. We note that we assumed that the fraction of the confinement region ($W_{\rm sh}/R_{\rm sh}\propto\kappa$) for MeV CRs is the same as that for GeV–TeV CRs. If this is not the case, the duration $t_{\rm pass}$ and the line intensity $I_{\rm 6.4keV}$ can change [equations (15) and (18)]. For example, if $W_{\rm sh}$ is doubled, $t_{\rm pass}$ is doubled and $I_{\rm 6.4keV}$ is halved.

5 Conclusion

We have shown that 6.4keV neutral iron line emissions and gamma-ray emissions from SNRs can be explained by an CR escaping scenario for SNRs. In this model, the SNRs with the iron line emissions are interacting with surrounding molecular clouds. We assume that CRs are accelerated at the SNR with a single power-law spectrum. When the SNR comes into contact with the clouds, MeV CRs are still confined in the SNR. They gradually leak into the clouds and produce the iron line emissions through interaction with irons in the clouds. On the the hand, the CRs with $E\gtrsim$ GeV have already escaped from the SNR at the contact.

We applied this model to the SNRs W28 and W44 and showed that both the observed iron line intensities and the gamma-ray spectra can be reproduced. These support a hadronic scenario for the gamma-ray emissions from the SNRs.

{ack}

We thank the anonymous referee, whose comments improved the clarity of this paper. This work was supported by MEXT KAKENHI No. 18K03647 (Y.F.), JP16J00548 (K.K.N.), JP15H02070 (H.M.), and 16K17702 (Y.O.).