Detection of gamma-rays around SNR HB9 and its implications for the diffusive shock-acceleration history

Abstract

We analyze the GeV gamma-ray emission data from the vicinity of the supernova remnant (SNR) HB9 (G160.9+2.6) from the Fermi-LAT 12 yr observations to quantify the evolution of diffusive shock acceleration (DSA) in the SNR. In the vicinity of HB9, there are molecular clouds whose locations do not coincide with the SNR shell in the line of sight. We detect significant gamma-ray emissions above 1 GeV spatially coinciding with the two prominent cloud regions, as well as emission from the SNR shell, the latter of which is consistent with the results of previous studies. The energy spectrum above 1 GeV in each region is fitted with a simple power-law function of dN/dE ∝ E^−Γ. The fitting result indicates harder spectra, with power-law indices of Γ = 1.84 ± 0.18 and 1.84 ± 0.14, than that at the SNR shell, with Γ = 2.55 ± 0.10. The observed spectra from the cloud regions are found to be consistent with the theoretically expected gamma-ray emissions originating in the protons that escaped from SNR HB9, where particles can be accelerated up to higher energies than those in the shell at present. The resultant diffusion coefficient in the vicinity of the SNR is comparable to that of the Galactic mean.

1 Introduction

The supernova remnant (SNR) is the most common class of cosmic ray (CR) acceleration sites with energies up to ∼ PeV in our Galaxy (see, e.g., Blasi 2013 for a review). However, there is no conclusive observational evidence of SNRs accelerating particles, specifically protons, up to PeV. Two scenarios have been suggested in an attempt to explain the lack of observational evidence of a PeV particle accelerator in SNRs against expectations: (i) an SNR is capable of accelerating particles up to PeV only while it is young, up to an age of 10³ yr (Ptuskin & Zirakashvili 2003) and (ii) such high-energy particles escape from the SNR at an early stage (Gabici et al. 2007). In order to scrutinize the scenario that an SNR has accelerated PeV CRs in the past, it is essential to quantify the evolution of diffusive shock acceleration (DSA) in the SNR. Systematic studies exploring the population of SNRs are expected to provide information on the evolution of the CR spectra of SNRs as a function of the SNR age (Ohira et al. 2011; Schure & Bell 2013; Gaggero et al. 2018). However, quantifying the evolution of DSA with this kind of method is challenging due to the large diversity of the surrounding environments of SNRs (Yasuda & Lee 2019; Suzuki et al. 2020).

Another way to study the evolution of DSA is simultaneous observation of a single SNR, specifically its shell part, and nearby massive molecular clouds (MCs). If a massive cloud exists in the vicinity of the SNR, protons that escape the SNR illuminate the cloud and generate gamma-ray emission via the π⁰-decay process (Aharonian & Atoyan 1996; Gabici et al. 2007). The delay of the timing of the gamma-ray emission from the cloud region from that of the incident proton escape depends on the propagation time and accordingly reflects the particle distribution in the SNR at a specific epoch in the past. Hence, a comparison of the spectra at the SNR shell and nearby clouds observed at roughly the same time enables us to quantify the evolution of the DSA in the SNR.

In observations of SNR W 28, such “delayed” gamma-rays have been detected from three cloud regions with HESS and Fermi-LAT (Aharonian et al. 2008; Hanabata et al. 2014). The observed gamma-ray spectra from these cloud regions, however, do not differ significantly from one another, and hence the maximum energy with the DSA at W 28 has never been successfully measured as a function of the SNR age. Since the angular distances from the SNR center to these individual clouds are almost the same in W 28, the emissions from the clouds should originate in protons accelerated at a similar epoch. The gamma-ray emission from another cloud (named HESS J1801−233) located to the north of W 28 was suggested to originate in protons that have been re-accelerated by a shock-cloud interaction (Cui et al. 2018), in which case the gamma-ray spectrum does not reflect the SNR age. Hanabata et al. (2014) reported detection of a fifth gamma-ray emitting region, named source W, located at the western part of W 28. Although there is a possibility that the emission can be explained by protons escaping from the SNR, no cloud counterpart has been found. Thus far, it is impossible to quantify the particle distribution at multiple epochs with the available observations of W 28.

One concern with the study using delayed gamma-ray spectra is a large dependence on the diffusion coefficient of CRs. The diffusion coefficient may be modified in the vicinity of an SNR due to the effect of self-confinement and/or magnetic-field amplification caused by the generation of turbulent plasma waves (Wentzel 1974; Fujita et al. 2011; D’Angelo et al. 2018). Although delayed gamma-ray emissions have been detected in the vicinity of several SNRs (e.g., Uchiyama et al. 2012), the derived diffusion coefficients had large uncertainties, and so did the estimated amount of particle acceleration at the SNR shell in the past. This is mainly because the SNRs discussed so far are relatively old. For older SNRs, MCs and SNR shells often interact directly and, as a result, the observed data do not provide much information on the diffusion coefficient. Furthermore, even if MCs are significantly separated from SNR shells, the delayed gamma-ray spectra do not differ significantly from those of the current SNR shells, and thus the data are not of much help for restricting the diffusion coefficients (Uchiyama et al. 2012).

Here, we focus on SNR HB9 (G160.9+2.6), which is relatively young (∼6.6 × 10³ yr; Leahy & Tian 2007) compared with other objects where delayed gamma-rays have been observed. HB9 has two additional advantages for this type of study in the DSA evolution. First, there are MCs in the vicinity of this SNR, but their locations do not coincide with the SNR in the line of sight (Sezer et al. 2019), enabling us to simply use the distance between the SNR and clouds to calculate the diffusion time. Secondly, HB9 has observable gamma-ray emission from the SNR shell (Araya 2014), which is essential to estimate the current maximum energy of the accelerated particles at the SNR shock.

HB9 has a large angular size (2° in diameter) and a shell-type radio morphology (e.g., Leahy & Roger 1991) with a radio spectral index of α = −0.47 ± 0.06 in frequencies between 408 and 1420 MHz (Leahy & Tian 2007). Leahy and Aschenbach (1995) estimated its distance and age to be 1.5 kpc and (8–20) × 10³ yr, respectively, using X-ray data observed with ROSAT. Its kinematic distance was estimated to be 0.8 ± 0.4 kpc from H i observations, yielding a Sedov age estimate of 6.6 × 10³ yr (Leahy & Tian 2007). Sezer et al. (2019) found a H i shell expanding toward the SNR with a velocity ranging between −10.5 and +1.8 km s⁻¹ and derived a kinematic distance of 0.6 ± 0.3 kpc, which is consistent with the estimate by Leahy and Tian (2007). In the gamma-ray band, spatially extended emission was detected from SNR HB9 along with its radio shell with the Fermi-LAT 5.5 yr observations in the energy band above 0.2 GeV (Araya 2014). The gamma-ray spectrum was characterized by a power law with a photon index of 1.44 ± 0.25 accompanied by an exponential cutoff at 1.6 ± 0.6 GeV. The spectrum is explained with emission from inverse Compton (IC) scattering of electrons with a simple power-law energy spectrum with a differential index of 2 and maximum electron energy of 500 GeV. Furthermore, Sezer et al. (2019) analyzed the Fermi-LAT 10 yr data in the energy range between 0.2 and 300 GeV and newly detected a point-like source near the SNR shell, named PS J0506.5+4546. The spectrum of the point source can be characterized by a simple power-law function with an index of 1.90 ± 0.19. The flux was (6.59 ± 3.47) × 10⁻¹⁰ cm⁻² s⁻¹. PS J0506.5+4546 was, however, suggested not to be related to the SNR shell, and its origin is unclear (Sezer et al. 2019).

