© L. Ambrogi et al. 2019

1 Introduction

Supernova remnants (SNRs) are believed to be sources of Galactic cosmic rays (CRs). According to the SNR paradigm, CRs are accelerated through diffusive shock acceleration (DSA) at SNR shocks and propagate diffusively in the Galaxy (Bell 1978; Blandford & Ostriker 1978). A confirmation of this acceleration scenario can be provided by gamma-ray observations: accelerated CRs might undergo interactions with the matter field in the region of the source, leading to the production of hadronic gamma-ray photons through π⁰-decay or leptonic gamma-ray photons through bremsstrahlung processes. Relativistic CR electrons might also produce leptonic gamma rays through inverse Compton scattering off the photon field. The spectral energy distribution (SED) of the hadronic gamma rays is expected to show a peculiar feature, the so-called pion bump, which is an abrupt turnover of the spectrum that rises steeply below 200 MeV due to the finite rest mass of the neutral pions. At the time of writing, this spectral feature is claimed to be detected in three bright SNRs: IC 443 (Ackermann et al. 2013), W44 (Giuliani et al. 2011; Ackermann et al. 2013), and W51C (Jogler & Funk 2016). All of these three objects are middle-aged SNRs interacting with molecular clouds, characterized by a shell-like radio morphology combined with a center-filled X-ray morphology.

Other middle-aged SNRs have been successfully detected so far. The SNR HB 21 (also known as G89.0+4.7) belongs to this class of objects and is a few tens of thousands years in age (Koo & Heiles 1991; Leahy & Aschenbach 1996). Similar to IC 443, W44, and W51C, HB 21 is also a mixed morphology SNR. Its radio shell has an angular size of 120′ × 90′, slightly tilted in the northwest-southeast direction (Hill 1974). Thermal X-ray emission is concentrated in the interior of the radio shell and appears elongated in the same direction (Lazendic & Slane 2006; Pannuti et al. 2010). Most previous studies assumed 0.8 kpc as the distance to HB 21, postulating an association with the molecular clouds in the Cyg OB7 complex on the eastern side of the SNR (Tatematsu et al. 1990; Koo et al. 2001). Indeed, on the eastern boundary of the SNR, HI observations revealed a rather uniformly distributed emission (the so-called wall) as well as three large molecular clouds (clouds A, B and C hereafter) detected through CO observations (Tatematsu et al. 1990). However, no direct kinematic evidence of interaction between the eastern clouds and the SNR shell was observed. By comparing the HI and CO column density with the X-ray absorbing column density, the distance was later estimated to be ~1.7 kpc (Byun et al. 2006). We take this as the distance to HB 21 throughout this paper. Broad line emissions have been detected from shocked molecular clumps in the northern and southern regions of the SNR (clouds N and S; Koo et al. 2001), as well as in the northwest (cloud NW; Byun et al. 2006), providing the most solid evidence for the SNR interaction with molecular clouds. The interaction with northern and southern clouds is confirmed by observations in other wavelengths: the lack of thermal X-ray emission in correspondence to these shocked molecular clumps suggests that the shock of the SNR could have gone radiative owing to interactions with molecular clouds (Lazendic & Slane 2006). A recent study based on X-ray observations of the inner part of the SNR interpreted the existence of recombining thermal plasma as due to interactions with molecular clouds, which causes weaker magnetic turbulence, permits the escape of CR protons, and cools down the thermal electron (Suzuki et al. 2018). Moreover, near- and mid-infrared images of the same regions showed shock-cloud interaction features (Shinn et al. 2009, 2010).

In the second Fermi LAT source catalog, three point-like sources were initially associated with the SNR, interpreted as local maxima of the extended G89.0+4.7 region; these sources are 2FGL J2041.5+5003, 2FGL J2043.3+5105, and 2FGL J2046.0+4954 (Nolan et al. 2012). Later studies, also considered a fourth point-like object lying very close to the northeast edge of the SNR shell (2FGL J2051.8+5054) as part of the SNR and determined that the gamma-ray emission is produced by a single extended source (Reichardt et al. 2012; Pivato et al. 2013). Consequently, the SNR has been included in the third Fermi LAT catalog (labeled as 3FGL J2045.2+5026e) as one of the 12 SNRs firmly identified as spatially extended objects (Acero et al. 2015). According to roughly four years of Pass-7 data, its extended emission is well modeled as uniform disk with radius 1°.19 and centered at l = 88°.75, b = + 4°.65 and has a SED best described by a log parabola (Pivato et al. 2013).

