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

Cosmic rays are mostly composed of protons accelerated to relativistic speeds. When those protons encounter interstellar material, they produce neutral pions, which in turn decay into gamma-rays. This offers a compelling way to identify the acceleration sites of protons. A characteristic hadronic spectrum, with a low-energy break around 200 MeV, was detected in the gamma-ray spectra of four supernova remnants (SNRs), IC 443, W44, W49B, and W51C, with the Fermi Large Area Telescope. This detection provided direct evidence that cosmic-ray protons are (re-)accelerated in SNRs. Here, we present a comprehensive search for low-energy spectral breaks among 311 4FGL catalog sources located within 5° from the Galactic plane. Using 8 yr of data from the Fermi Large Area Telescope between 50 MeV and 1 GeV, we find and present the spectral characteristics of 56 sources with a spectral break confirmed by a thorough study of systematic uncertainty. Our population of sources includes 13 SNRs for which the proton–proton interaction is enhanced by the dense target material; the high-mass gamma-ray binary LS I+61 303; the colliding wind binary η Carinae; and the Cygnus star-forming region. This analysis better constrains the origin of the gamma-ray emission and enlarges our view to potential new cosmic-ray acceleration sites.


INTRODUCTION
The acceleration site of protons, the main components of cosmic rays, is one of the most fundamental topics of high energy astrophysics.The strong shocks associated with supernova remnants (SNRs) are widely believed to accelerate the bulk of Galactic cosmic rays (E < 10 15 eV) through the diffusive shock acceleration mechanism (e.g.Drury 1983).Indeed, accelerated cosmic rays interact with surrounding matter and produce π 0 mesons which usually quickly decay into two gamma rays, each having an energy of 67.5 MeV in the rest frame of the neutral pion.In turn, the gamma-ray number spectrum F (E) is symmetric at this same energy in log-log representation (Stecker 1971) which then leads to a γ-ray spectrum in the usual E 2 F (E) representation rising below 200 MeV and approximately tracing the energy distribution of parent protons at energies greater than a few GeV.This characteristic spectral feature, often referred to as the "pion-decay bump", uniquely identifies proton acceleration since leptonic γ-ray production mechanisms such as bremsstrahlung and inverse Compton (IC) emission require fine tuning to produce a similar feature.Esposito et al. (1996) explored this hypothesis by studying the γ-ray emission from SNRs, and potential associations of γ-ray sources with five radio-bright shell-type SNRs were reported using data taken by the EGRET instrument on board the Compton Gamma Ray Observatory.More recently, this signature of protons was detected in five SNRs interacting with molecular clouds (MCs) and detected at gamma-ray energies by Fermi -LAT: IC 443 and W44 (Giuliani et al. 2011;Ackermann et al. 2013), W49B (H.E. S. S. Collaboration et al. 2018a), W51C (Jogler & Funk 2016) and HB 21 (Ambrogi et al. 2019), although in this last source both the leptonic and hadronic processes are able to reproduce the γ-ray emission.Finally, the young SNR Cassiopeia A was also analyzed at low energy and Yuan et al. (2013) derived an energy break at 1.72 +1.35  −0.89 GeV which is better reproduced by a hadronic scenario.More details on this characteristic feature observed in the gamma-ray emission are provided in Appendix A, showing a stronger signature for a soft proton injection index (Γ = 2.5) than for a hard index (Γ = 1.5).
Electrons can also radiate at gamma-ray energies via the inverse Compton scattering and bremsstrahlung processes.It has been demonstrated, for the supernova remnants interacting with molecular clouds cited above, that the large gamma-ray luminosity is difficult to explain via inverse Compton scattering.In addition, the steep gamma-ray spectrum detected at low energy requires additional breaks in the electron spectrum if we consider a model in which electron bremsstrahlung is dominant.Accurate estimation of the spectral characteristics of a γ-ray source at low energy is therefore crucial since it probes the nature of the particles (electrons or protons) emitting these gamma rays.However, the analysis of sources below 100 MeV is complicated due to large uncertainties in the arrival directions of the gamma rays, which lead to confusion among point sources and difficulties in separating point sources from diffuse emission.Thus, catalogs released by the Fermi -LAT Collaboration have focused on energies greater than 100 MeV until the 4FGL catalog (Abdollahi et al. 2020) expanded the lower bound to 50 MeV.This allows better constraint of low-energy spectra, but since the 4FGL upper energy bound is 1 TeV, the spectral model for most sources is dominated by data with energies above a few hundred MeV.In addition, the spectral representation of sources in the 4FGL catalog considered three spectral models: power law (PL), power law with sub-exponential cutoff, and log-normal (or log-parabola, hereafter called LP).This means that any source presenting a spectral break will be represented by a log-normal shape which may not adequately represent the low-energy behavior.Similarly, sources presenting two spectral breaks, as it is the case for W49B (H.E. S. S. Collaboration et al. 2018a) will be represented with a log-normal shape that better describes the high-energy interval due to the better angular resolution and increased effective area at these high energies.This directly implies that the description of the low-energy spectral parameters of a source requires a dedicated spectral analysis.
In this paper we use 8 years of Pass 8 data to analyse 311 Galactic sources detected in the 4FGL catalog and search for significant spectral breaks between 50 MeV and 1 GeV.The paper is organized as follows: Section 2 describes the LAT and the observations used, Section 3 presents our systematic methods for analyzing LAT sources in the plane at low energy, Section 4 discusses the main results and a summary is provided in Section 5.