In this paper, we study the gamma-ray morphology of SNR HB9 and the spectra of the SNR shell and the nearby cloud regions, using 12 yr observations with Fermi-LAT, to quantify the evolution of DSA. In addition, since we suspect that the gamma-ray emission from PS J0506.5+4546 originates in an MC located outside the radio shell, we search for spatial correlation between the ¹²CO (J = 1–0) line and gamma-ray emissions. In section 2, we describe the Fermi-LAT observations and show the results of our analysis. We present the modeling study to interpret the origin of the observed gamma-ray in section 3. The summary and conclusions are given in section 4.

2 Gamma-ray observations

2.1 Data selection and analysis

Fermi-LAT is capable of detecting gamma-rays in an energy range from ∼30 MeV to >300 GeV (Atwood et al. 2009). We analyze its 12 yr data from 2008 August to 2020 August in the vicinity of SNR HB9. The standard analysis software, Science Tools (version v11r5p3),¹ is used. The “Source” selection criteria and instrument responses (P8R2_SOURCE_V6)² are chosen considering a balance between the precision and photon-count statistics. The zenith-angle threshold is set to 90° to suppress the contamination of the background from the Earth rim. We employ the tool gtlike (in the binned mode), using a standard maximum likelihood method (Mattox et al. 1996), for spatial and spectral analyses. We choose a square region of 15 × 15° aligned with the galactic coordinate grid with the center coinciding with that of HB9 (RA = 75|${_{.}^{\circ}}$|25, Dec = 46|${_{.}^{\circ}}$|73) as the region of interest (ROI) for the (binned) maximum likelihood analysis based on Poisson statistics. The pixel size is 0|${_{.}^{\circ}}$|1.

The source spatial-distribution model includes all the sources in the fourth Fermi catalog (4FGL; Abdollahi et al. 2020) within the ROI and the two diffuse backgrounds, the Galactic (gll_iem_v7.fits) and extragalactic (iso_P8R3_SOURCE_V2_v1.txt) diffuse emissions. Regarding the emission of the SNR shell, Araya (2014) and Sezer et al. (2019) concluded that the radio template produced with the 4850 MHz radio continuum data from the Green Bank Telescope (Condon et al. 1994) is the best spatial model. Accordingly, we use the radio template provided in Science Tools as the spatial model. In the fitting of the maximum likelihood analysis, all spectral parameters of HB9 SNR itself, 4FGL sources (Abdollahi et al. 2020) located within 5° from the center of HB9, and the two diffuse backgrounds are allowed to vary freely. We do not use data below 1 GeV in this analysis since the fitting results of delayed gamma-ray emission in this band suffer from systematic uncertainty in the Galactic diffuse background model (Ackermann et al. 2012). Note that the flux of the Galactic diffuse emission in the MC region is larger in this energy range than that of the delayed gamma-ray emission from HB9, which will be estimated in section 3.

The significance of a source is represented in this analysis by a test statistic (TS) defined as −2 log (L₀/L), where L₀ and L are the maximum likelihood values for the null hypothesis and a model including additional sources, respectively (Mattox et al. 1996). The detection significance of the source can be approximated as |$\sqrt{\mathrm{TS}}$| when the number of counts is sufficiently large.

2.2 Morphology study

Figures 1a and 1c show the background-subtracted gamma-ray TS map created from the Fermi-LAT 12 yr data above 1 GeV, where the background model consists of the Galactic and isotropic extragalactic emissions and contributions from known Fermi sources (see the previous subsection). The map is overlaid with cyan contours of the ¹²CO (J = 1–0) line emission from the Dame survey data (Dame et al. 2001), which are integrated over a velocity range between −10.4 and +2.6 km s⁻¹, and also green contours of 1420 MHz radio continuum emission obtained from the CGPS survey with DRAO (Taylor et al. 2003). In figure 1c, no significant emission from the SNR shell was found (green contours in the figure), which is consistent with the fact that the SNR spectrum has a cutoff below 10 GeV (Araya 2014). Figure 1b is a similar gamma-ray TS map to figure 1a, but with the contribution from the SNR shell additionally subtracted, for which we use the 4850 MHz radio-template model as a spatial model and assume a simple power-law spectrum. The gamma-ray excess in figures 1b and 1c appears to be more extended than the point source J0506.5+4546 reported in the previous study and rather spatially coincident with the CO line emission.

Fig. 1.

Gamma-ray TS maps in the vicinity of SNR HB9 observed with Fermi-LAT. All maps are given with square bins of 0.05°, and Gaussian smoothing with a kernel ****σ = 0|${_{.}^{\circ}}$|1 is applied. The energy ranges are 1–500 GeV for panels (a) and (b) and 10–500 GeV for panel (c). The subtracted background consists of (a) and (c) the Galactic and extragalactic diffuse emissions and known gamma-ray source contribution and (b) additionally the estimated gamma-ray emission from the HB9 SNR shell (see text for details). Green contours show the radio emission of SNR HB9 at 1420 MHz with DRAO (Taylor et al. 2003) and are linearly spaced in increments of 0.5 K from 5.5–10.0 K. Cyan contours show the ¹²CO (J = 1–0) line intensity integrated over a velocity range between −10.4 and +2.6 km s⁻¹ and are linearly spaced in increments of 1.0 K km s⁻¹ from 4.5–10.5 K km s⁻¹. The two apparent CO-emission regions (R1 and R2) are indicated. The magenta cross in each panel indicates the position of PS J0506.5+4546 (Sezer et al. 2019). (Color online)