In this paper we present an analysis based on more than nine years of the public Fermi LAT data from SNR HB 21 and we discuss the nature of its emission mechanism. In Sect. 2, we describe the data used for the study and their preparation.In Sect. 3, we report the morphological and spectrometric analysis performed on the data. In Sect. 4, we discuss the results and their interpretation. Finally, conclusions are summarized in Sect. 5.

2 Data analysis

The analysis presented is based on Pass-8 data collected from the beginning of the science operations of the instrument, on August 4 2008, to October 13 2017, corresponding to ~ 9.2 yr of data. The analysis of the data sample was performed using Fermipy (Wood et al. 2017), an open-source python package built on the Fermi Science Tools, which provides a high-level interface for analyzing LAT data. The LAT Science Tools exploited in this work is version v10r0p5, available at the Fermi Science Support Center¹. We selected class 128 events with reconstructed energy from 60 MeV to 870 GeV. We did not perform any selection in terms of point spread function (PSF) quality and energy reconstruction quality and we considered both front and back photon-to-pair conversion layer (evtype=3). We selected good time intervals requiring the instrument configuration to be in the science mode (LAT_CONFIG= =1) and removing data with bad quality flag (DATA_QUAL>0). Moreover, to filter out photons coming from the Earth’s limb, we rejected events with a zenith angle larger than 90°. We considered a region of interest (ROI) defined by a circle centered at the position of 3FGL J2045.2+5026e and with a radius of 10°. To compensate for the large PSF of the detector at the low energies, we account for the exposure of nearby objects by considering the contribution of gamma-ray emitters up to 20° from the ROI center.

We study the morphology and the spectral properties of HB 21 by means of a three-dimensional binned likelihood analysis: for the two spatial dimensions, the data are binned using square bins of 0°.1 side, while for the third dimension we assume ten logarithmically-spaced energy bins per decade. We verified the validity of the morphological results using PSF event types in a joint likelihood analysis. The initial model includes thecurrent Galactic diffuse emission model (gll_iem_v06.fits) and the corresponding model for the Extragalactic isotropic diffuse emission (iso_P8R2_SOURCE_V6_v06.txt), both of which are provided in the Science Tools², as well as the gamma-ray sources in the 3FGL catalog. We left the normalization of the Galactic and extragalactic diffuse models as a free parameter. The same was also done for 3FGL objects located within 5° from the source center, while for those at distances smaller than 3° we free all the spectral parameters. The resulting test statistic (TS) map reveals a number of significant excesses in the ROI, as shown in Fig. 1. We localized these gamma-ray emissions, which are not included in the 3FGL model, and estimated their extension (see Table 1). In a second iteration step the new hotspots have been added to the background model. We removed 3FGL J2045.2+5026e from such a model and we refer to it as the null hypothesis, as shown in Fig. 1.

Fig. 1 Test statistic map of the selected ROI with the new gamma-ray excesses found and modeled within this analysis. Red crosses indicate point-like sources, while red circles quantify the estimated extensions of the excesses (see Table 1). The position and the extension of the SNR HB 21 as defined in the 3FGL catalog is shown in green.

3 Results

3.1 Morphology

Several scenarios describing the morphology of HB 21 were considered and tested as follows:

• 1.

  a disk-type geometry with fixed radius of 1°.19, as estimated in Pivato et al. (2013);

• 2.

  a disk-type geometry with free radius;

• 3.

  a radial Gaussian shape with free width;

• 4.

  a ring-type geometry defined on the basis of the radio contours of the SNR at 4850 MHz (see Fig. 2b), i.e., an elliptical ring 40° tilted in the northwest-southeast direction with inner radii of 0°.3 and 0°.45 and outer radii of 0°.9 and 1°.05.

