Search for New Cosmic-Ray Acceleration Sites within the 4FGL Catalog Galactic Plane Sources

Fermi-LAT is a gamma-ray telescope that detects photons with energies from 20 MeV to more than 500 GeV by conversion into electron-positron pairs, as described in Atwood et al. (2009). The LAT is composed of three primary detector subsystems: a high-resolution converter/tracker (for measurement of the direction of the incident gamma-rays), a CsI(Tl) crystal calorimeter (for energy measurement), and an anticoincidence detector to identify the background of charged particles. Since the launch of the spacecraft in 2008 June, the LAT event-level analysis has been upgraded several times to take advantage of the increasing knowledge of how Fermi-LAT functions as well as the environment in which it operates. Following the Pass 7 data set, released in 2011 August, Pass 8 is the latest version of Fermi-LAT data. Its development is the result of a long-term effort aimed at a comprehensive revision of the entire event-level analysis and comes closer to realizing the full scientific potential of the LAT (Atwood et al. 2013). The current version of LAT data is Pass 8 P8R3 (Atwood et al. 2013; Bruel et al. 2018). It offers 20% more acceptance than P7REP (Bregeon et al. 2013). We used the SOURCE class event selection, with the instrument response functions (IRFs) P8R3_SOURCE_V3.

We used exactly the same data set as that used to derive the 4FGL catalog of sources, namely, 8 yr (2008 August 4–2016 August 2) of Pass 8 SOURCE class photons. This means that similarly to the 4FGL data set, our data were filtered removing time periods when the rocking angle was greater than 90° and intervals around solar flares and bright gamma-ray bursts (GRBs) were excised.

Pass 8 introduced a new partition of the events, called PSF event types, based on the quality of the angular reconstruction, with approximately equal effective area in each event type at all energies. Due to the very low signal-to-noise ratio at low energy, the angular resolution is critical to distinguish point sources from the background and we decided to use only PSF3 events (the best-quality events) below 100 MeV. We add PSF2 events between 100 MeV and 1 GeV. This high-energy bound was selected since middle-aged SNRs commonly exhibit a high-energy spectral break at around 1–10 GeV, which would then bias our low-energy analysis (Uchiyama et al. 2010). For both PSF3 and PSF2 events, we excised photons detected with zenith angles larger than 80° to limit the contamination from gamma-rays generated by cosmic-ray interactions in the upper layers of the atmosphere. That procedure eliminates the need for a specific Earth limb component in the model.

The data reduction and exposure calculations are performed using the LAT fermitools version 1.2.23 and family (Wood et al. 2017) v0.19.0. We used only binned likelihood analysis because the unbinned mode is much more CPU intensive and does not support energy dispersion.

We accounted for the effect of energy dispersion (reconstructed event energy not equal to the true energy of the incoming gamma-ray), which becomes significant at low energies (see below). To do so, we used edisp_bins = −3, which means that the energy dispersion correction operates on the spectra with three extra bins below and above the threshold of the analysis. ⁷⁸

Our binned analysis includes three logarithmically spaced energy bins between 50 and 100 MeV, and 9 energy bins between 100 MeV and 1 GeV. The Galactic diffuse emission was modeled by the standard file gll_iem_v07.fits and the residual background and extragalactic radiation were described by an isotropic component (depending on the point-spread function (PSF) event type) with the spectral shape in the tabulated model iso_P8R3_SOURCE_V3_PSF(3/2)_v1.txt. The models are available from the Fermi Science Support Center. ⁷⁹ In the following, we fit the normalizations of the Galactic diffuse and the isotropic components.

A crucial point that needs to be considered when analyzing LAT data at low energies is the effect of energy dispersion. For Pass 8, the energy resolution is < 10% between 1 and 100 GeV but it worsens below 1 GeV. It is ∼20% at 100 MeV and ∼28% at 30 MeV. The combination of energy dispersion and the rapidly changing effective area below 100 MeV could result in biased measurements of flux and spectral index of the source under study. In order to quantify the effects of energy dispersion, 200 simulations of the spectrum of IC 443 as published by Ackermann et al. (2013) were performed for an 8 yr observation time using the gtobssim tool included in the LAT fermitools. For these simulations, we assumed a point-source spatial model located at (R.A., decl. (J2000): 9451, 2266) and a smooth broken PL spectral model of the form:

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}={N}_{0}{(E/{E}_{0})}^{-{{\rm{\Gamma }}}_{1}}{\left(1+{\left(E/{E}_{\mathrm{br}}\right)}^{({{\rm{\Gamma }}}_{2}-{{\rm{\Gamma }}}_{1})/\alpha }\right)}^{-\alpha },\end{eqnarray} \tag{ 1 }$$