We test whether the gamma-ray (>1 GeV) spatial distribution is correlated with CO line emissions using the likelihood method. The ¹²CO (J = 1–0) line image exhibits two distinctive regions (thus two MCs), designated as R1 and R2 (figure 1). For the test, we create a CO template for each of R1 and R2, which is made from the ¹²CO (J = 1–0) line image (Dame et al. 2001) integrated over a velocity range between −10.4 and +2.6 km s⁻¹ and cut with a threshold of >4.5 K km s⁻¹ (thick cyan contours in figure 1). We assume the gamma-ray spectra of the cloud regions and HB9 SNR shell to follow a simple power-law function of dN/dE∝E^−Γ. Here, (i) the background model (the null hypothesis) consists of the radio template for the HB9 SNR shell as well as the Galactic and extragalactic diffuse emissions and 4FGL sources (Abdollahi et al. 2020). To compare spatial models, we use the Akaike information criterion (AIC; Akaike 1974) defined as 2k − 2 log (L), where k is the number of estimated parameters in the model. Consequently, we find that the model with the smaller AIC value is favored and that the better model gives a larger ΔAIC [= (AIC)₀ − (AIC)_m], where (AIC)₀ and (AIC)_m are the AIC values for the null hypothesis and a model including additional sources, respectively. When we apply the CO template models to the two cloud regions (ii), ΔAIC is found to be 44.4, indicating a significant correlation between the gamma-ray and ¹²CO (J = 1–0) line emissions (table 1). Since the gamma-ray emission from R1 spatially coincides with PS J0506.5+4546 reported by Sezer et al. (2019), we also apply a point-source model to the same position as PS J0506.5+4546, instead of the CO template model of R1 (iii). The resultant AIC is marginally (1.3σ level) improved from the case with the CO template. As a result, the gamma-ray emission from R1 is consistent with PS J0506.5+4546, while R2 is newly detected with a statistical significance of 6.1σ in this study.

Table 1.

Likelihood test results for spatial models and the resultant spectral parameters of the HB9 SNR shell and the two cloud regions using an energy range of 1–500 GeV.*

		|${\rm {R1}}$|			|${\rm {R2}}$|			|${\rm {HB9\,SNR\,itself}}$|
Model	ΔAIC	Flux	Γ	|$\sqrt{\rm {TS}}$|	Flux	Γ	|$\sqrt{\rm {TS}}$|	Flux	Γ	|$\sqrt{\rm {TS}}$|
(i)	0	—	—	—	—	—	—	24.9 ± 2.1	2.54 ± 0.06	13.2
(ii)	44.4	2.0 ± 0.8	1.84 ± 0.18	4.5	4.8 ± 1.2	1.84 ± 0.14	6.1	24.3 ± 2.1	2.55 ± 0.10	12.8
(iii)	58.2	1.3 ± 0.5	1.77 ± 0.18	5.8	4.8 ± 1.2	1.84 ± 0.14	6.1	24.5 ± 2.1	2.54 ± 0.10	13.0

*

The following three models are considered: (i) the background model (null hypothesis); (ii) additionally including gamma-ray emission from R1 and R2 estimated by the CO template model (see text); and (iii) assuming a point-source model, where the position of the source is the same as PS J0506.5+4546 (Sezer et al. 2019), for R1 instead of the CO template model, while the CO template is as in (ii) for R2. All sources of interest are fitted with a simple power-law function of dN/dE ∝ E^−Γ. The “Flux” and “|$\sqrt{\rm {TS}}$|” columns refer to the integral flux in the energy band of 1–500 GeV in 10⁻¹⁰ cm⁻² s⁻¹ and the statistical significance (σ) of each source approximated as |$\sqrt{\mathrm{TS}}$|⁠, respectively.

2.3 Spectral results

We extract Fermi-LAT energy spectra from the radio SNR shell region and the two regions R1 and R2, individually, using the CO template model for the cloud regions. Figure 2 shows the resultant spectra for an energy range between 1 and 500 GeV. We find that the obtained Fermi-LAT spectrum of the SNR shell is consistent with the gamma-ray spectrum reported by Araya (2014) (figure 2). The respective best-fitting parameters are summarized in table 1. Here, a potential concern about the spectrum of R1 is that our assumption of diffuse gamma-ray spatial distribution may not be appropriate, given the spatial proximity with the already identified point source PS J0506.5+4546 (figure 1). However, the discrepancy of the determined spectral properties (flux and index) between the results of the two spatial-distribution models is smaller than the 1σ uncertainty (table 1). Therefore, we conclude that the difference in the results due to the difference between the assumed spatial distributions is not significant. We note that the power-law index of R1 is consistent with that of PS J0506.5+4546 reported by Sezer et al. (2019).

Fig. 2.

Spectral energy distributions measured with the Fermi-LAT (data points) and fitting results in the gamma-ray band. Top: Black and open green stars show the data of the SNR shell obtained in this work and Araya (2014), respectively. Bottom: Red and blue data points represent the spectra of cloud regions R1 and R2, respectively. Also shown by the shaded black region is the fitting result of the SNR shell. (Color online)

3 Spectral modeling and its implications

In this work, we explore the possibility that R1 and R2 are attributed to delayed gamma-rays from MCs illuminated by the CRs accelerated in HB9. In the following section, we calculate the delayed gamma-ray spectra in the MC regions, R1 and R2, using the method described in subsection 3.1, and try to fit the observed spectra of the cloud regions and the SNR shell simultaneously. According to Araya (2014), once the gamma-ray spectrum of the SNR shell is modeled with the hadronic emission, it would require the enormous explosion energy of a supernova or a dense density of the interstellar medium (ISM).³ Hence, we assume that the gamma-ray emission in the SNR shell originates mainly from leptonic processes. On the other hand, we consider two cases for the origin of gamma-ray emissions from the cloud regions: one where the leptonic emission dominates (leptonic-dominated model), and the other where the hadronic emission dominates (hadronic-dominated model) in the GeV band. We present the results of the modeling in subsection 3.2 and discuss the implications for the model parameters in subsection 3.3.

3.1 Modeling delayed gamma-ray emission

The energy spectra of the gamma-ray emissions from the MCs around SNR HB9 are calculated following the method proposed by Gabici, Aharonian, and Casanova (2009) and Ohira, Murase, and Yamazaki (2011). They considered the escape process of the CR protons, which are accelerated by an SNR, from a shock front into interstellar space during the Sedov phase. Escaped CR protons emit high-energy gamma-rays via π⁰ production through p–p collision in MCs in the vicinity of the SNR. In addition, we also discuss the contribution from inverse Compton scattering and non-thermal bremsstrahlung from CR electrons escaping from the SNR. We evaluate later in this subsection the expected gamma-ray fluxes from the R1 and R2 regions in conjunction with the energy spectra of the escaped CR protons and electrons.

The radio continuum contours (figure 1) show a circular symmetric morphology of SNR HB9, suggesting it has probably maintained a spherically symmetric structure throughout its evolution. Therefore, we assume that the SNR shell itself and the escaped CR distribution are spherically symmetric. Based on X-ray observations, the energy of the supernova explosion (E_SN) was estimated to be (0.15–0.30) × 10⁵¹ erg, and the relation between the explosion energy and the hydrogen density of ISM (n_ISM) was also estimated as |$(E_{\rm SN}/10^{51}\, {\rm{erg}}) (n_{\rm ISM}/1\,{\rm cm}^{-3})^{-1} = 5$| (Leahy & Tian 2007). In this paper, we adopt 0.3 × 10⁵¹ erg as the explosion energy, in which case the density is given as 0.06 cm⁻³. If the initial velocity of the blast wave was u_sh,0 = 10⁹ cm s⁻¹, HB9 entered the Sedov phase at the age t_Sedov = 3.6 × 10² yr when its radius was 3.6 pc. Given that the observed radius of HB9 is ∼10 pc, which is larger than 3.6 pc, HB9 must already be in the Sedov phase. Leahy and Tian (2007) estimated the Sedov age of SNR HB9 to be 6.6 × 10³ yr, which we adopt as the age of the SNR in this discussion.