The morphology of the source was studied under a certain assumption on the shape of its spectral profile. We consider a log parabola energy distribution, which corresponds to the best description of the spectrum reported in the catalog (Acero et al. 2015; Pivato et al. 2013). The disk-type morphology, with the radius as free parameter, is the preferred model and has a significance of 44σ with respect to the null hypothesis.

This model predicts a radius smaller than that previously published (Pivato et al. 2013). Indeed, the TS map revealed the presence of a hotspot very close to HB 21, placed northeast with respect to HB 21 center and hereafter called hs_ne. The hotspot is well described as a point-like emitter located at l = 89°.72 ± 0°.02, b = + 4°.11 ± 0°.02 and with a power-law spectrum , with E₀ = 1000 MeV and estimated normalization N₀ = (0.18 ± 0.02) ⋅ 10⁻¹² MeV⁻¹ cm⁻² s⁻¹ and photon index α = 2.28 ± 0.07. The position of hs_ne is in good agreement with the gamma-ray source 2FGL J2051.8+5054 (l = 89°.77, b = 4°.15), as shown in Fig. 2. This 2FGL source was removed in the following catalogs because it was included as part of the extended SNR HB 21 (Reichardt et al. 2012; Pivato et al. 2013). Conversely, we found that the significance of the hypothesis of hs_ne to be separated from the SNR is ~9σ. Once the hotspot is included in the model as an independent source, the extension of HB 21 is eventually estimated to be 0°.83 ± 0°.02 and centered at l = 88°.76 ± 0°.02, b = +4°.74 ± 0°.02, as shown in Fig. 2.

We further investigate the morphology of the source by testing the hypothesis of two extended objects, one located in the northand one in the south with respect to the center of the source. This configuration is disfavored at 3σ level with respectto the single disk-type morphology, providing no indication for the existence of multiple sources within the region of SNR HB 21. In Fig. 2d, the residuals of the modeling with the single disk-type geometry are shown. An excess of gamma-ray emission is left in the southern part of the source, interestinglyclose to the region of the shocked CO clumps (Koo et al. 2001). To assess the significance of this excess, we postulate the existence of an additional component placed on top of the whole disk and located at l = 88°.47, b = 4°.45, which is the mean position of the three clumps S1, S2, S3. Its extension is estimated to be 0°.24 ± 0°.05 (red region in Fig. 2d) and has a power-law spectrum with parameters E₀ = 1000 MeV, N₀ = (0.76 ± 0.21) × 10⁻¹² MeV⁻¹ cm⁻² s⁻¹ and α = 2.32 ± 0.12. Only marginal evidence (4σ) of additional emission is found in correspondence to this cloud region, which does not allow us to establish firmly the presence of this additional component. These results, together with the fact that the gamma-ray source is produced by a single extended source, suggest that the bulk of the detected gamma-ray emission originates from the remnant and the contribution from the shocked clouds can be neglected.

3.2 Spectral analysis

We compute the SED of the source considering four independent logarithmic bins per decade of energy. In each energy bin, the spectral points are obtained performing a fit according to a power-law profile with spectral index α = 2. The corresponding SED of the SNR is shown in Fig. 3.

Our spectral analysis confirms that the spectrum shows clear evidence of curvature, as previously found in Reichardt et al. (2012) and Pivato et al. (2013), and provides for the first time a solid estimation of the statistical significance of the falloff. To account for the curved gamma-ray emission over the whole energy domain, we tested several parameterizations of the spectrum in this analysis: (i) a log parabola; (ii) a smooth broken power law; and (iii) a power law with an exponential cutoff. In Table 2 the results of the fit are reported. Since the spectral models considered are not nested, we use the Akaike information criterion as an estimator of the relative quality of the models: , with k is the number of estimated parameters in the model and is the maximum value of the likelihood function (Akaike 1974). The favored model, which corresponds to the minimum AIC value, is the log parabola, defined as , where E₀ = 1000 MeV and estimated parameters α = 2.35 ± 0.03, β = 0.39 ± 0.03 and N₀ = (0.157 ± 0.004) × 10⁻¹⁰ MeV⁻¹ cm⁻² s⁻¹. The analytical description of the spectrum according to the log parabola function makes it possible to reproduce locally the observed steepening of the gamma-ray energy distribution, which moves from α ~ 0.5 between 60 and 200 MeV, to α ~ 1 from 200 to 400 MeV, then to α ~ 2 from 400 to 1000 MeV and then to α ~ 3 above 1000 MeV. Globally, with respect to the catalog values (Pivato et al. 2013), our results show a spectrum with a slightly harder slopeand weaker flux; this difference is explained by the different extraction regions used for the estimation of the gamma-ray emission from the SNR.