where α = 0.1, the break energy E_br = 245 MeV, and the spectral indices Γ₁ = 0.57, Γ₂ = 1.95. These simulations include the effect of energy dispersion. The analysis of these simulations was performed with the exact same configuration (region size, PSF components used, spatial bin size, energy interval) as the one used for real data. The only two parameters that have been varied in this study are the number of energy intervals and the value of the parameter edisp_bins as discussed in Section 2.2. For each combination of (energy bins, edisp_bins), we analyzed the 200 simulations, plotted the distributions of the reconstructed values of the break energy, Γ₁ and Γ₂, and fitted a Gaussian on each distribution.

The centroid of the Gaussian fit together with their size are reported in Figure 1 for the four tests performed: (12 energy bins, edisp_bins = 0), (12 energy bins, edisp_bins = −1), (12 energy bins, edisp_bins = −3), and (10 energy bins, edisp_bins = −3). As can be seen in this figure, the main effect is on Γ₁, as expected. If the energy dispersion is not taken into account (edisp_bins = 0), the spectrum falls less steeply at low energy and the spectral index Γ₁ is reconstructed with a value 0.1 higher than the simulated value set in the simulations. This is also true if the energy dispersion is taken into account with only one extra bin (edisp_bins = −1), which is not sufficient to properly take into account the effect of energy dispersion at these low energies even if this configuration has the advantage to reproduce slightly more accurately the value of the break energy. Even with a configuration using edisp_bins = −3, if the number of bins is too small, the reconstructed value of Γ₁ will be biased toward a lower value, which will artificially create a stronger break at low energy. This is directly due to the fact that the energy resolution varies with energy. This requires choosing an energy binning that is fine enough to capture this energy dependence. The best compromise that was found between good reconstruction and computation time (since higher values of edisp_bins or of the number of energy bins increase the CPU time) was obtained for a configuration using edisp_bins = −3 and 12 energy bins between 50 MeV and 1 GeV. This configuration was used for all results presented in the following.

Zoom In Zoom Out Reset image size

This analysis intends to find new cosmic-ray acceleration sites in our Galaxy. When cosmic-ray protons accelerated by a source penetrate high-density clouds, the gamma-ray emission is expected to be enhanced relative to the interstellar medium because of the more frequent proton–proton interactions. Targeting the presence of such clouds, we restricted our search to sources within 5° from the Galactic plane. In addition, we removed from our list all identified pulsars and active galactic nuclei (AGNs). For AGNs, we removed all subclasses, namely, flat-spectrum radio quasars, BL Lac-type objects, blazar candidates of uncertain type, radio galaxies, narrow-line Seyfert 1, steep spectrum radio quasars, Seyfert galaxies, or simply AGNs. Finally, to ensure that the source is significant in the low-energy domain covered by our analysis, we removed all sources with a significance below 3σ between 300 MeV and 1 GeV as reported in the 4FGL catalog. In the end, these selection criteria provide us with the list of 311 candidates reported in Appendix B.

We perform an independent analysis of the 311 candidates selected in Section 3.1. The procedure followed is inspired by the Fermi High-Latitude Extended Sources Catalog (Ackermann et al. 2018), which already used the functions provided by fermipy.

For each source of interest, we define a 20° × 20° region and include in our baseline model all 4FGL sources located in a 40° × 40° region centered on our source of interest (SOI). We model each 4FGL source using the same spectral parameterization as used in the 4FGL. For extended sources, we use the spatial models from the 4FGL and keep them as fixed parameters since the angular resolution between 50 MeV and 1 GeV does not allow us to perform a morphological analysis. Similarly, the positions of all point sources are fixed at their 4FGL values.

Starting from the baseline model, we proceed to optimize the model using the optimize function provided by fermipy. In this optimization step, we first fit the spectral parameters of the Galactic interstellar emission model and residual background together with the normalization of the five brightest sources.

Then, we individually fit the normalizations of all sources inside the region of interest (ROI) that were not included in the first step in the order of their total predicted counts in the model (Npred) down to Npred = 1. The optimization is concluded by individually fitting the index and normalization parameters of all sources with a test statistic (TS) value above 25 starting from the highest TS sources. This TS value is determined from the first two steps of the ROI optimization by $\mathrm{TS}=2(\mathrm{ln}{{ \mathcal L }}_{1}-\mathrm{ln}{{ \mathcal L }}_{0})$, where ${{ \mathcal L }}_{0}$ and ${{ \mathcal L }}_{1}$ are the likelihoods of the background (null hypothesis) and the hypothesis being tested (source plus background). This optimization is followed by a second one where the number of bright sources fit together with the diffuse backgrounds is increased to 10. This allows a better convergence for complex regions containing a large number of bright sources.