First of all, let us model the energy spectrum of CR protons. During the Sedov phase, the maximum energy of the CR protons accelerated at the SNR shock is determined by the timescale of escape from the acceleration region, and decreases with time (Ptuskin & Zirakashvili 2005; Caprioli et al. 2009; Ohira et al. 2010). The temporal evolution of the maximum energy depends on non-linear processes, such as amplifying the magnetic field at the shock, which are still theoretically unclear. Thus, we adopt a phenomenological power-law dependence of the cutoff energy E_esc of the proton spectrum at the SNR shock on the age t of the SNR, as discussed by Gabici, Aharonian, and Casanova (2009) and Ohira et al. (2010):

$$\begin{equation} E_{\rm esc}(t)=E_{\max }\!\left(\frac{t}{t_{\rm Sedov}}\right)^{-\alpha }, \end{equation}$$

(1)

where E_max is the maximum energy of the CR protons at t_Sedov and is set to 3 PeV, and α is the power-law index that is determined so that the current maximum energy (E_now) is equal to E_esc(t) at the current age t_age = 6.6 × 10³ yr (Leahy & Tian 2007). Here, we adopt E_now = 300 GeV, which is the same value as that of the electrons in the SNR shell at present (see also subsection 3.2). However, in general, the maximum energy of CRs in the SNR shell depends on the particle species and is determined by the balance between acceleration and escape from the shell or energy loss due to radiative cooling. Araya (2014), assuming that radiative cooling limits the maximum CR electron energy, estimated that the gyro factor would be too large, ∼660, compared to the standard SNR value. In such a case, as shown in Ohira et al. (2012), the escape process determines both the highest energies of CR electrons and protons. Therefore, according to Ohira et al. (2012), we can assume that the current maximum energy of electrons is the same as that of the protons. The timescale for a particle with the energy E, t_esc(E), to escape into interstellar space from the supernova explosion as a function of the particle energy E is given by

$$\begin{equation} t_{\rm esc}(E)=t_{\rm Sedov}\!\left(\frac{E}{E_{\max }}\right)^{-1/\alpha }. \end{equation}$$

(2)

The distribution function f_p,out per unit energy per unit volume of the CR protons escaping into interstellar space can be obtained by solving the transport equation,

$$\begin{equation} \frac{\partial f_{\rm p,out}}{\partial t}(t, r, E)-D_{\mathrm{ISM}}(E) \Delta f_{\rm p,out}(t, r, E)=q_{\rm p,s}(t, r, E), \end{equation}$$

(3)

where t is the age of the SNR, r is the distance from the SNR center, D_ISM is the diffusion coefficient in the ISM, and q_p,s is the injection rate of the CR protons from the SNR shock into interstellar space per unit energy, unit volume, and unit time. We adopt the following form of the diffusion coefficient D_ISM:

$$\begin{equation} D_{\mathrm{ISM}}(E)=D_0 \left(\frac{E}{10\,{\rm GeV}}\right)^{\delta }, \end{equation}$$

(4)

where D₀ is the diffusion coefficient of CRs at E = 10 GeV. We assume D₀ = 3 × 10²⁸ cm² s⁻¹ and the index δ = 1/3, the latter of which is consistent with the Galactic mean expected in the CR propagation model (e.g., Blasi & Amato 2012). Following Ohira, Murase, and Yamazaki (2011), we assume that CRs with the energy E are injected from the SNR shell R_esc = u_sh,0t_Sedov[t_esc(E)/t_Sedov]^2/5 at t = t_esc(E). We also assume that the energy spectrum of the injected CRs is monochromatic with the energy E and that particles start to escape at any given time [see equation (2)]. Taking account of these assumptions, the injection rate q_s is given by

$$\begin{equation} q_{\rm p,s}=\frac{N_{\mathrm{esc}}(E)}{4\pi r^2} \delta [r-R_{\rm esc}(E)] \delta \lbrace [t-t_{\rm esc}(E)]\rbrace , \end{equation}$$

(5)

where r is the displacement from the center of the SNR and N_esc is the spectrum of all the CRs that have escaped from the SNR up to the present time. In this paper, we assume that the total energy of the escaped CR protons is proportional to the explosion energy of the supernova explosion E_SN, namely

$$\begin{equation} \int _{E_{\rm now}}^{E_{\max }} E N_{\rm esc}(E) dE =\eta E_{\rm SN}, \end{equation}$$

(6)

where η is the acceleration efficiency coefficient. Here, we assume that N_esc is a power-law function of E with an index p_esc. According to Ohira et al. (2010), it is shown that the index p_esc of CRs escaping from the SNR is steeper than the index p_SNR of CRs confined in the SNR shell. While the value of p_esc is determined by the time evolution of the maximum energy of CRs at the SNR shell and the CR production rate, we treat p_esc as a parameter. By assuming a power-law form of N_esc, the integral of equation (6) yields

$$\begin{eqnarray} && N_{\rm esc}(E)= \nonumber\\ && \left\lbrace \begin{array}{ll}\displaystyle {\frac{\eta \left(2-p_{\rm esc}\right) E_{\rm SN}}{E_{\rm now}^2}} \left[\left(\frac{E_{\rm max}}{E_{\rm now}}\right)^{2-p_{\rm esc}}-1\right]^{-1}\left(\frac{E}{E_{\rm now}}\right)^{-p_{\rm esc}} & p_{\rm esc}\ne 2\\ \displaystyle {\frac{\eta E_{\rm SN}}{\ln \left(E_{\max }/E_{\rm now}\right)}E^{-2}} & p_{\rm esc}=2 \end{array} \right. . \end{eqnarray}$$

(7)

The propagation models (Obermeier et al. 2012; Yuan et al. 2017) expect the spectral index of particles with a Galactic CR origin to be ∼2.4, which is adopted as a fiducial parameter of p_esc in this modeling. The effect of this parameter on the modeling will be discussed in subsection 3.3. The energy spectrum of the escaped CRs is obtained by combining the transport equation (3) and equations (5) and (7). Specifically, the solution of equation (3) for E > E_esc(t) is given by, according to Ohira et al. (2010),

$$\begin{eqnarray} f_{\rm {p,out }}(t, r, E)&=& \frac{N\,_{\rm {esc }}(E)}{4 \pi ^{3 / 2} r R\,_{\rm {esc}} R\,_{\mathrm{d,p}}} \nonumber \\ &&\times \left\lbrace \exp\! \left[-\frac{\left(r-R\,_{\rm {esc}}\right)^{2}}{R\,_{\rm d,p}^{2}}\right] -\exp \left[-\frac{\left(r+R\,_{\rm {esc }}\right)^{2}}{R\,_{\rm d,p}^{2}}\right]\right\rbrace \nonumber \\ \end{eqnarray}$$