In orderto assess the significance of the observed low energy falloff of the spectrum, we perform a likelihood fit of HB 21in the energy range between 60 MeV and 1 GeV assuming: (i) a simple power-law in the form , where E₀ = 200 MeV; and (ii) a smooth broken power law parametrized as , where E₀ = 200 MeV and β = 0.1. The improvement in the likelihood fit achieved with the smooth broken power-law profile is estimated through theTS of the two nested models. We prove that the break is detected with a significance of 5.4σ and its characteristic energy is estimated to be E_b = 389 ± 95 MeV, with a change in photon index from α₁ = 1.04 ± 0.28 to α₂ = 2.13 ± 0.20, and a normalization N₀ = (2.53 ± 0.25) × 10⁻¹⁰ MeV⁻¹ cm⁻² s⁻¹.

Fig. 2 Test statistic map of the selected ROI around the SNR HB 21. The map in panel a is computed in the whole energy interval considered for the analysis, i.e., above 60 MeV. Overlaid to this map, the 3FLG sources are reported; 3FGL J2045.2+5026e is also included and defined in the catalog as a disk-type geometry with a radius of 1°.19 (green markers). The best-fit model for the SNR HB 21 as found in this work is also shown (disk-type geometry with a radius of 0°.83), together with the nearby hotspot hs_ne (white regions). The position of the point source 2FGL J2051.8+5054 is presented as well (red cross). Panel b: TS map above 500 MeV is presented, together with the contours of radio continuum emission at 4850 MHz (white) and the ROSAT X-ray (dark blue) contours. The ROSAT data are available at https://skyview.gsfc.nasa.gov/current/cgi/titlepage.pl. Panel c: for energies above 1 GeV. On top of the TS map, the CO contours are reported in white. The CO distribution are taken from Dame et al. (2001) and integrated between −20 and 10 km s−1: six contour levels are linearly spaced from 1.5 to 28 K km s−1. We also show the position of the shocked CO clumps S1, S2, S3, N1, N2 (black circles; Koo et al. 2001) and clouds A, B, C (Tatematsu et al. 1990) and NW (Byun et al. 2006; white crosses). Panel d: TS map of the residual emission of the best modeling of the SNR and the hotspot (white regions) is shown above 1 GeV. In addition, the position of the shocked CO clumps is reported, together with the estimated cloud region (red region).

4 Discussion

4.1 Modeling of the nonthermal emission

The gamma-ray emission from remnants can originate from different mechanisms, induced either by proton interactions or by electron interactions. To investigate the nature of the parent particles, we fit the obtained gamma-ray SED of the SNR by both hadronic- and leptonic-induced gamma-ray spectrum. We use the Naima package (version v0.8.1) to compute the nonthermal radiation (Zabalza 2015). The gamma-ray emission resulting from pion-decay and bremsstrahlung processes are shown in Fig. 4. Different particle parent spectra are considered for the two scenarios. The Bayesian information criterion (BIC) enables us to select a model on a statistic basis; it is defined as , where n is the number of data points, k the number of parameters estimated by the model, and the maximum value of the likelihood function (Schwarz 1978). The model with the lowest BIC is preferred. In Tables 3 and 4 the BIC values are shown together with the estimated parameters of the parent particle spectra for pion-decay and bremsstrahlung mechanism, respectively. In the following we report a discussion for both the two emission processes.

Fig. 3 SNR HB 21 SED. The black dashed curve indicates the global best-fit model, i.e., the log parabola spectral shape. The red dashed curve up to 1 GeV is the broken power-law fit.