Fermi-LAT
The Fermi -LAT is a γ-ray telescope which 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 direction measurement of the incident γ rays), a CsI(Tl) crystal calorimeter (for energy measurement), and an anti-coincidence detector to identify the background of charged particles.Since the launch of the spacecraft in June 2008, the LAT event-level analysis has been upgraded several times to take advantage of the increasing knowledge of how the Fermi -LAT functions as well as the environment in which it operates.Following the Pass 7 data set, released in August 2011, Pass 8 is the latest version of the 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 the 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.

Data selection and reduction
We used exactly the same dataset as that used to derive the 4FGL catalog of sources, namely 8 years (2008 August 4 to 2016 August 2) of Pass 8 SOURCE class photons.This means that similarly to the 4FGL dataset, our data were filtered removing time periods when the rocking angle was greater than 90 • and intervals around solar flares and bright 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 γ 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 f ermitools version 1.2.23 and f ermipy (Wood et al. 2017) version 0.19.0.We used only binned likelihood analysis because 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 γ 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 1 .Our binned analysis includes three logarithmically spaced energy bins between 50 MeV 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 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 (FSSC)2 .In the following, we fit the normalizations of the Galactic diffuse and the isotropic components.

Effect of the energy dispersion
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 GeV 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 a 8-year observation time using the gtobssim tool included in the LAT f ermitools.For these simulations, we assumed a point source spatial model located at (RA, Dec (J2000): 94.51 • , 22.66 • ) and a smooth broken power-law spectral model of the form: where α = 0.1, the break energy E br = 245 MeV and the spectral indices Γ 1 = 0.57, Γ 2 = 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, Γ 1 and Γ 2 and fitted a Gaussian on each distribution.The centroid of the Gaussian fit together with their size 1 https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Pass8_edisp_usage.html 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), (10 energy bins, edisp bins = −3).As can be seen on this Figure, the main effect is on Γ 1 , as expected.If the energy dispersion is not taken into account (edisp bins = 0), the spectrum is less steeply falling at low energy and the spectral index Γ 1 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 Γ 1 will be biased towards lower value which will create artificially a stronger break at low energy.This is directly due to the fact that the energy resolution varies with energy.This imposes to choose 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.

List of candidates
This analysis intends to find new cosmic-ray acceleration sites in our Galaxy.When cosmic-ray protons accelerated by a source penetrate into 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.For AGNs, we removed all subclasses, namely flat-spectrum radio quasars (FSRQs), BL Lac-type objects (BLLs), blazar candidates of uncertain type (BCUs), radio galaxies (RDGs), narrowline Seyfert 1 (NLSY1s), steep spectrum radio quasars (SSRQs), Seyfert galaxies (SEYs) 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 a list of 311 candidates reported in the Appendix B.