After optimizing the parameters of the baseline model components, we then perform a fit of the region by leaving free the normalization of all sources within 2° of the SOI, their spectral shape if their TS value is above 16, the normalization of all sources with TS > 100 in the ROI and the spectral shape of all sources in the ROI with TS > 200. If the number of degrees of freedom (Ndof) is above 100, we increase the two last TS criteria by 100 until the Ndof becomes smaller than 100.

Once this complete fit of the ROI is performed, we further refine the model by identifying and adding new point-source candidates. We identify candidates by generating a TS map for a point source that has a PL spectrum with an index of Γ = 2. When generating the TS map, we fix the parameters of the background sources and fit only the amplitude of the test source. We add a source at every peak in the TS map with TS > 16 that is at least 2° from a peak with higher TS due to the poor angular resolution at these low energies. New source candidates are modeled with a PL whose normalization and index are fit in this procedure. We then generate a new TS map after adding the point sources to the model and repeat the procedure until no candidates with TS > 16 are found. Though we do not expect to find a large number of additional sources, this step is crucial since we are not using the weights, first introduced for the 4FGL Catalog (Abdollahi et al. 2020). Indeed, the generation of the candidate list for the catalog is done above 100 MeV (instead of 50 MeV for our analysis) and each candidate is kept if the TS value obtained via a weighted maximum likelihood fit is above 25. These weights mitigate the effect of systematic errors due to our imperfect knowledge of the Galactic diffuse emission. As a consequence, the TS value of soft sources decreases.