4.1.1 Pion decay

The hadronic models shown in Fig. 4a provide a good agreement with the gamma-ray data. In case of proton–proton interactions, the low energy break of the photon spectrum comes naturally from the kinematics of the process. However, a simple power-law parent spectrum fits the data only marginally, as it does not describe properly the high energy part of the spectrum. Consequently, a high energy steepening in the protons spectrum has to be introduced, either in the form of a cutoff or in the form of a break. The available statistics does not allow us to discriminate between a power law with an exponential cutoff or a broken power law. Under the hypothesis of hadronic-dominated scenario, this underlying proton spectrum has a spectral index on the order of α_1,p ~ 2.5. Similar properties are found also in the SNRs IC 443, W44 (Ackermann et al. 2013), and W51C (Jogler & Funk 2016), for which a proton index steeper thanα_p = 2 is established. Compared to the other three SNRs, the proton and the photon indexes are slightly steeper in HB 21, while the energy of the break in the gamma-ray spectrum is higher; this is still within statistical uncertainties however. The steepening of the spectrum can be explained in terms of confinement time of particles in the source, as the escape of the most energetic particles might lead to a quasi-stationary spectrum that is significantly softer than the acceleration spectrum. Therefore, depending on the age of the accelerator and on its distance and diffusion coefficient, different gamma-ray spectra are expected from source to source (Aharonian & Atoyan 1996).

In case of hadronic-dominated models, the CRs released by the SNR can escape the remnant and interact with the surrounding gas environment. Consequently, molecular clouds in the vicinity of the SNR can be illuminated by such escaping CRs and an enhanced gamma-ray emission is expected as due to the increased rate of π⁰ decays (Gabici et al. 2009). In this scenario, the gamma-ray brightness of hotspot hs_ne might be explained in terms of CRs escaping from HB 21 and interacting with the large molecular clouds along the eastern boundary of the remnant, which seem to be associated with the Cyg OB7 complex (Tatematsu et al. 1990). In particular, the northern cloud A represents the best candidate among the three eastern clouds, given the positional agreement with the hotspot. In Tatematsu et al. (1990), the density of cloud A is quoted to be of about 100 cm⁻³. However, the impossibility to resolve the extension of hotspot, together with uncertainties on its distance, prevent us from assessing its total mass and therefore the interaction of the SNR with the northeastern cloud cannot be firmly established, as also argued in Koo et al. (2001) and Byun et al. (2006) from CO observations.

Interestingly, in the southern part of the SNR, shocked molecular clumps with densities on the order of 10² –10⁴ cm⁻³ have been observed (Koo et al. 2001). This is so far the most convincing evidence for the interaction of the SNR with molecular clouds (Byun et al. 2006). As previously reported, we find a hint of a gamma-ray excess in correspondence to this region. Whether the shocked molecular clouds are inside the SNR or lie along the same line of sight and appear located in the interior of the shell as a projection effect, they can be the sites of efficient gamma-ray production through enhanced π⁰ decays. In the first scenario the gamma-ray emission could be explained in terms of trapped CRs (Uchiyama et al. 2010) and in the latter by escaping CRs (Aharonian & Atoyan 1996; Gabici et al. 2009). In the context of shocked clouds inside the SNR, preexisting CRs might be shock-accelerated by the overrun of the blast waves and the reacceleration in the dense clouds might result in a hardening of the spectrum of existing particles. So far, the spectral properties of the southern excess cannot be solidly constrained owing to limited statistics. Therefore, more data are needed to highlight possible spectral variations with respect to the remnant.

4.1.2 Bremsstrahlung

We hereafter consider models in which electron bremsstrahlung is dominant in the Fermi-LAT energy band. In order to reproduce the low energy drop of the gamma-ray emission, a break in the electron spectrum has to be postulated, as shown in Fig. 4b. If a cutoff at lower energies is also included in the electron spectrum (i.e., when E_e,min = 400 MeV), the Fermi data are satisfactorily well fitted in terms of bremsstrahlung emission. The shapes of the two components of the electron spectrum are estimated as α_1,e = 1.6 ± 0.3 and , with a break around 2 GeV. Such a spectrum notably evokes the spectrum of Galactic CR electrons in the interstellar medium (ISM), for which a low energy falloff with a break at ~ 2 GeV followed by a second component with slope ~3.1 are measured (Strong 2016; Strong et al. 2011). Remarkably, in this scenario the SNR would contribute to the acceleration of the diffuse component of the Galactic CR electrons.