(8)

where R_d,p is the diffusion length of a CR proton with the energy E, defined as

$$\begin{equation} R_{\rm d,p}\equiv \sqrt{4D_{\mathrm{ISM}}(E)\left[t-t_{\rm esc}(E)\right]}. \end{equation}$$

(9)

In order to obtain the energy spectrum of CR electrons f_e,out, it is necessary to solve the following transport equation with radiative cooling:

$$\begin{eqnarray} \frac{\partial f_{\rm e,out}}{\partial t}(t, r, E) -D_{\mathrm{ISM}}(E) \Delta f_{\rm e,out}(t, r, E) \nonumber \\ \quad+\,\frac{\partial }{\partial E}\left[P(E)f_{\rm e,out}(t, r, E)\right] =q_{\rm e,s}(t, r, E), \end{eqnarray}$$

(10)

where P(E) is the energy loss rate of CR electrons, and q_e,s is the injection rate of CR electrons, which can be written as

$$\begin{equation} q_{\rm e,s}=K_{\rm ep}q_{\rm p,s}, \end{equation}$$

(11)

where K_ep is the ratio of the total energy of CR electrons to CR protons. As a cooling process, we consider only the synchrotron radiation, which is the most dominant effect for electrons with energies above |$\mathcal {O}({\rm GeV})$| (e.g., Berezinskii et al. 1990).

The temporal evolution of the maximum energy of CR electrons in the SNR shell is different from that of protons due to radiative cooling. After entering the Sedov phase, the maximum energy of CR electrons is first determined by the balance between radiative cooling and acceleration, and, after that, by the balance between acceleration and escape from the SNR shell, similar to the proton (Ohira et al. 2012). When these two phases switch, CR electrons start to escape from the SNR shell and be injected into the interstellar space. In this paper, we parametrize this switching time t_e as follows:

$$\begin{equation} t_{\rm e}=\xi _{\rm e} t_{\rm Sedov}, \end{equation}$$

(12)

where ξ_e is the time in units of t_Sedov at which the electron starts to escape and, in the limit of ξ_e → 1, the injection is identical to that of the proton [see equation (5)]. After t_e, since the time evolution of the maximum energy of electrons is expected to be the same as that of the proton (Ohira et al. 2012), we use equations (1) and (2). Originally, ξ_e is given based on the assumed environment (Ohira et al. 2012), but we set ξ_e = 1 to consider the limit where the gamma-ray flux of leptonic emissions is maximized in this modeling. The solution to equation (10) can be derived using the method described in the Appendix.

Considering only the synchrotron radiation as the cooling process [i.e., P(E) = Q_synE²; see equation (A7) of the Appendix for details], the solution for E_esc(t) < E < E_esc(t_e) can be written as:

$$\begin{eqnarray} f_{\rm {e,out}}(t, r, E) &=& \frac{K_{\rm ep}N_{\mathrm{esc}}\left(E_{c}\right)}{4 \pi ^{3 / 2} r R_{\mathrm{c}} R_{\rm d, e}} \frac{E_{c}^{2}}{E^{2}} \frac{1}{1-Q_{\rm syn}t_{c}E_{\rm c} / \alpha } \nonumber \\ && \times \left\lbrace \exp\!\! \left[-\frac{\left(r-R_{\mathrm{c}}\right)^{2}}{R_{\rm d, e}^{2}}\right] \!\!-\!\exp\!\! \left[-\frac{\left(r+R_{\mathrm{c}}\right)^{2}}{R_{\rm d, e}^{2}}\right]\right\rbrace . \nonumber \\ \end{eqnarray}$$

(13)

Once the energy spectrum of the CRs has been obtained, the flux of the gamma-ray emission is calculated based on the neutral pion-decay process in MCs for CR protons, and inverse Compton scattering and the relativistic bremsstrahlung for CR electrons. Here we assume that the MC is “optically” thin for the CRs and that the cloud consists of a spatially uniform gas with the total mass M_cl. For simplicity, we also assume that the MC is a sphere of radius R_cl. The spectrum of CRs in the MC can be written as follows:

$$\begin{equation} N_{\rm cl}(E)=\int _{d_{\mathrm{cl}}-R_{\mathrm{cl}}}^{d_{\mathrm{cl}}+R_{\mathrm{cl}}} \frac{\pi r}{d_{\mathrm{cl}}} \left[R_{\mathrm{cl}}^{2}-\left(r-d_{\mathrm{cl}}\right)^{2}\right] f_{s,\mathrm{out}}(r) d r, \end{equation}$$

(14)

where d_cl is the distance to the center of the MC from the center of the SNR, and s is species (⁠|$\rm e$| = electron, |$\rm p$| = proton). Note that, for the parameters used in this paper, a spatial gradient of the CR distribution is not important, and only a few % difference occurs even if |$N_{\rm cl} \sim 4\pi R_{\rm cl}^3f_{s,{\rm out}}(d_{\rm cl})/3$| is used.

The energy spectrum of CRs in the SNR shell, N_s,shell, is calculated consistently with the distribution function of escaped particles. By using the normalization equation (6) in the same way as equation (7), the spectra at E < E_esc of CR protons and electrons in the SNR shell can be written as

$$\begin{equation} N_{p,{\rm shell}}(E)=N\,_{\rm esc}(E_{\rm now})\left(\frac{E}{E_{\rm now}}\right)^{-p_{\rm SNR}} \end{equation}$$

(15)

and

$$\begin{equation} N\,_{{\rm e},{\rm shell}}(E)=K_{\rm ep}N\,_{\rm esc}(E_{\rm now})\left(\frac{E}{E_{\rm now}}\right)^{-p_{\rm SNR}}, \end{equation}$$

(16)

respectively, where p_SNR is the index of CRs that have not yet escaped from the SNR shell. We fix p_SNR = 2 in this paper because the index p_SNR of HB9 has been well determined by observations of the radio continuum (see section 1).

To calculate the spectra of non-thermal emissions, we use the radiative code naima (Zabalza 2015). As for the seed photon fields in the IC process, we assume the cosmic microwave background and Galactic far-infrared (FIR) radiation, the latter of which was not considered in the modeling by Araya (2014). The energy density of the FIR radiation is estimated to be 0.099 eV cm⁻³ at T = 27 K, using the package GALPROP (Vladimirov et al. 2011).

3.2 Application to the cloud regions around HB9

In order to calculate the gamma-ray flux at the cloud regions, R1 and R2, the distance between the MCs and SNR is required. By fitting the CO intensity map integrated over the velocity range between −10.4 and +2.6 km s⁻¹ with a two-dimensional symmetric Gaussian, we determine the center positions of clouds R1 and R2 to be (l, b) = (161|${_{.}^{\circ}}$|8 ± 0|${_{.}^{\circ}}$|1, 2|${_{.}^{\circ}}$|8 ± 0|${_{.}^{\circ}}$|1) and (162|${_{.}^{\circ}}$|6 ± 0|${_{.}^{\circ}}$|1, 1|${_{.}^{\circ}}$|6 ± 0|${_{.}^{\circ}}$|1), respectively, in the galactic coordinates. Then, the projected distances between the centers of the SNR and clouds R1 and R2 are calculated to be 17.8 and 39.4 pc, respectively, using the distance to the SNR from the Earth of 0.8 kpc (section 1 and table 2). We treat these projected distances as the actual distances in our discussion, although they should be considered as the lower limit of the true distances in reality due to the uncertainty of their locations along the line of sight. We also obtain the radii of the MCs from the standard deviation of the fitted Gaussian, which are 3.6 pc for R1 and 7.0 pc for R2.