In the final pass of the analysis, a second general fit of the ROI is performed using the same criteria as above to free the spectral parameters of all sources. If sources added previously by using the TS map fall below TS > 16, they are removed from the model. If their TS value is above 16 and they are located at a distance smaller than 5° from our SOI, we test iteratively for each of them the improvement of the log-normal representation with respect to the PL model. The LP model is defined as

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}={N}_{0}{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{-\left(\alpha +\beta \mathrm{log}\left(E/{E}_{0}\right)\right)},\end{eqnarray} \tag{ 2 }$$

where N₀ is the overall normalization factor to scale the observed brightness of a source, E₀ is a fixed scale energy (kept at 300 MeV in our analysis), and α, β are left free, which adds 1 degree of freedom with respect to the PL representation. The improvement of the LP model with respect to the PL one is performed by determining ${\mathrm{TS}}_{\mathrm{LP}}=2(\mathrm{ln}{{ \mathcal L }}_{\mathrm{LP}}-\mathrm{ln}{{ \mathcal L }}_{\mathrm{PL}})$. If TS_LP is above 9 (which corresponds to a 3σ improvement for one additional degree of freedom), we switch to the log-normal representation. The spectral parameters of all added sources located within 5° of a candidate are reported in Table 1. As can be seen, the curvature index β is hard to constrain for the additional faint sources, even if the LP model significantly improves the fit. It is also clear that several added sources are located within the Vela and Cygnus regions for which the morphological templates used for the Vela-X PWN or the Cygnus Cocoon are not precise enough to properly characterize the region. Because the 4FGL Catalog rejects most point sources found inside extended sources, this leaves many residuals that translate into sources.

Once the ROI is well characterized, we first test the TS value of our SOI in our energy range (50 MeV−1 GeV). If it is below 25, we stop the analysis for this source since it is not significantly detected in our pipeline. It is the case for 64 sources among the 311 selected and their TS values are reported in Table B1. If the TS of the SOI is above 25, we move on and we fix all sources located more than 5° from the SOI and we test the spectral curvature of our SOI.

To ensure that the curvature is real and affects several energy bins as would be the case for a pion-decay bump signature, we first test a log-normal representation for the source as defined in Equation (2). If TS_LP is below 9, we consider that no significant curvature is detected by our pipeline, we report this value in Table B1 and we stop the analysis of this source. This is the case for 167 sources in our sample.

If the value is above 9, we then test a smoothly broken PL following Equation (1), where N₀ is the differential flux at E₀ = 300 MeV and α = 0.1, as was done previously for the cases of IC 443 and W44 (Ackermann et al. 2013). This adds two additional degrees of freedom with respect to the PL model (the break energy E_br and a second spectral index Γ₂). The improvement with respect to the PL one is determined by ${\mathrm{TS}}_{\mathrm{SBPL}}=2(\mathrm{ln}{{ \mathcal L }}_{\mathrm{SBPL}}-\mathrm{ln}{{ \mathcal L }}_{\mathrm{PL}})$. Since this test requires the addition of 2 degrees of freedom to the fit and diffusive shock acceleration predicts Γ₂ ∼ 2, we also test the improvement of the smooth broken PL with the second index fixed at 2 with respect to the PL one ${\mathrm{TS}}_{\mathrm{SBPL}2}=2(\mathrm{ln}{{ \mathcal L }}_{\mathrm{SBPL}2}-\mathrm{ln}{{ \mathcal L }}_{\mathrm{PL}})$. We require TS_SBPL > 12 or TS_SBPL2 > 9 (implying a 3σ improvement for two and one additional degrees of freedom, respectively) to keep the source in the significant energy break list reported in Table 2. We switch to the SBPL parameterization for all sources detected in this list. This means that, when a source located within 5° shows a significant energy break, we re-optimize the ROI and we redo the whole process as illustrated in the flowchart in Figure 2. This procedure allowed the detection of 77 sources presenting a significant energy break in their low-energy spectrum. The values of TS_LP, TS_SBPL, and TS_SBPL2 for each of them are reported in Table B1.

Zoom In Zoom Out Reset image size

The primary source of systematic error in this low-energy analysis is the Galactic interstellar emission model (IEM). Our nominal Galactic IEM is the recommended one for PASS 8 source analysis. It represents the first major update to the LAT Collaboration's IEM since the model for the 3FGL catalog analysis gll_iem_v05.fits, developed for Pass 7 Source class, and later rescaled for Pass 8 Source as gll_iem_v06.fits. The development of the new model is described in more detail (including illustrations of the templates and residuals) online. ⁸⁰ The new model has higher resolution and correspondingly greater contrast but some shortcomings in the new Galactic IEM have been recognized when producing the 4FGL catalog.

To quantify the impact of diffuse systematics, we repeated our analysis for the 77 sources listed in Table 2 using the old diffuse rescaled for Pass 8 Source gll_iem_v06.fits. This alternative analysis means that the whole flowchart in Figure 2 was performed again, from the optimization of the ROI to the source finding algorithm, up to the determination of the spectral curvature. Performing the same complete analysis with the eight alternate diffuse models from Acero et al. (2016) would have become extremely CPU time-consuming and this is why the old diffuse model is only tested. Here, since we already know that the source presents a break with the new model, we directly tested if TS_SBPL > 12 or TS_SBPL2 > 9. If it was not the case, then this source was discarded from the final list of sources presenting significant energy break.

The second source of systematic error in our analysis is the instrument response functions (IRFs) and especially the inaccuracies in the effective area. Following the standard method (Ackermann et al. 2012), we estimated the systematic error associated with the effective area by calculating uncertainties in the IRFs, which symmetrically bracket the standard effective area and flip from one extremum to the other at the measured value of the break energy. Here we started from the best-fit model obtained with the standard IRF, which is optimized before running the final spectral fit of each candidate with each of the two bracketing IRFs. The source finding algorithm was not relaunched in this case since these changes mainly affect the spectral parameters of the source and will not produce extra sources in the field of view.

A third source of systematic error that can affect the presence or absence of a spectral break for the SOI is related to the inaccuracy of the emission models of nearby point sources. A thorough investigation of this effect is beyond the scope of this paper but we included in Table 3, for each candidate, the distance of the nearest source as well as the relative contribution of photons from the neighboring sources and the diffuse backgrounds. Those values show that the diffuse background impacts the sources much more than their neighbors, with the exception of 4FGL J2021.0+4031e around the bright PSR J2021+4026.

Overall, 56 sources among the 77 sources detected with the standard IEM and IRFs are confirmed by our systematic studies. The 21 candidates rejected are all sources that do not meet the TS_SBPL or TS_SBPL2 criteria when using the old diffuse model, while the inaccuracy in the effective area has a minor effect in our analysis as can be seen in Table 2. The spectral parameters of the confirmed sources are listed in Table 4. As can be seen in this table, even if the old diffuse background detects a significant energy break, the energy of this break can be significantly different than with the standard IEM, leading to large systematics as well on Γ₁. However, the value of Γ₂ is much more robust.

In addition to performing a spectral fit over the entire energy range, we computed an SED by fitting the flux of the source independently in 10 energy bins spaced uniformly in log space from 50 MeV–1 GeV. During this fit, we fixed the spectral index of the source at 2 as well as the model of background sources to the best fit obtained in the whole energy range except the normalizations of the Galactic diffuse and isotropic backgrounds. We determined the flux in an energy bin when TS ≥ 1 and otherwise computed a 95% confidence level Bayesian flux upper limit, assuming a uniform prior on flux following Helene (1983). The systematic studies with the old diffuse and bracketing IRFs were also computed on all spectral energy distribution (SED) points for the 56 confirmed sources and the two uncertainties were added in quadrature. When an upper limit was derived, the maximal and minimal upper limits derived in this energy interval are plotted to indicate the systematics related to this data point.

We detected 56 4FGL gamma-ray sources showing a significant energy break in their spectrum between 50 MeV and 1 GeV confirmed by our studies of systematics. As can be seen in Figure 3, the distribution of sources showing a significant break in their low-energy spectrum is more uniform in both latitude and longitude than the parent distribution even if there remains a peak at latitude 0 and in the Galactic Ridge.

Zoom In Zoom Out Reset image size

The sources that we detect significantly with our analysis (TS > 25) follow the same trend except for the region at ∼ 300° longitude, which contains more faint sources than the other regions of the plane. Figure 4 clearly shows that the sources that we do not detect with TS > 25 in our pipeline have predominantly low significance in the 4FGL catalog in the 300 MeV–1 GeV energy band, which is reassuring. However, there is no correlation between the significance value in the 4FGL catalog and the detection of a break with our pipeline. It can be seen in this same Figure since the distribution for the sources presenting a significant break is uniform.

Zoom In Zoom Out Reset image size

The association summary is given in Table 5 and is illustrated by the pie charts in Figure 5. Out of 311 candidates, 210 are unidentified, representing 67.5% of the sources analyzed. It is striking to see that only 26 UNIDs show a spectral break confirmed by our systematic studies (which represents 46.4% of the sources with significant breaks). The 30 remaining candidates out of 56 confirmed cases present an association reported by the 4FGL Catalog listed in Table 6.

Zoom In Zoom Out Reset image size

On the other hand, the fraction of sources associated with SNRs increases from 7.4% (23 out of 311 sources) to 23.2% (13 out of 56 sources). This makes SNRs the dominant class of sources with significant low-energy spectral break. Similarly, the fraction of sources associated with binaries increases from 1.6% (five out of 311) to 7.1% (four out of 56), showing that almost all binaries except 4FGL J1826.2-1450 (also known as LS 5039), show a significant spectral break. Despite their small fractions, binaries could contribute significantly to our population of sources with low-energy spectral breaks; however, it should be noted here that the spectral analysis is performed over 8 yr and these sources often present variable gamma-ray emission. A more thorough analysis of these sources would need to be done. Finally, only one SFR is analyzed (and confirmed) which prevents us from drawing a firm conclusion on this source class.

Figure 6 illustrates the distribution over the sky of the 56 4FGL gamma-ray sources showing a significant energy SFR break. The lack of these sources at latitudes smaller than −2° appears clearly. One can also note a large fraction of unidentified sources at longitude comprised between −50° and 50°. These sources are part of the large fraction of 4FGL unassociated sources located less than 10° away from the Galactic plane with a wide latitude extension hard to reconcile with those of known classes of Galactic gamma-ray sources.

Zoom In Zoom Out Reset image size

Looking now at the spectral parameters of the 56 confirmed sources, the distribution of the energy of the breaks detected by our analysis is relatively uniform between 70 and 700 MeV, with no breaks detected below and above this energy interval (as a direct consequence of the energy interval analyzed here) and a higher proportion of breaks at ∼400 MeV as illustrated by Figure 7. Interestingly, no low-energy spectral breaks (< 140 MeV) are detected for the 13 sources associated with SNRs. As can be seen on the top panel of this figure, the large error bars on this parameter prevent us from drawing any firm conclusion or even rejecting any candidate by a comparison with the standard value expected for proton–proton interaction indicated by the green line. On the other hand, there is a trend concerning the distributions of Γ₁ with a peak at ∼0.2 and ∼1.0. The peak at 0.2 is expected by proton–proton interaction (as indicated by the green line presenting the results of the simulations carried out in Appendix A) but the peak at 1 is not predicted, though it is present for a large number of SNRs interacting with MCs. It might be due to some confusion by the Galactic and isotropic diffuse background. A double-peaked distribution is also visible in Figure 8 for Γ₂ at ∼2.1 and ∼3.6. For this parameter, the distribution restricted to SNRs contains a single peak at ∼2.1. Looking now at the distribution of Γ₂ − Γ₁ in Figure 8 (right), a peak at ∼0.9 is highly pronounced for SNRs. This tends to show that the values obtained on Γ₂ and Γ₂ − Γ₁ could be used in the future to probe the type of particles radiating in a gamma-ray source.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The most famous sources with pion bump signature are the middle-aged remnants IC 443 and W44. Figure 9 presents the residual TS maps of the region of IC 443 (4FGL J0617.2+2234e) and W44 (4FGL J1855.9+0121e) as well as their SEDs, showing the overall agreement with the 4FGL SED points superimposed. This figure also illustrates the advantages of using a restricted energy range and different spectral shape than the 4FGL to better reproduce the significant energy break at low energy since we are not dominated here by photons at high energies. The spectral parameters reported in Table 4 for these two sources are in reasonable agreement with those published by Ackermann et al. (2013), knowing that this first analysis did not take into account the effect of energy dispersion and no systematic uncertainties were evaluated at that time.

Zoom In Zoom Out Reset image size

Among the 56 sources with significant breaks, one can see from the 4FGL Classification column listed in Table B1 that 10 sources are firm SNR identifications and three are associated with SNRs. Among the three SNR associations, 4FGL J1911.0+0905 (Figure 17 top right) is associated with W49B and thus can be safely identified as an SNR since it is one of the few other sources for which a pion-decay bump signature was published with W51C (4FGL J1923.2+1408e, Figure 17, middle left) and HB 21 (4FGL J2045.2+5026e, Figure 18, middle right). The only missing source for which a low-energy break has been published is Cassiopeia A (4FGL J2323.4+5849) but the break energy reported by Yuan et al. (2013) is at ${1.72}_{-0.89}^{+1.35}$ GeV, which seems consistent with our non-detection in the 50 MeV–1 GeV energy interval. The five sources confirmed by our analysis are all SNRs interacting with MCs. These MCs are excellent targets for cosmic-ray interactions and subsequent pion decay.

The hadronic scenario was also preferred for other LAT-detected SNRs interacting with MCs, though their gamma-ray analysis starting above a few hundred megaelectronvolts did not allow rejection of a leptonic scenario: the SNR HB 3 and the W3 H ii complex (Katagiri et al. 2016a), S147 (Katsuta et al. 2012), HB 9 (Araya 2014), the SNR G326.3-1.8 (Devin et al. 2018) and the SNR W28 (Hanabata et al. 2014). Our low-energy analysis presents a rapid turnover of the spectrum at low energy, which confirms the conclusions of the previous publications for 4FGL J0222.4+6156e (W3, see Figure 10, top left), 4FGL J0500.3+4639e (HB 9, see Figure 10, bottom right), 4FGL J0540.3+2756e (S147, Figure 11, top left) 4FGL J1552.9-5607e (G326.3-1.8, Figure 14, middle left) and 4FGL J1801.3-2326e (W28, see Figure 15, middle left). No significant curvature is detected for the SNR HB 3 but it should be noted that its gamma-ray emission is much fainter than the adjacent MC W3 (TS value of 75.9 with respect to 1307.1 for W3) and more data would be needed to constrain the low-energy spectrum of the SNR. A hadronic scenario was also invoked for the SNR Monoceros Loop (Katagiri et al. 2016b). In this case, the brightest gamma-ray peak is spatially correlated with the Rosette Nebula, a young stellar cluster and MC complex located at the edge of the southern shell of the SNR which has a role similar to W3 for the HB 3/W3 complex. The interaction between the SNR and the MC provides the target to naturally produce gamma-rays via proton–proton interaction and it is not a surprise that we confirm a spectral break at low energy for the Monoceros SNR (4FGL J0639.4+0655e, see Figure 11, bottom left) and for the Rosette complex (4FGL J0634.2+0436e, see Figure 11, middle right). More recently, modeling of the nonthermal emission of the gamma Cygni SNR (Fraija & Araya 2016; Fleischhack 2019), associated with the source 4FGL J2021.0+4031e, also suggested that the gamma-ray emission (analyzed above 100 MeV) might be of hadronic nature with enhanced gigaelectronvolt emission spatially coincident with the teraelectronvolt source VER J2019+407. Here again, our low-energy analysis detects a low-energy break in the spectrum of this SNR (see Figure 17, bottom right) but it should be noted that the bright gamma-ray emission from the pulsar PSR J2021+4026, lying near the center of the remnant, is very difficult to disentangle from the signal of the SNR at these low energies, which could lead to some contamination in the SNR spectrum. A follow-up study in the off pulse of the pulsar would therefore be needed to confirm the results obtained with our pipeline. This applies not only to supernova remnants but also to all sources coincident with (or very close to) a bright gamma-ray pulsar. It is even more clear for 4FGL J1514.2-5909e associated with the pulsar wind nebula MSH 15-52 and coincident with the soft gamma-ray pulsar PSR B1509-58. The very high low-energy flux visible in Figure 13 (bottom right) is most likely to the associated pulsar PSR B1509-58, which is not included in the 4FGL Catalog and would be hard to disentangle from the PWNe at these energies.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As discussed in Section 4.2, gamma-ray observations are suggestive of hadron acceleration in a number of SNRs: the young SNRs Tycho and Cassiopeia A, and the middle-aged remnants with pion-decay signature cited above. However, definite proof of proton acceleration, especially at petaelectronvolt energies, is still missing and alternative Galactic sources of cosmic rays could play a significant role.

The shocks generated by the stellar winds of massive stars or star-forming regions are among these cosmic-ray accelerators. In this respect, the detection of gamma-rays of the Cygnus region by the LAT (Ackermann et al. 2011) opened new perspectives by revealing the presence of a cocoon of freshly accelerated cosmic rays over a scale of ∼50 pc. Our analysis revealed a spectral break for the SFR analyzed, 4FGL J2028.6+4110e (see Figure 18, top left), which is associated with the cocoon. A very hard index Γ₁ = 1.00 ± 0.02_stat ± 0.37_syst is detected up to break energy at 383 ± 13_stat ± 138_syst MeV, followed by a spectral index Γ₂ = 2.23 ± 0.06_stat ± 0.24_syst, similar to those observed for the population of identified SNRs as can be seen in Figure 8. Complete modeling of the source at gamma-ray energies is beyond the scope of this paper but our results tend to favor the hadronic scenario, thus reinforcing the long-standing hypothesis that massive SFRs house particle accelerators.

Gamma-ray binaries, microquasars, and colliding wind binaries could also contribute to the sea of Galactic cosmic rays and at least contribute significantly to the population of sources with significant breaks as reported in Section 4.1. Spectral breaks have been detected for these three types of sources with 4FGL J0240.5+6113 associated with the high-mass gamma-ray binary (HMB) LS I+61 303 (Figure 10, top right), the HMB 4FGL J1018.9−5856 (Figure 12, bottom left), 4FGL J1045.1-5940 associated with the colliding wind binary η Carinae (Figure 12, bottom right) and 4FGL J2032.6+4053 associated with the microquasar Cyg X-3 (Figure 18, top right). However, this last source presents the highest value of spectral index Γ₁ (1.90 ± 0.16_stat ± 0.07_syst) among the 56 candidates, which does not really look like the standard pion bump signature observed for interacting SNRs. Finally, the source 4FGL J1405.1-6119 was recently identified as a high-mass gamma-ray binary using Fermi-LAT observations (Corbet et al. 2019), and should therefore be added to the small set of gamma-ray binaries detected in our analysis. Since significant variability was detected by the LAT for these five gamma-ray sources, an individual analysis taking into account their orbital period would be needed to see if the spectral break detected is a signature of proton–proton acceleration.

Among the sources for which a significant spectral break is detected with our pipeline, several are classified as SPP, UNK, or even unassociated as can be seen in Table 5. Among these three source classes, SPP is the only one for which the fraction of sources with significant break is similar to the analyzed fraction (11.9% versus 10.7%), while UNK and UNIDs both show a clear decrease between the analyzed fraction and the confirmed one (see Figure 5). The SPP are sources of unknown nature but overlapping with known SNRs or PWNe and thus candidates to these classes, while UNK are sources associated with counterparts of unknown nature. Unassociated, SPP and UNK represent 29.7% of the 4FGL sources: revealing the mystery of the nature of these unidentified gamma-ray sources might shed new light on the problem of the origin of Galactic cosmic rays. In this respect, three sources detected by our pipeline are of special interest since they are coincident with SNRs and/or dense molecular clouds.

This is the case for 4FGL J1601.3-5224 (Figure 14, middle right) coincident with the SNR G329.7+00.4, which presents a diffuse shell at radio energies (Whiteoak & Green 1996) but is not detected at any other wavelength. Our analysis indicates a soft spectrum Γ₂ = 3.78 ± 0.32_stat ± 0.89_syst with large systematics due to the diffuse background. The same systematics affect the value of the energy break showing that our results may suffer from contamination.

Similarly, the source 4FGL J1934.3+1859 (Figure 17, bottom left) is coincident with SNR G054.4-00.3 detected as a nearly circular shape and angular diameter of ∼40ʹ at radio energies (Junkes et al. 1992) while Swift and Suzaku X-ray observations allowed for the detection of the X-ray counterpart (Karpova et al. 2017) of the gamma-ray pulsar PSR J1932+1916 (Pletsch et al. 2013) located near the edge of the SNR. Suzaku observations also revealed diffuse emission with an extent of about 5ʹ whose spectral properties are compatible with those of PSR+PWN systems. Interestingly, large-scale CO structures across the SNR were observed, indicating the SNR interaction with the ambient molecular gas, which is an important ingredient to enhance the gamma-ray emission due to proton–proton interaction. Our analysis reveals a spectral index above the break energy Γ₂ = 3.13 ±0.27_stat ± 0.12_syst, which may again indicate that the association with an SNR is spurious or that our low-energy analysis suffers from contamination from other neighboring sources in this crowded region.

Finally, the unidentified source 4FGL J1931.1+1656 (Figure 17, middle right) is coincident with the SNR candidate G52.37-0.70 detected in a recent THOR+VGPS analysis (Anderson et al. 2017). However, the spectral index of α = 0.3 ± 0.3 using Very Large Array observations (Driessen et al. 2018) seems to indicate that this candidate is unlikely to be an SNR. The gamma-ray spectrum derived by our analysis resembles that of other SNRs and is not affected by large systematics especially the break energy 203 ± 8 ± 19 and the spectral index above the break Γ₂ = 2.64 ± 0.10 ± 0.04. It is the best candidate for proton acceleration among these three potential SNR associations.

These three regions are extremely complex and would deserve a dedicated analysis at higher energy with Fermi to constrain their location and their association with the corresponding SNR, as well as a spectral analysis over a larger energy interval to definitively constrain the type of radiating particles.

Even more care should be taken for the extended sources 4FGL J1633.0-4746e (Figure 15, top left) and 4FGL J1813.1-1737e (Figure 15, bottom right) since their disk radii of 061 and 06, respectively, in confused Galactic plane regions adds to the complexity of such analysis at low energy. With its large extension, 4FGL J1633.0-4746e overlaps with both the teraelectronvolt PWN candidate HESS J1632-478 and the unidentified source HESS J1634-472, both detected at gigaelectronvolt energies but not included in our list of selected candidates due to their low significance at low energy. This implies that the region contains three sources: a point-like source coincident with HESS J1634-472, an extended source coincident with HESS J1632-478 but with an extension of 0.256° almost twice as large as the teraelectronvolt size, and the very extended source 4FGL J1633.0-4746e overlapping them detected above 10 GeV with a spectral index of 2.25 ± 0.01_stat ± 0.10_syst (Ackermann et al. 2017). Interestingly, our spectral analysis indicates a break at 517 ±18_stat ± 252_syst MeV followed by an index of Γ₂ = 2.11 ±0.15_stat ± 0.12_syst in agreement with the index detected above 10 GeV (though with very large systematics on the break energy due to the diffuse background). The break detected at low energy by our analysis, the hard spectral index Γ₂ consistent with the one detected at higher energy (which seems to indicate a flat spectrum over a large energy range) and the presence of dense clumps in this region traced NH₃(1,1) emission (de Wilt et al. 2017) make this source a very interesting proton accelerator. A dedicated analysis would therefore be very valuable in this case.

The disk radius of $0.60\,\pm \,0.{06}_{\mathrm{stat}}^{^\circ }$ of the Fermi source 4FGL J1813.1-1737e, coincident with the compact teraelectronvolt PWN candidate HESS J1813-178 (Gaussian size of 0.049 ± 0.04° in H.E.S.S. Collaboration et al. 2018b), was first detected by Araya (2018). The authors reported a hard index of 2.07 ± 0.09_stat above 500 MeV compatible with the teraelectronvolt index. This spectrum is compatible with the spectral index Γ₂ = 2.17 ± 0.03_stat ± 0.42_syst derived in our analysis. With such a large extension in the Galactic plane, several sources could contribute to the gigaelectronvolt signal: the PWN powered by PSR J1813-1749 thought to emit at teraelectronvolt energies as seen by H.E.S.S. and HAWC (Abeysekara et al. 2017), the SNR G12.82-0.02 whose contribution to the teraelectronvolt signal was explored by Funk et al. (2007) and the giant SFR W33 that comprises a region of 15' at a distance of 2.4 kpc (Immer et al. 2013). This last hypothesis was considered by Araya (2018), showing that the energetics, extended morphology, and spectrum of the gigaelectronvolt emission are similar to those of the other gamma-ray detected SFR, the Cygnus Cocoon. To firmly establish the presence of protons radiating at gamma-ray energies, such a complex region definitively is worth an individual analysis above 1 GeV to constrain the morphology and a spectral analysis over a larger energy range to model the broadband emission.

Finally, several sources detected by our analysis are completely unassociated and follow-up observations at teraelectronvolt energies and X-rays would be needed to constrain their nature. They all present values of Γ₂ much softer than those of the identified SNRs discussed in Section 4.2. Similarly, the values of Γ₂ − Γ₁ obtained in our analysis are much larger ( ≥ 2.96) than those of the identified SNRs and dense MC regions. This tends to indicate that these sources are not associated with SNR shock acceleration.