In the leptonic scenario, we need to account for the radio emission of synchrotron origin. We assume the same parent spectral shape since in the easiest scenario the two emissions are produced by the same population of electrons, although the radio-shell morphology does not find a strict correspondence with the gamma-ray disk emission. In Fig. 5 the radio data (taken from Kothes et al. 2006; Pivato et al. 2013) are shown together with the synchrotron modeling performed with the Naima package (Zabalza 2015). Assuming a broken power-law spectrum and postulating E_e,min = 400 MeV, the electron parameters are estimated as α_1,e = 1.60 ± 0.17, , E_b,e = 6 ± 3 GeV, and . The energy of the break is slightly higher than that found for the bremsstrahlung. To reproduce the observed radio break with an electron break at E_b,e = 1.8 GeV (see Table 4 for the spectral parameters of the bremsstrahlung electrons), a magnetic field on the order of ~ 170 μG has to be postulated because of the scaling of the characteristic energy of the synchrotron radiation. However, within the uncertainties, the synchrotron and bremsstrahlung energy breaks, as well as the other parameters of the electron spectrum, are in agreement. The missing correspondence between the radio the gamma-ray emission morphologies could eventually be explained in terms of different cooling times: in case of bremsstrahlung interactions the cooling timescale depends on the ambient density as yr, while in case of synchrotron losses it scales as yr. The ratio of the two characteristic times cannot be significantly constrained owing to the uncertainties on the relevant quantities. It can vary from region to region of the SNR (e.g., due to the observed nonuniform gas distribution) and it is expected to be further amplified in the region of the shell because of the enhanced magnetic field. Higher resolution gamma-ray data might help to further investigate this scenario and potentially confirm the hypothesis of same production regions.

Fig. 4 Emission modeling of the gamma-ray radiation from HB 21. Panel a: pion-decay process. Panel b: bremsstrahlung mechanism.

4.2 Energetics

The total energy content in protons required to explain the Fermi-LAT emission is on the order of erg. A bremsstrahlung origin of the gamma rays implies for the parent electrons erg. The energy budget in protons is high, although we cannot exclude a higher local gas density, which would reduce the requiredbudget. For what concerns the electrons, from the modeling of the radio emission we obtained erg. The comparison of and predicts a proton density of ~10 cm⁻³, which is reasonable for a SNR. We note that as long as the simple one-zone model is assumed and the same electrons produce both the radio and the gamma-ray photons, then the product B² × n⁻¹ should be fixed. Consequently, for a magnetic field as high as 170 μG, a density above hundreds of protons per cubic-centimeter is expected, which is unlikely within a SNR. In this sense, more precise radio measurements of the SNR are needed to better contain the position of the break, and therefore the magnetic field. This would eventually allow us to rule out the hypothesis of same parent electrons, together with a deeper understanding of the SNR environment.

Fig. 5 Modeling of the radio emission from HB 21. Radio data points are from the literature: up to 4 × 10−5 eV are from Kothes et al. (2006), above are WMAP data taken from Pivato et al. (2013).

5 Conclusions

We analyzed more than nine years of Fermi-LAT Pass-8 data from the region of the SNR HB 21. The substantially richer statistics with respect to the previous analysis (based on roughly four years of Pass-7 data) led toa revision of the published results with relevant implications. The source is centered at l = 88°.76 ± 0°.02, b = + 4°.74 ± 0°.02 with an extension of 0°.83 ± 0°.02, that is, smaller than the extension previously published. The revised morphology allowed us to resolve a point-like source at the edge of the SNR, previously considered as part of HB 21. We proved that the bulk of the gamma-ray emission comes from the remnant and there is a hint of additional excess in correspondence to the region of the shocked molecular clouds observed at 4σ. The spectral analysis confirmed the presence of a break in the gamma-ray emission. We demonstrated the significant detection of the break at the level of 5.4σ and an energy of E_b = 389 ± 95 MeV.