Table 2.

Fiducial parameters used to calculate the model spectra.

SNR parameters	Symbol
SN explosion energy	ESN	0.3 × 1051 erg
Initial shock velocity	ush	109 cm s−1
Age of the SNR	tage	6.6 × 103 yr
Distance to the SNR		0.8 kpc
Acceleration efficiency*	η	0.1 (0.003)
Electron to proton flux ratio*	Kep	0.02 (1)
Maximum CR energy at t = tSedov	Emax	3 PeV
Current maximum energy of CRs	Enow	300 GeV
Magnetic field in the SNR	BSNR	|$8\,{\rm \mu G}$|
Particle index in the SNR	pSNR	2.0
Particle index after escaping from SNR	pesc	2.4
ISM parameters	Symbol
Number density	nISM	0.06 cm−3
Diffusion coefficient at E = 10 GeV	D0	3 × 1028 cm2 s−1
Index of dependence on E of diffusion	δ	1/3
Magnetic field in ISM	BISM	|$3\,{\rm \mu G}$|
Molecular cloud (MC) parameters	Symbol	R1	R2
Distance to the MC from SNR	dcl	17.8 pc	39.4 pc
Radius of the MC	Rcl	3.6 pc	7.0 pc
Average hydrogen number density*	nH	150 (1000) cm−3	200 (1000) cm−3
Mass of the MC*	Mcl	730 (4800) M⊙	7100 (36000) M⊙
Reflected SNR age	tref	−1.7 × 102 yr	−7.5 × 102 yr
Magnetic field in MC	BMC	|$3\,{\rm \mu G}$|

*

The values in parentheses are used in the leptonic-dominated model.

For the spectral modeling, the data points in the radio band for the HB9 SNR shell are obtained from the literature (Dwarakanath et al. 1982; Roger et al. 1999; Reich et al. 2003; Leahy & Tian 2007; Gao et al. 2011), while the radio flux of the mean local background including the cloud region at a frequency of 865 MHz (Reich et al. 2003) is used as an upper limit for the cloud regions. In the X-ray band, only the thermal emission from the hot gas inside the shell has been detected, while non-thermal emission has not been measured (e.g., Leahy & Aschenbach 1995; Sezer et al. 2019). We adopt the 0.1–2.5 keV flux of the thermal emission from HB9 (Leahy & Aschenbach 1995) as an upper limit for the non-thermal emission in the energy band.

We fit the spectra of R1 and R2 under two assumptions, the leptonic-dominated and hadronic-dominated models, and simultaneously reproduce the spectrum of the SNR shell with the leptonic emission. The electron to proton flux ratio, K_ep, is assumed to be 1 and 0.02 for the leptonic-dominated and hadronic-dominated models, respectively, the latter of which is consistent with the ratio in the local CR abundance (e.g., Aguilar et al. 2016).⁴ Figures 3 and 4 show the results of the leptonic-dominated and hadronic-dominated models, respectively. The total electron energy and the magnetic field in the SNR shell are in agreement with those estimated by Araya (2014), while the electron maximum energy at the shell (corresponding to E_now) is determined to be 300 GeV. Our obtained value of E_now is slightly lower than that derived in the previous study (Araya 2014) because the FIR radiation as a seed photon in the IC process is newly taken into account in this work. In the leptonic-dominated model, the gamma-ray emissions via a bremsstrahlung process of relativistic electrons dominate in the cloud regions, and thus the delayed gamma-ray spectra for 1–500 GeV have a hard index of ∼1.3, which contradicts the observed one at the cloud regions (figure 3). In the hadronic-dominated model (figure 4), the hadronic emission reproduces well the observed spectra even though the assumed parameters in the calculation are typical ones for an SNR and ISM (table 2). We note that the gamma-ray flux from the background CRs is lower than the observed data at higher energy than ∼1 GeV (see the middle and right panels of figure 4). By calculating the energy of the proton corresponding to the peak of photon spectra in each cloud region with the fiducial parameters, we estimate the epoch at which those protons escaped from the SNR shock by using equation (2); they are t_ref ≲ −1.7 × 10² yr for R1 and ≲ −7.5 × 10² yr for R2 dating back from now.

Fig. 3.

Broad-band spectral energy distributions of the non-thermal emission from the HB9 shell and the cloud regions with the leptonic-dominated model, using the parameters given in table 2. The left, middle, and right panels show the results for the HB9 shell and R1 and R2, respectively. The radio and X-ray data are taken from Dwarakanath, Shevgaonkar, and Sastry (1982), Reich, Zhang, and Fürst (2003), Leahy and Tian (2007), Roger et al. (1999), Gao et al. (2011), and Leahy and Aschenbach (1995), respectively (see text for details). The filled squares and open circles show the Fermi-LAT data derived in subsection 2.3 and Araya (2014), respectively. The lines represent each component of the emission models: synchrotron (blue), electron bremsstrahlung (magenta), IC (green), neutral pion decay (red), and the total gamma-ray emissions (black). In order to demonstrate how much of the Galactic diffuse background is present in the cloud regions, which may affect the systematic uncertainties in the Fermi analysis, we show the hadronic emission from the background CRs with the energy spectrum J_CR(E) = 2.2(E/GeV)^−2.75 cm⁻² s⁻¹ GeV⁻¹ sr⁻¹ (e.g., Dermer 1986) as the shaded gray region. (Color online)

Fig. 4.

Same as figure 3, but with the hadronic-dominated model. (Color online)

3.3 Implications for the diffusion coefficient

We investigate in the following procedure how the hadronic-dominated model curve varies depending on the input parameters and summarize the results in figure 5. First, we evaluate the dependence of the model curve on D₀, which is the value of the diffusion coefficient at E = 10 GeV (figures 5a and 5b). The Fermi-LAT spectra obtained in this work are found to be well reproduced with the Galactic mean of D₀ = 3 × 10²⁸ cm² s⁻¹. Orders of magnitude smaller D₀, in particular D₀ = 3 × 10²⁶ cm² s⁻¹, however, clearly fail to reproduce the observed spectra. In the previous studies for other SNRs (e.g., Fujita et al. 2009; MAGIC Collaboration 2020), the estimated values of D₀ were ∼10 times smaller than the Galactic mean, which were explained in conjunction with self-confinement caused by the generation of turbulent plasma waves (Wentzel 1974; Fujita et al. 2011; D’Angelo et al. 2018). Our diffusion coefficient value, which is close to the Galactic mean, indicates that the excitation of such turbulent plasma waves at the distances to R1 and R2 is inefficient or that the wave damping has a significant effect. Secondly, we find that the model with either δ = 1/3 or 1/2 is preferred over that with δ = 0 (figures 5c and 5d). Thirdly, we also find that E_max above 10 TeV still explains the data points (figures 5e and 5f). Given that there is a trend for a larger difference between the model curves in the higher-energy band, future observations in the TeV band will provide results more sensitive to determine the δ and E_max parameters.