Input source model construction
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 f ermipy.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.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 res-olution 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 f ermipy.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 TS = 2(ln L 1 − ln L 0 ), where L 0 and 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 source of interest, 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 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 Γ = 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 down 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 source of interest, we test iteratively for each of them the improvement of the log-normal representation with respect to the power-law model.The log parabola model is defined as: where N 0 is the overall normalisation factor to scale the observed brightness of a source, E 0 is a fixed scale energy (kept at 300 MeV in our analysis), and α, β are left free, which adds one degree of freedom with respect to the power-law representation.The improvement of the log-parabola model (LP) with respect to the power law one (PL) is performed by determining TS LP = 2(ln L LP − ln L 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 logparabola 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 characterize properly the region.Because the 4FGL Catalog rejects most point sources found inside extended sources, this leaves many residuals which translate into sources.

Spectral energy breaks
Once the ROI is well characterized, we first test the TS value of our source of interest in our energy range (50 MeV to 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 value is reported in Table B2.If the TS of the source of interest is above 25, we move on and we fix all sources located more than 5 • from the source of interest and we test the spectral curvature of our source of interest.To ensure that the curvature is real and affects several energy bins as it 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 B2 and we stop the analysis of this source.It is the case for 167 sources in our sample.If the value is above 9, we then test a smoothly broken power law following Equation 1 where N 0 is the differential flux at E 0 = 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 power-law model (the break energy E br and a second spectral index Γ 2 ).The improvement with respect to the power law one is determined by TS SBPL = 2(ln L SBPL − ln L PL ).Since this test requires the addition of two degrees of freedom to the fit and diffusive shock acceleration predicts Γ 2 ∼ 2, we also test the improvement of the smooth broken power law with the second index fixed at 2 with respect to the power law one TS SBPL2 = 2(ln L SBPL2 − ln L PL ).We require TS SBPL > 12 or TS SBPL2 > 9 (implying a 3σ improvement for 2 and 1 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 re-do the whole process as illustrated in the flow chart 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 , TS SBPL2 for each of them is reported in Table B2.

Diffuse and IRF Systematics
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 3 .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 flow chart 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 only is 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.Note-Columns 2 and 3 are obtained with the galactic diffuse background rescaled for Pass 8 Source (gll iem v06.fits) and provide values of the improvement of the smooth broken power-law representation with respect to the power-law model TSSBPL and the improvement of the smooth broken power-law representation when fixing Γ2 = 2 called TSSBPL2 as defined in Section 3.2.Columns 4, 5, 6 and 7 provide the same values of TSSBPL and TSSBPL2 for the two bracketing IRFs.Stars denote spectral breaks that are robust to all tests.See Section 3.4 for more details.
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. 2012a), 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 extrema 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 source of interest 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 from 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+4031earound the bright PSR J2021+4026.Overall, 56 sources among the 77 sources detected with the standard IEM and IRFs are confirmed with 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 reported 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 Γ 1 .However, the value of Γ 2 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 uni-formly in log space from 50 MeV to 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 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.Note-Results of the maximum likelihood spectral fits for sources showing significant breaks confirmed by the systematic studies.These results are obtained using a smooth broken power law representation.Columns 2, 4, 6 and 8 report the integrated flux, the break energy and the photon indices Γ1 and Γ2 of the source fit in the energy range from 50 MeV to 1 GeV following Equation 1. Columns 3, 5, 7 and 9 report the statistic and systematic uncertainties on these spectral parameters.Note-For the source classes SNR, PWN, SPP, SFR, BIN and HMB, we add both the firm identifications reported in the 4FGL catalog as well as the associations (capital and lower case letters as seen in Column 6 of Table B2) 4. DISCUSSION

Population study
We detected 56 4FGL γ-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.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.
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 with 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.On the other hand, the fraction of sources associated with supernova remnants (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% (5 out of 311) to 7.1% (4 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 years and these sources often present variable γ-ray emission.A more thorough analysis of these sources would need to be done.Finally, only one star-forming region 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 γ-ray sources showing a significant energy break.The lack of these sources at latitude 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 γ-ray sources.
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 MeV 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 Γ 1 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 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 Γ 2 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 Γ 2 − Γ 1 in Figure 8 (right), a peak at ∼ 0.9 is highly pronounced for SNRs.This tends to show that the values obtained on Γ 2 and Γ 2 − Γ 1 could be used in the future to probe the type of particles radiating in a γ-ray source.Note-Columns 2 and 3 are derived from the Assoc1 and Assoc2 columns of the 4FGL Catalog.The latter provides an alternate designation or an indicator as to whether the source is inside an extended source.

Supernova remnants and molecular clouds
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 spectral energy distributions, 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.Among the 56 sources with significant breaks, one can see from the 4FGL Classification column listed in Table B2 that ten 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 to W49B and thus can be safely identified as a 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 +1.35  −0.89 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 supernova remnants interacting with molecular clouds (MCs).These molecular clouds are excellent targets for cosmic-ray interactions and subsequent piondecay.The hadronic scenario was also preferred for other LATdetected SNRs interacting with MCs, though their γray analysis starting above a few hundred MeV did not allow rejection of a leptonic scenario: the SNR HB3 and the W3 HII complex (Katagiri et al. 2016b), Distribution of the 4FGL significance between 300 MeV and 1 GeV for the 311 sources selected (black line), the 247 sources with T S > 25 in our pipeline (black dashed line), the 77 sources with significant breaks (blue line) and the 56 confirmed cases by our studies of systematics (red line).
S147 (Katsuta et al. 2012), HB9 (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 turn-over 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(HB9, 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 HB3 but it should be noted that its γ-ray emission is much fainter than the adjacent molecular cloud 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. 2016a).In this case, the brightest gamma-ray peak is spatially correlated with the Rosette Nebula, a young stellar cluster and molecular cloud complex located at the edge of the southern shell of the SNR which has a role similar to W3 for the HB3/W3 complex.The interaction between the SNR and the molecular cloud 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 non-thermal emission of the gamma Cygni SNR (Fraija & Araya 2016;Fleischhack 2019), associated with the source 4FGL J2021.0+4031e, also suggested that the γ-ray emission (analyzed above 100 MeV) might be of hadronic nature with enhanced GeV emission spatially coincident with the TeV 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 γ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−5909eassociated 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 PWN at these energies.

Constraints on other identified sources
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 PeV 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 cosmicray accelerators.In this respect, the detection in 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 CRs over a scale of ∼ 50 pc.Our analysis revealed a spectral break for the star-forming region analyzed, 4FGL J2028.6+4110e(see Figure 18 top left), which is associated with the cocoon.A very hard index Γ 1 = 1.00 ± 0.02 stat ± 0.37 syst is detected up to a break energy at 383 ± 13 stat ± 138 syst MeV, followed by a spectral index Γ 2 = 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.A complete modeling of the source at gamma-ray energies is beyond the scope of this paper but our results tend to favour the hadronic scenario, thus reinforcing the long-standing  hypothesis that massive-star-forming regions 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+6113associated with the high-mass γ-ray binary (HMB) LS I +61 303 (Figure 10 top right), the HMB 4FGL J1018.9−5856(Figure 12 bottom left), 4FGL J1045.1-5940associated with the colliding wind binary η Carinae (Figure 12 bottom right) and 4FGL J2032.6+4053associated with the microquasar Cyg X-3 (Figure 18 top right).However, this last source presents the highest value of spectral index Γ 1 (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 ob- servations (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 γ-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.4.4.Interesting new cases: potential proton accelerators ?
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% vs 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 to 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 CRs.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 Γ 2 = 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.3detected as a nearly circular shape and angular diameter of ∼ 40 arcmin at radio energies (Junkes et al. 1992) while Swift and Suzaku X-ray observations allowed 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 supernova remnant.Suzaku observations also revealed diffuse emission with extent of about 5 arcmin 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 Γ 2 = 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 VLA observations (Driessen et al. 2018) seems to indicate that this candidate is unlikely to be an SNR.The γ-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 = 2.64 ± 0.10 ± 0.04.It is the best candidate for proton acceleration among these three potential SNR association.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 0.61 • and 0.6 • respectively in confused Galactic plane regions adds to the complexity of such analysis at low energy.With its large extension, 4FGL J1633.0−4746eoverlaps with both the TeV PWN candidate HESS J1632−478 and the unidentified source HESS J1634−472, both detected at GeV 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 TeV size, and the very extended source 4FGL J1633.0−4746eoverlapping 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 = 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 Γ 2 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 3 (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 ± 0.06  Araya (2018).The authors reported a hard index of 2.07 ± 0.09 stat above 500 MeV compatible with the TeV index.This spectrum is compatible with the spectral index Γ 2 = 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 GeV signal: the PWN powered by PSR J1813−1749 thought to emit at TeV energies as seen by H.E.S.S. and HAWC (Abeysekara et al. 2017), the SNR G12.82−0.02whose contribution to the TeV signal was explored by Funk et al. (2007) and the giant starforming region (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 GeV 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 broad-band emission.Finally, several sources detected by our analysis are completely unassociated and follow-up observations at TeV energies and X-rays would be needed to constrain their nature.They all present values of Γ 2 much softer than those of the identified SNRs discussed in Section 4.2.Similarly, the values of Γ 2 − Γ 1 obtained in our analysis is much larger (≥ 2.96) than those of the identified SNRs and dense molecular cloud regions.This tends to indicate that these sources are not associated with SNR shock acceleration.

SUMMARY
Using 8 years of Pass 8 LAT data between 50 MeV and 1 GeV, we have analyzed 311 4FGL sources located within 5 • from the Galactic plane and detected 77 sources with significant spectral breaks.We carried out a thorough study of the systematics associated with the diffuse Galactic background and with the effective area for each of them and we confirmed the spectral break for 56 of them.With 13 SNRs identified within this sample of 56 sources, SNRs are the dominant class of sources showing significant breaks at low energy.Only five binaries are included in the sample of 311 sources analyzed but four of them show a significant break at low energies.This seems to indicate that binaries could also have a significant contribution.The spectral char-acteristics were also evaluated for these 56 sources.The break energy of the sources ranges uniformly between 100 MeV and 550 MeV.However, a clear pattern is detected in the spectral index Γ 2 of the sources which tends to center at 2.3 for the population of 13 identified SNRs.Similarly, the value of Γ 2 − Γ 1 tends to center at ∼ 1 for the same population of sources.This provides an interesting way to constrain the nature of the radiating particles.Our analysis also provides three interesting new proton accelerator candidates: 4FGL J1931.1+1656 is coincident with the SNR candidate G52.37−0.70 detected in a recent THOR+VGPS analysis, the extended source 4FGL J1633.0−4746eoverlapping the TeV PWN candidate HESS J1632−478 and the unidentified source HESS J1634−472, and the extended source 4FGL J1813.1−1737ecoincident with the compact TeV PWN candidate HESS J1813−178 and the star-forming region W33.The current and future observations of the LAT are thus crucial to probe the spectral characteristics of a source at low energy, providing excellent targets of proton acceleration for current and future Cherenkov telescopes such as CTA.
The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis.These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden.Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Etudes Spatiales in France.Work at NRL is supported by NASA.MLG acknowledges support from Agence Nationale de la Recherche (grant ANR-17-CE31-0014).

A. THE PION-DECAY BUMP SIGNATURE
As already discussed in the main text, when accelerated protons interact with the interstellar matter, they produce neutral pions which in turn decay into gamma rays.This will create a characteristic signature at low energy in the gamma-ray spectrum called the "pion-decay bump signature".To better understand how this signature is characterized in our energy interval of interest (50 MeV -1 GeV), we have used the python package naima (Zabalza 2015) to derive LAT residual TS maps in equatorial coordinates and significance units (left) and spectral energy distributions (right) of IC 443 (top) and W44 (bottom) between 50 MeV and 1 GeV.In the residual TS maps, all white crosses indicate the 4FGL sources included in the model of the region.For the SEDs, the blue points and butterflies are obtained in this analysis while the black points and dashed line are from the 4FGL catalog.The red lines take into account both the statistical and systematic errors added in quadrature.A 95% C.L. upper limit is computed when the TS value is below 1.
the gamma-ray emission produced by proton-proton interaction.To do so, naima uses an implementation of the analytical parametrizations of the energy spectra and production rates of gamma rays from Kafexhiu et al. (2014), which is accurate within 20%.The inclusive π 0 production cross section is included as a combination of the experimental data cross sections at low energies, the Geant 4.10.0cross section at intermediate energies and at higher energies the hadronic model Pythia 8.18 as default.We applied naima to three different power-law distributions of protons with spectral index Γ 1 varying between 1.5 and 2.5.The results presented in Figure A1 (Left) are in perfect agreement with those published in Ackermann et al. (2013) and show that a very steep spectrum is expected below the break energy at ∼ 200 MeV.This Figure also highlights that the pion-decay bump signature might be more difficult to detect for a hard proton distribution (red curve) than for steep injection spectra.This might in turn increase the systematic errors on the derived break energy.Finally, this Figure demonstrates that the restricted energy interval of our analysis does not allow constraining the spectral index of the parent distribution since the gamma-ray spectra trace the energy distribution of parent protons at energies greater than 1 GeV.The upper bound of our energy interval was chosen since middle-aged SNRs commonly exhibit a high energy spectral break at around 1-10 GeV (see the case of W28 with a break at 1 GeV reported by Abdo et al. 2010) and a simple broken power-law model would not apply anymore above 1 GeV.To test for this effect, we applied naima to the same power-law distributions of protons adding a break at 1 GeV in their distributions.We assumed that Γ 2 = Γ 1 + 1. Figure A1 (Right) demonstrates that the break energy of the gamma-ray emission detected in our energy interval is not affected.However, this break significantly impacts the spectral index derived assuming a simple broken power-law model.In a second step, we used the gamma-ray emission obtained with naima assuming a power-law distribution of protons with Γ 1 = 2.0 and 2.5 to produce 200 Fermi simulated data files for each index over 8 years using the gtobssim tool included in the LAT f ermitools.We then analyzed these simulations using the f ermitools following the same procedure as with the real data, assuming a smooth broken power-law spectral model with α = 0.1 (see Equation 1), except that the SED was produced for 12 energy bins instead of 6 to reflect the high statistics of our simulations (which mimics the flux of IC 443).We did not introduce any diffuse background in our simulations to clearly show how a proton spectrum will be reconstructed at low energies with the LAT in the absence of systematic errors.Figure A2 presents the gamma-ray spectrum derived for one of these simulations demonstrating that the spectral index derived above the break does not trace the parent proton distribution.The error bars are extremely small due to the high flux of the simulated source and the absence of diffuse background.This Figure also shows that the smooth broken power-law model used in our analysis with α = 0.1 reproduces the gamma-ray spectrum well.A smoothness parameter α = 0.5 was also tested but it does not significantly improves the likelihood of the fit.Finally, one can see that the injection spectral index does not seem to impact the break energy and the index Γ 1 of our smooth broken power-law fit.The only parameter affected is the index Γ 2 .To confirm this trend, we plotted the distributions of the 200 reconstructed values of the break energy, Γ 1 and Γ 2 for each injection spectral index, and fitted a Gaussian on each distribution as can be seen in Figure A2 (Right) for the case of the break energy.The results are presented in Table A1 confirming that the only parameter affected by the different injection spectral index is Γ 2 .This study demonstrates that no steep spectrum is predicted for a standard injection spectral index.The only way to produce the steep spectra observed for some of the candidates detected in our analysis would be to include an energy break at (or below) 1 GeV in the injection proton spectrum as shown in Figure A1 (Right).

B. LIST OF THE 311 GALACTIC PLANE SOURCES ANALYZED
Table B2 provides the list of all candidates analyzed.Columns 2 and 3 provide the Galactic longitude and latitude of the 311 candidates.Columns 4, 5 and 6 give the curvature significance, the significance between 300 MeV and 1 GeV and the source class reported in the 4FGL catalog.Columns 7, 8, 9, 10 provide the values obtained in our analysis concerning the TS of each source, the improvement of the log-normal representation with respect to the power-law model TS LP as defined in Section 3.2, the improvement of the smooth broken power-law representation with respect to the power-law model TS SBPL and the improvement of the smooth broken power-law representation when fixing Γ 2 = 2 called TS SBPL2 .Left: Gamma-ray spectrum derived using our analysis pipeline from one simulation of a power-law distribution of protons with spectral index of 2.0 predicted by naima (Zabalza 2015).The fit assuming a smooth broken power-law model (see Equation 1) with α = 0.1 is presented with the blue curve.The fit assuming α = 0.5 is presented with the red curve.The best fit of a similar simulation assuming a proton injection index of 2.5 is presented with the green line.Right: Distribution of the 200 values of energy break fitted for a proton injection spectral index of 2.0.The black line represents our best Gaussian fit.
Table A1.Results of the Gaussian fits of the 200 naima simulations assuming a proton injection spectral index of 2.0 and 2.5.Note-The first column indicates the parameter fitted (energy break, Γ1 and Γ2) , the second and third column presents the mean and sigma of the distribution obtained for a proton injection of 2.0 and 2.5, respectively.

Figure 1 .
Figure 1.Effect of the number of energy bins and value of edisp bins on the reconstructed values of the spectral index Γ1 (left), Γ2 (middle) and the break energy (right) of the broken power-law model of IC 443.Four configurations are tested: 12 energy bins and edisp bins = 0 (blue circle), 12 energy bins and edisp bins = −1 (orange triangle), 12 energy bins and edisp bins = −3 (green diamond), 10 energy bins and edisp bins = −3 (red square).

Figure 2 .
Figure 2. Flow chart illustrating the individual analysis procedure of each source of interest (SOI) located in a 20 • × 20 • region of interest.See text for further details.

Figure 3 .
Figure 3. Latitude (top) and longitude (bottom) distributions of the 311 sources selected (black line), the 247 sources with T S > 25 in our pipeline (black dashed line), the 77 sources with significant breaks (blue line) and the 56 confirmed cases by our studies of systematics (red line).

Figure 5 .
Figure 5. Pie charts showing the classes of sources analyzed (Left) and those for which a significant break is detected (Right).The class names are those used in the 4FGL catalog: SNR stands for Supernova remnant, PWN for pulsar wind nebula, SFR for star-forming region, BIN for binary, HMB for high-mass binary.The designation SPP indicates potential association with SNR or PWN.The UNK class includes low-latitude blazar candidates of uncertain type associated solely via the Likelihood-Ratio (LR) method.

Figure 6 .
Figure 6.Distribution of sources in Galactic coordinates.Light gray markers indicate the 311 sources analyzed in this paper.Coloured markers indicate the position of the 56 sources for which a significant break is detected: yellow for UNIDs, blue for SNRs, orange for PWNe, green for SPPs, red for SFR, pink for UNKs and purple for BIN/HMB.The boundary of the latitude selection is 5 • .

Figure 7 .
Figure 7. Break energy for the 56 sources confirmed by our studies of systematics (black line) and for the identified SNRs and/or crushed molecular clouds (dotted line, see Section 4.2).The green line indicates the value of the break energy obtained using simulations based on the naima package for a proton injection index of 2.0 and the two green dotted-dashed lines indicate the one sigma confidence interval derived (more details in the Appendix A). (Top) Individual values; (Bottom) Corresponding histograms.

Figure 8 .
Figure 8. Γ1 (blue line, left), Γ2 (red line, left) and Γ2 − Γ1 (right) distributions for the 56 sources confirmed by our studies of systematics.In all cases, the dotted line corresponds to the same distribution presented with the solid line but restricted to SNRs (see Section 4.2).The green line indicates the value of Γ1 obtained using simulations based on the naima package for a proton injection index of 2.0 and the two green dotted-dashed lines indicate the one sigma confidence interval derived (more details in the Appendix A).

Figure 9 .
Figure9.LAT residual TS maps in equatorial coordinates and significance units (left) and spectral energy distributions (right) of IC 443 (top) and W44 (bottom) between 50 MeV and 1 GeV.In the residual TS maps, all white crosses indicate the 4FGL sources included in the model of the region.For the SEDs, the blue points and butterflies are obtained in this analysis while the black points and dashed line are from the 4FGL catalog.The red lines take into account both the statistical and systematic errors added in quadrature.A 95% C.L. upper limit is computed when the TS value is below 1.

Figure A1 .
Figure A1.Left: Gamma-ray spectra produced by a power-law distribution of protons with spectral index of 1.5 (red), 2.0 (black), 2.5 (blue) as predicted by naima(Zabalza 2015).Right: Same figure assuming that an energy break is present in the particle distribution at 1 GeV.The energy interval analyzed in this work is defined by the green area.

Figure A2.
Figure A2.Left: Gamma-ray spectrum derived using our analysis pipeline from one simulation of a power-law distribution of protons with spectral index of 2.0 predicted by naima(Zabalza 2015).The fit assuming a smooth broken power-law model (see Equation1) with α = 0.1 is presented with the blue curve.The fit assuming α = 0.5 is presented with the red curve.The best fit of a similar simulation assuming a proton injection index of 2.5 is presented with the green line.Right: Distribution of the 200 values of energy break fitted for a proton injection spectral index of 2.0.The black line represents our best Gaussian fit.

Table 3 .
Fractions of photons from neighboring sources and diffuse background affecting all confirmed sources showing a significant break Note-Column 1 indicates the distance (in degrees) of the nearest neighboring source.Columns 2 and 3 report the ratio, in the pixel at the source position, between the predicted number of photons from the source of interest with respect to those of the galactic and isotropic diffuse background (Nsoi/N diff ), and to those of all neighboring sources Nsoi/Nsrcs, respectively.

Table 4 .
Spectral parameters of all confirmed sources showing a significant break

Table 5 .
Summary of source classes

Table 6 .
Candidates with firm associations reported in the 4FGL Catalog

Table B2 .
List of selected Galactic plane candidates