The discovery of a spectral break at low energies is crucial for the understanding of the origin of the radiation from the SNRs. We showed that the SED of the SNR HB 21 can be fitted by both leptonic and hadronic gamma-ray production processes. A good agreement of the LAT data with the pion-decay mechanism is found, predicting a proton spectrum with a slope α ~2.5 and requiring a steepening at higher energies. The spectral and morphological properties of the high energy radiation from SNR HB 21 can be also interpreted as bremsstrahlung emission from a population of electrons accelerated within the SNR, whosespectrum resembles very closely the spectrum of CR electrons in the Galactic plane. This novel and remarkable result of our investigation shows that the SNR HB 21 contributes to the Galactic background electrons and, more in general, it can be treated as an indication of the major contribution of SNRs to relativistic electrons in the ISM. Within the uncertainties, the same population of electrons (with same spectral slopes and break) also explains theradio emission from the source. These results suggest that the break in the electron spectrum is not due to cooling processes, but rather reflects an intrinsic property of the source that is likely related to acceleration processes.

Acknowledgements

The authors would like to thank T. M. Dame and V. Zabalza for useful discussions. We thank the Max-Planck-Institut für Kernphysik as hosting institution. S.C. acknowledges the Polish Science Centre through grant agreement DEC-2017/27 / B / ST9 / 02272. R.Z. acknowledges the Alexander von Humboldt Foundation for financial support.

References

All Tables

All Figures

	Fig. 1 Test statistic map of the selected ROI with the new gamma-ray excesses found and modeled within this analysis. Red crosses indicate point-like sources, while red circles quantify the estimated extensions of the excesses (see Table 1). The position and the extension of the SNR HB 21 as defined in the 3FGL catalog is shown in green.
In the text

Fig. 2 Test statistic map of the selected ROI around the SNR HB 21. The map in panel a is computed in the whole energy interval considered for the analysis, i.e., above 60 MeV. Overlaid to this map, the 3FLG sources are reported; 3FGL J2045.2+5026e is also included and defined in the catalog as a disk-type geometry with a radius of 1°.19 (green markers). The best-fit model for the SNR HB 21 as found in this work is also shown (disk-type geometry with a radius of 0°.83), together with the nearby hotspot hs_ne (white regions). The position of the point source 2FGL J2051.8+5054 is presented as well (red cross). Panel b: TS map above 500 MeV is presented, together with the contours of radio continuum emission at 4850 MHz (white) and the ROSAT X-ray (dark blue) contours. The ROSAT data are available at https://skyview.gsfc.nasa.gov/current/cgi/titlepage.pl. Panel c: for energies above 1 GeV. On top of the TS map, the CO contours are reported in white. The CO distribution are taken from Dame et al. (2001) and integrated between −20 and 10 km s−1: six contour levels are linearly spaced from 1.5 to 28 K km s−1. We also show the position of the shocked CO clumps S1, S2, S3, N1, N2 (black circles; Koo et al. 2001) and clouds A, B, C (Tatematsu et al. 1990) and NW (Byun et al. 2006; white crosses). Panel d: TS map of the residual emission of the best modeling of the SNR and the hotspot (white regions) is shown above 1 GeV. In addition, the position of the shocked CO clumps is reported, together with the estimated cloud region (red region).

In the text

	Fig. 3 SNR HB 21 SED. The black dashed curve indicates the global best-fit model, i.e., the log parabola spectral shape. The red dashed curve up to 1 GeV is the broken power-law fit.
In the text

	Fig. 4 Emission modeling of the gamma-ray radiation from HB 21. Panel a: pion-decay process. Panel b: bremsstrahlung mechanism.
In the text

	Fig. 5 Modeling of the radio emission from HB 21. Radio data points are from the literature: up to 4 × 10−5 eV are from Kothes et al. (2006), above are WMAP data taken from Pivato et al. (2013).
In the text