Fig. 5.

Input parameter dependence of the model corresponding to regions R1 and R2. (a), (b) Dependence on D₀ (diffusion coefficient). Dotted, dashed, solid, and dot–dashed curves indicate the models with D₀ of 3 × 10²⁶, 3 × 10²⁷, 3 × 10²⁸, and 3 × 10²⁹ cm² s⁻¹, respectively. (c), (d) Dependence on δ (degree of dependence of D₀ on the energy). Dashed, solid, and dot–dashed curves indicate the models with δ of 0, 1/3, and 1/2, respectively. (e), (f) Dependence on E_max (maximum energy of the accelerated particles in the Sedov phase). Dashed, dot–dashed, and solid curves indicate the models with E_max of 1 × 10¹³, 1 × 10¹⁴, and 3 × 10¹⁵ eV, respectively. (g), (h) Dependence on δ with D₀ of 3 × 10²⁶ cm² s⁻¹. Solid, dot–dashed, dashed, and dotted curves indicate the models with δ of 1/3, 0.7, 1.1, and 1.5, respectively.

In order to verify that the Galactic mean value of the diffusion coefficient (D₀) is appropriate to explain the gamma-ray spectra, we also try to model the gamma-ray spectra, fixing D₀ at two orders of magnitude smaller D₀ = 3 × 10²⁶ cm² s⁻¹. As can be seen in figures 5a and 5b, for smaller D₀, the spectrum shifts to the higher-energy side. Also, as shown in figures 5e and 5f, the spectrum below |$\mathcal {O} \left(10^{11}\right)$| eV does not depend much on E_max. Therefore, the deviation of the spectrum in the lower-energy band due to the small D₀ should be countered by δ dependence. Once we assume δ > 1.0 while other parameters are kept the same as in table 2, the model curves are roughly consistent with the observed gamma-ray spectra (figures 5g, 5h). However, δ > 1.0 is inconsistent with the CR propagation model (e.g., Blasi & Amato 2012). Thus, a small diffusion coefficient (i.e., D₀ ∼ 10²⁶ cm² s⁻¹) is unlikely.

We have assumed p_esc = 2.4 in the modeling so far, but the discussion on the diffusion coefficient may be affected by this parameter as well as δ. Here, we attempt to fit the observed spectra by considering two cases where a flatter (p_esc = 2.0) or steeper (p_esc = 2.7) index is assumed. Note that these assumptions are not favored by the propagation model as mentioned in subsection 3.2. The results are shown in figures 6 and 7. Although the model parameters except for p_esc and η are fixed at the values in table 2, the observed spectra can be reproduced by giving a feasible value of η = 0.3 (0.05) for p_esc = 2.0 (2.7) (figure 6). Figure 7 shows the dependence on D₀ and δ for R2 as in figures 5b and 5h. To explain the observed spectrum with the diffusion coefficient two orders of magnitude smaller, D₀ = 3 × 10²⁶ cm² s⁻¹, than the Galactic mean, an unrealistic assumption that δ > 1.0 is required for both cases of p_esc = 2.0 and 2.7. This result is consistent with the case where p_esc = 2.4; i.e., the diffusion coefficient of the Galactic mean value is preferred regardless of the assumption on p_esc.

Fig. 6.

Energy spectra with the hadronic-dominated model under different assumptions on p_esc. The dashed, solid, and dotted lines represent the model emissions with p_esc of 2.0, 2.4, and 2.7, respectively. The model with p_esc = 2.4 is the same as in figure 4, while for the models with p_esc = 2.0 and 2.7, the acceleration efficiency is given as η = 0.3 and 0.05, respectively.

Fig. 7.

Input parameter dependence of the model in the R2 region when p_esc = 2.0 (left) and 2.7 (right) are assumed. (a), (b) Dependence on D₀ (diffusion coefficient). (c), (d) Dependence on δ fixing at D₀ = 3 × 10²⁶ cm² s⁻¹.

4 Conclusions

We have analyzed the GeV gamma-ray emissions in the vicinity of SNR HB9 with the Fermi-LAT data spanning 12 yr, aiming to quantify the evolution of DSA. We detected significant gamma-ray emission spatially coinciding with two MCs in the vicinity of the SNR. We found that the gamma-ray spectra above 1 GeV at the cloud regions could be characterized by a simple power-law function with indices of 1.84 ± 0.18 and 1.84 ± 0.14, which are flatter than that at the SNR shell of 2.55 ± 0.10. By modeling the diffusion of the CRs that escaped from SNR HB9, we found that the gamma-ray emission could be explained with the scenario that both CR protons and electrons accelerated by the SNR illuminate MCs rather than the leptonic scenario (see figures 3 and 4). This result implies that the maximum energy with DSA in younger SNRs is likely to be higher than that in older ones. We also found that the Galactic mean value of the diffusion coefficient (D₀) is appropriate to explain the observed gamma-ray spectra, indicating that self-confinement by turbulent plasma waves is not effective in the vicinity of SNR HB9. Future observations in the TeV band will provide more sensitive results to determine the maximum energy of the SNR in the past.

Acknowledgements

We would like to thank Toshihiro Fujii, Yutaka Fujita, Hidetoshi Kubo, Takaaki Tanaka, and Kenta Terauchi for helpful discussion and comments. We thank the participants and the organizers of the workshops with the identification number YITP-W-20-01 for their generous support and helpful comments. This work is supported by a Grant-in-Aid for JSPS Fellows Grant No. 20J21480 (TO).

Appendix. Derivation of the spectrum of CR electrons

The transport equation for CR electrons is

$$\begin{eqnarray} && \frac{\partial f_{{\rm e},\rm out}}{\partial t}(t, r, E) -D_{\mathrm{ISM}}(E) \Delta f_{{\rm e},\rm out}(t, r, E) \nonumber\\ &&\quad +\frac{\partial }{\partial E}\left[P(E)f_{{\rm e},\rm out}(t, r, E)\right] =q_{\rm e,s}(t, r, E). \end{eqnarray}$$

(A1)

According to Atoyan, Aharonian, and Völk (1995), the Green’s function of equation (A1) is

$$\begin{eqnarray} G(t, \mathbf {r}, E;t_0,\mathbf {r}_0;E_0) &=& \frac{P\!\left(E_{t-t_0}\right)}{\pi ^{3 / 2} P(E) r_{\mathrm{dif}}^{3}} \nonumber\\ && \times\,\exp \left(-\frac{\left|\mathbf {r}-\mathbf {r}_0\right|^{2}}{r_{\mathrm{dif}}^{2}}\right) \delta \left(E_{t-t_0}-E_0\right), \nonumber\\ \end{eqnarray}$$

(A2)

where |$E_{t-t_0}$| is defined by the following implicit relation:

$$\begin{equation} t-t_0=\int _{E}^{E_{t-t_0}} \frac{d E^{\prime }}{P(E^{\prime })}, \end{equation}$$

(A3)

 

$$\begin{equation} r_{\mathrm{dif}}=\sqrt{4\Delta u}, \end{equation}$$

(A4)

and

$$\begin{equation} \Delta u=\int _{E}^{E_{t-t_0}} \frac{D_{\rm ISM}(E^{\prime })}{P(E^{\prime })} d E^{\prime }. \end{equation}$$

(A5)

From linearity, the solution to equation (A1) can be calculated as follows:

$$\begin{eqnarray} f_{\rm {e, out}}(\mathrm{r}, t, E) &=& \int d E_0 \int d^{3} \mathrm{r}_{0} \int d t_{0} G(t, \mathbf {r}, E;t_0,\mathbf {r}_0;E_0) q_{\rm e,s} \nonumber\\ && \times\, \left(\mathrm{r}_{\mathbf {0}}, t_{0}, E_0\right). \end{eqnarray}$$

(A6)

In this paper, only synchrotron radiation is considered as the most important process in radiative cooling, so the cooling function P(E) can be written as follows:

$$\begin{equation} P(E)=\frac{4}{3}\sigma _{\rm T} c \left(\frac{E}{m_{\rm e} c^2}\right)^2 \frac{B_{\rm ISM}^2}{8\pi }\equiv Q_{\rm syn}E^2, \end{equation}$$

(A7)

where σ_T is the Thomson scattering cross-section, m_e is the electron mass, and B_ISM is the value of the magnetic field in the interstellar space. Evaluating equation (A3) using equation (A7), we obtain

$$\begin{equation} E_{t-t_0}=\frac{E}{1-Q_{\rm syn}E(t-t_0)}. \end{equation}$$

(A8)

For equation (A5), we can use equations (4) and (A7) as follows:

$$\begin{equation} \Delta u=\frac{r_{\rm dif}^2}{4} = \frac{D_{\rm ISM}(E)}{\left(1-\delta \right)Q_{\rm syn}E} \left[ 1-\left(\frac{E}{E_{t-t_0}}\right)^{1-\delta } \right]. \end{equation}$$

(A9)

By using equations (1), (5), (11), (A2), (A8), and (A9) and performing the integration of equation (A6), the solution to equation (A1) can be obtained as follows:

$$\begin{eqnarray} f_{\rm {e,out}}(t, r, E) &=& \frac{K_{\rm ep}N\,_{\mathrm{esc}}\left(E_{\rm c}\right)}{4 \pi ^{3 / 2} r R_{\mathrm{c}} R_{\rm d, e}} \frac{E_{\rm c}^{2}}{E^{2}} \frac{1}{1-Q_{\rm syn}t_{\rm c}E_{\rm c} / \alpha } \nonumber\\ && \times\, \left\lbrace \! \exp\! \left[\!-\frac{\left(r-R_{\mathrm{c}}\right)^{2}}{R_{\rm d, e}^{2}}\!\!\right]\! -\exp\! \left[\!\!-\frac{\left(r+R_{\mathrm{c}}\right)^{2}}{R_{\rm d, e}^{2}}\right]\right\rbrace , \nonumber\\ \end{eqnarray}$$

(A10)

where t_c = t_c(t, E) is the time when the electron whose energy is E at the current time t is injected and determined as a solution to an algebraic equation |$t_{\rm c}=t_{\rm esc}(E_{t-t_{\rm c}})$|⁠, or, equivalently, the below equation:

$$\begin{equation} \frac{E}{1-Q_{\rm syn}E\left(t-t_{\rm c}\right)}=E_{\rm esc}(t_{\rm c}); \end{equation}$$

(A11)

E_c(t, E) = E_esc(t_c), R_c(t, E) ≡ R_esc(E_c), and R_d,e(t, E) is the diffusion length of CR electrons:

$$\begin{equation} R_{\rm d, e}=\sqrt{\frac{4D_{\rm ISM}(E)}{\left(1-\delta \right)Q_{\rm syn}E} \left[ 1-\left(\frac{E}{E_{\rm c}}\right)^{1-\delta } \right]}. \end{equation}$$

(A12)

Here, we have assumed that there is only one solution to equation (A11); however, this is not generally true. The condition for there to be only one solution can be written as follows:

$$\begin{equation} t_{\rm Sedov} \ge t_{{\rm c}, \rm { crit }} \equiv t_{\rm {Sedov }} \left(\frac{Q_{\rm syn} t_{\rm {Sedov }} E_{\rm max}}{\alpha }\right)^{1 /(\alpha -1)} \end{equation}$$

(A13)

For the range of parameters that we use in this paper, this condition is satisfied; thus the assumption of only one solution is justified.

Footnotes

1

〈https://fermi.gsfc.nasa.gov/ssc/data/analysis/software〉.

2

〈https://fermi.gsfc.nasa.gov/ssc/data/Cicerone/〉.

3

Ergin et al. (2021) explored the possibility of a hadronic origin for gamma-rays from the HB9 shell using the 10 yr data of Fermi-LAT but nevertheless found no hint. However, since they performed the modeling study without radio and X-ray data, they could not restrict the electron distribution from synchrotron radiation and thus could not judge whether the hadronic or leptonic model is preferred.

4

Here K_ep is the flux ratio of CR protons and electrons in the sources (i.e., in the SNR shell), which does not necessarily correspond to the measurements in the solar system. In our model, for example, by integrating equation (13) over the whole space, we can obtain the energy spectrum N_e,esc(t, E) of the total escaped CR electrons, namely

$$\begin{eqnarray*} N\,_{\rm e,esc}(t,E)=\int 4\pi r^2 dr f_{\rm e,out}=\frac{K_{\rm ep}N\,_{\rm esc}(E_{\rm c})}{1-Q_{\mathrm{syn}} t_{c} E_{c} / \alpha }\frac{E_{\rm c}^2}{E^2}. \end{eqnarray*}$$

At the time t = t_age [i.e., E_esc(t) = E_now], the value of N_e,esc is N_e,esc(t_age, E_now) = K_epN_esc(E_now)/(1 − Q_synt_ageE_now/α), which is not identical to the value at the SNR shell, i.e., N_{e, shell}(t_age, E_now) = K_epN_esc(E_now) ≠ N_e,esc(t_age, E_now). This difference is due to the radiative cooling effect in the interstellar space; hence, the flux at the source does not necessarily match one of the escaped particles. On the other hand, the propagation effect should be taken into account if comparing the value of K_ep with the measurements in the solar system. However, the factor 1/(1 − Q_synt_ageE_now/α) is almost unity in the parameter range in this paper; thus here we compare K_ep directly to the observed value on the Earth.