Incremental Fermi Large Area Telescope Fourth Source Catalog

The data for the 4FGL-DR3 catalog were taken during the period 2008 August 4 (15:43 UTC) to 2020 August 2 (8:33 UTC) covering 12 years. During most of this time, Fermi was operated in sky-scanning survey mode, with the viewing direction rocking north and south of the zenith on alternate orbits such that the entire sky is observed every ∼3 hr. Starting on 2018 March 16, the Fermi spacecraft was put in safe hold after one solar array drive became stuck. Scientific operations of the LAT were interrupted for more than 3 weeks, by far the longest missing time interval since 2008, and restarted on April 8 in partial sky-scanning mode. ⁸⁸ The Sun more rarely enters the field of view, and during some phases of the ∼53-day precession period of the orbit, the entire sky is not covered every 3 hours. This hampers the solar observations and short-term light curves, but has little impact on the integrated sky coverage.

As in 4FGL, intervals around solar flares and bright gamma-ray bursts (GRBs) were excised. During the additional 4 years, 108 ks were cut due to two successive bright solar flares in 2017 September, and 2.4 ks around five new bright GRBs. The current version of the LAT data is Pass 8 P8R3 (Atwood et al. 2013; Bruel et al. 2018). For DR2 and DR3, we used the P8R3_V3 IRFs instead of P8R3_V2 used to build DR1. The relative effective areas for the point-spread function (PSF) event types were recalibrated on bright sources, ⁸⁹ reducing the systematics from 20% to 5%. The energy range remains 50 MeV to 1 TeV. In 4FGL, the event types were not treated separately above 10 GeV. We now decompose the events into four components up to 30 GeV, resulting in 19 components in total (Table 1), for better precision but at the cost of more disk space and CPU.

The interstellar emission is a very important ingredient when characterizing the point sources, because it is structured and 90% of the photons in the Galactic plane are of diffuse origin. We used exactly the same model as in 4FGL. ⁹⁰ The isotropic spectrum was updated to P8R3_V3. The model for the emission of the Sun and the Moon was also kept the same, neglecting the modulation of the emission along the solar cycle. It became necessary, however, to generate a different effective model for each year (for the light curves). After the solar array drive anomaly during the tenth year, the coverage of the Sun changed significantly (Figure 1). For intervals of about 2 weeks, the Sun does not enter the field of view at all.

Zoom In Zoom Out Reset image size

Several extended sources (existing and new) were the subject of publications since the 4FGL paper. They were not considered in DR2, but were considered for DR3:

• 1.  
  We adopted the Gaussian model proposed by Devin et al. (2020) for the broad supernova remnant (SNR) G150.3+4.5, after we confirmed that it is a better fit than the disk model used in DR1.
• 2.  
  Tibaldo et al. (2018) have studied the Vela X pulsar wind nebula (PWN) and reported that the radio template was superior to the geometric templates. We found that it does not provide a better fit than the DR1 model, but it eliminates six dubious point sources, so it is simpler and we adopted it. The same authors reported a spatial evolution as a function of energy, but we found that their high-energy disk is not significant when using the radio template.
• 3.  
  Araya (2020) has shown that the cluster of DR1 sources around the SNR G279.0+1.1 is actually an extended LAT source. We adopted their ring model.
• 4.  
  HESS J1640−465 is seen as extended by H.E.S.S., but was considered a point source in the DR1 catalog. We adopted the H.E.S.S. Gaussian template, after it was shown by Mares et al. (2021) that it is a better fit to the LAT data as well.
• 5.  
  Two papers discussed the Fermi-LAT emission from the PWN HESS J1825−137. The morphology considered by Araya et al. (2019) was too complex for our catalog. We tried the broader Gaussian proposed by Principe et al. (2020) but found it worsened the fit. We did change the template, however, because we noticed that the one that we used since 2FGL was not the final version of Grondin et al. (2011).
• 6.  
  Ambrogi et al. (2019) have shown that the previous template for the SNR HB 21 was too broad, because it mistakenly included a separate point source. We adopted their smaller template.
• 7.  
  Xin et al. (2019) reported the detection with Fermi-LAT of the extended VERITAS source VER J2227+608, corresponding to the SNR G106.3+2.7 next to the much brighter pulsar PSR J2229+6114. We have included it in the DR3 model.

Table 2 lists the source name, origin, spatial template, and the reference for the dedicated analysis. These sources are tabulated with the point sources, with the only distinctions being that no position uncertainties are reported and their names end in e (see the Appendix). The new and modified extended sources have DataRelease set to 3. Unidentified point sources inside extended ones are indicated as “xxx field” in the ASSOC2 column of the catalog, where “xxx” is the name of the extended source. A total of 78 extended sources are considered in DR3.

The source detection procedure followed the same lines as in 4FGL. It used pointlike (Kerr 2010) and a specialized diffuse model in which the nontemplate features are estimated differently. It started from the previous list of sources (DR1 for DR2, DR2 for DR3), relocalized them over 12 yr of data, looked for peaks in the residual TS maps generated for several spectral shapes, introduced those in the model, refit, and iterated over the full procedure. The resulting list (called uw1216) contained 11,128 seed sources at TS > 10.

That procedure naturally resulted in changing the positions of all sources. Since we wanted the catalog to be incremental, we forced most DR1 and DR2 sources back to their original positions (consistent with their names). This requires associating the uw1216 seed sources with the DR2 catalog. To do that, we started by defining pairs in which each member is the nearest counterpart of the other, within at most 1°. We repeated that procedure after removing the pairs found at the first step, in order to handle faint sources close to much brighter ones (in which the nearest counterpart can be the bright source). We rejected 58 pairs in which the distance between the two members was larger than the largest of the two 99.9% position errors (approximately 1.5 times the 95% error) and 0°.15. We eliminated the uw1216 seeds in a pair with a DR2 source, and formed the input list from the DR2 sources (including 139 with no counterpart in uw1216) and the remaining uw1216 seeds.

In DR2 we actually used the original positions of all DR1 sources strictly, and kept all of them. However a number of sources became much brighter in recent years than they were during the first 8 yr, so they are now much better localized. In recognition of this, we adopted for DR3 the new error ellipse when its area was at least 10 times smaller and TS at least 10 times larger than the same quantities in DR1. We used those very strict conditions to preserve the incremental nature of the catalog. This selected 10 sources, clear outliers when looking at the distribution of ratios of ellipse areas and TS (Figure 2). All are associated with blazars (very variable sources). One of them was in uw8606 (list of seeds for DR1) but too faint to be included in DR1 (it appeared in DR2), so only nine DR1 sources were updated. They have a 4FGL counterpart in the ASSOC_FGL column, but have DataRelease set to 3 because they have moved.

Zoom In Zoom Out Reset image size

We have also adopted the new error ellipses of one DR1 source close to the modified extended source HESS J1825−137 and three faint soft sources (two from DR1, one from DR2) whose 95% semimajor axes R₉₅ were larger than 1° in DR1 or DR2. Three DR1 sources with broad error ellipses R₉₅ ∼ 04 were split into two better localized seeds inside the original error ellipse. We replaced the DR1 sources by the two seeds. The same applied to two DR2 sources. Conversely we replaced two DR1 sources (4FGL J1831.5−0935 and J1830.2−1005) by a single uw1216 seed compatible with PSR J1831−0952.

We applied the same procedure as in 4FGL to eliminate seeds too close to a bright source and inside extended sources. Since we changed the extended sources (Section 3.1), this resulted in deleting several DR1 sources as well, not only new seeds. More precisely, 4FGL J1640.6−4632 was replaced by HESS J1640−465, five DR1 and two DR2 sources were deleted in or close to Vela X, three DR1 and one DR2 sources in the new SNR G279.0+1.1. We also deleted the faint 4FGL J1814.1−1710 next to the extended source HESS J1813−178 (even though this one did not change) because its TS decreased from 30 in DR1 to 16 in DR3, and it was not spectrally distinct from HESS J1813−178.

Following the catXcheck verification (Section 3.6), we have applied several modifications to fit the data better. We replaced 4FGL J0857.7−4507 (in one of the clusters described in Section 6.2) by a stronger uw1216 source about 0°.5 away. We replaced 4FGL J1750.0−3849 just outside the extended source FHES J1741.6−3917 by a stronger point source just inside it.

We have looked at individual faint (TS < 25) sources with a very soft (Γ > 3.1) spectrum, a large error ellipse (R₉₅ > 0°.5), or a stronger close neighbor in uw1216 and deleted three DR1 (4FGL J0533.9+2838, J0313.6−7508, and J2326.5+8555) and three DR2 (4FGL J0750.0+7140, J1112.0+1021, and J1752.2−3002) sources. Nearby seeds appeared in DR3 and fit the data better. Similarly, we have deleted a faint very hard (Γ < 1.3) DR2 source (4FGL J1830.3−1601) that was superseded by a brighter nearby DR3 seed.

As before we manually added 11 known LAT pulsars that could not be localized by the automatic procedure without phase selection. Only one of those (PSR J1731−4744) reached TS > 25 in DR3. It appears with a NULL error ellipse. The 10 others were discarded. A total of 5029 new seeds were entered to the gtlike source characterization for DR3 in addition to the 5690 remaining DR1 and DR2 point sources.

We reassessed the systematic corrections to localization. Up to 4FGL, they were defined globally over the entire sky, so their values were dominated by high-latitude sources (92% of associated point sources in DR1 are at ∣b∣ > 5°). However, confusion and strong interstellar emission make localization much more difficult in the Galactic plane. In DR2 we estimated localization systematics separately at ∣b∣ < 5°, based on 110 pulsars and low-latitude active galactic nuclei (AGN). We found that they should indeed be increased there, to 27″ for the absolute 95% error, and 1.37 for the systematic factor. We did not change the systematics in the high-latitude sky (25″ and 1.06). The DR2 systematics are used again for DR3.

We used the Fermi Tools 1.4.7 analysis suite. ⁹¹ Compared to the Science Tools v11r7p0 that were used in 4FGL, the main difference is the new edisp_bins setting for energy dispersion, which allows for broadening of the energy interval considered in the model space. The Science Tools setting was equivalent to edisp_bins = −1, not enough to fully fill the first and last data bins for the small energy bins that we use (10 per decade). We used edisp_bins = −2 instead, considering two bins for energy dispersion on each side beyond the data bins. The acceptance for the PSF3 event type increases by a factor 2.4 between 50 and 100 MeV, so to account for this rapid change, we used six energy bins in the first component (Table 1) instead of three and set edisp_bins to −4 accordingly.

The likelihood weights were recomputed over 12 yr of data, resulting in slightly smaller weights everywhere. We reoptimized all RoIs, resulting in 1988 RoIs containing up to 10 sources in their core. The maximum of eight sources in the RoI core enforced in 4FGL was not optimal (it led to too many small neighboring RoIs and too many iterations).

For DR1 and DR2, the normalization of one of the diffuse components sometimes became very large (close to the maximum, set to 2) or very small (close to 0) in regions where that component contributed little. This was particularly obvious for the isotropic component in the Galactic ridge. This could happen because of missing features in the interstellar emission model (IEM), leading to a higher isotropic component in the fit. This was not satisfying since the isotropic emission is supposed to be approximately constant over the sky, and the stuck parameters caused trouble in the error estimates. In order to stabilize the diffuse parameters, we replaced the hard limits with Bayesian priors. We chose the priors from the distributions of the parameters (in DR2) in regions where they are dominant (∣b∣ < 10° for the Galactic parameters and ∣b∣ > 30° for the isotropic normalization). This led to the values reported in Table 3. The intrinsic reasons for the large or small values did not go away, of course, so the resulting isotropic normalizations were still larger than one in the Galactic ridge, but they rarely exceeded 1.5. The worst RoIs around (l, b) ∼ (315°, 3°) reached 1.7.

A major difference with the 4FGL procedure is that we retained all previous 4FGL sources (DR1 and DR2, except those discussed in Section 3.2), so that they would not be deleted from the model even if they have TS < 25. The resulting catalog contains 6659 entries, among which 1610 are new (1607 point sources plus the three extended sources discussed in Section 3.1), and 181 have TS < 25 (112 DR1 and 69 DR2 sources).

The light curves over 1 yr bins were recomputed (even the first 8 yr because adding new sources changes the fluxes of DR1 sources). The procedure was very similar to 4FGL, except we considered events up to 1 TeV, and we used a different model for the Sun and Moon at each year (see Section 2.2). Neither change has a major effect. Recomputing light curves over 2 month bins would have been very time consuming, and we showed in the 4FGL paper that the finer time binning was not critical to detecting variability. Monthly (and shorter intervals) light curves are provided in the Fermi-LAT Light Curve Repository.

The threshold on Variability_Index is now 24.725 (same 1% probability of false detection, applied to 12 bins). The number of significantly variable sources increased from 1443 to 1695 (among which 158 are DR2 or DR3 sources). The fraction of variable sources decreased because part of the variability information in 4FGL came from the 2 month light curves. The distribution of fractional variability (Frac_Variability in the FITS file) is similar to 4FGL, peaking between 50% and 90%. The blazars all have fractional variability larger than 10% except 4C +55.17 (4FGL J0957.6+5523). The flux of that source, first discussed by McConville et al. (2011), decreased by about 13% over 10 yr, not showing any flare.

The spectral analysis proceeded along the same lines as in 4FGL, with five important differences detailed below: we changed the parameterization for pulsars, we fit more sources with four parameters, we changed the TS_curv threshold for using a curved spectral shape, we added new columns reporting the peak energy in ν F_ν , and we added one spectral bin to the spectral energy distributions (SEDs).

Most significantly curved sources are still fit with a lognormal function (LogParabola under SpectrumType in the FITS table, hereafter LP):

The significantly curved pulsars are still fit with a subexponentially cutoff power law (PLSuperExpCutoff under SpectrumType in the FITS table, hereafter PLEC). However we noticed that the parameters of the parameterization used in 4FGL (PLSuperExpCutoff2 in the Fermi Tools) were strongly correlated, particularly when the exponential index (b in Equation (2)) was free. This hampered convergence and precluded freeing b except in the very brightest sources. To work around this, we devised a new parameterization (PLSuperExpCutoff4 in the Fermi Tools):

in which the normalization K is directly the flux density at the reference energy E₀, and the shape parameters are the spectral slope and the spectral curvature at E₀.

The development in Equation (3) shows explicitly that for b = 0 the expression converges to LP with d = 2β. For b < 0 the shape (in log) is reversed (it decreases exponentially toward low energies and as a PL toward high energies). The parameters K, Γ_S , d, and b appear as PLEC_Flux_Density, PLEC_IndexS, PLEC_ExpfactorS, and PLEC_Exp_Index in the FITS table, respectively. The reference energy E₀ is chosen so that the error on K is minimal and appears as Pivot_Energy in the FITS table. This parameterization proved indeed much more stable, and the correlation between parameters is considerably reduced.

Figure 3 shows that on average Γ_S is much better defined (error ratio below 1) than the low-energy asymptotic index Γ = Γ_S − d/b used in 4FGL. On average the error on α (Δα) decreased by only 2.3% between 10 and 12 yr, whereas ΔΓ_S is on average 0.626 ΔΓ. Note that 15 pulsars had no ΔΓ at all in DR2 because Γ was stuck at 0. This problem does not exist in the new parameterization. The small error ratio at significance below 40 corresponds to PSR J0622+3749. It had Γ = 0.05 in DR2, very close to the boundary. At low significance, a fraction of the pulsars have ΔΓ in DR2 that is similar or less than ΔΓ_S in DR3. The DR2 errors may have been underestimated because of the large correlation between parameters.

Zoom In Zoom Out Reset image size

Since the new PLEC parameters describe the same mathematical quantities (spectral slope and curvature) as the LP parameters (except the curvature is 2β), it is interesting to compare them. Both spectral models are fit to all sources, so this information is readily available. Figure 4 shows that, as expected, the correlation between α and Γ_S is quite tight. The correlation between 2β and d is not as tight, because the curvature has larger error. On average d < 2β, because the curvature in the PLEC model depends on energy, and the pivot energy (where d is defined) is below the energy at which the two curvatures would be the same. This effect is stronger for bright sources that have lower pivot energy for the same spectral shape.

Zoom In Zoom Out Reset image size

The new PLEC parameterization allows a 4-parameter fit (with free b) in fainter sources. We have applied it to all pulsars with TS >10,000 (28, up from 6 in DR1 and DR2). At TS <10,000 Δb becomes larger than 0.15, which is the natural scatter on b in the brightest pulsars, so freeing b is no longer beneficial. For all other significantly curved pulsars, b is fixed to 2/3 as in 4FGL. For comparison, the median b over the 28 pulsars with free b is 0.51, its weighted (with 1/σ²) average is 0.55, and its intrinsic dispersion is 0.16. The Small Magellanic Cloud has b = 1, as before (Abdo et al. 2010).

Besides 3C 454.3, the five other AGN with TS > 80,000 were modeled with PLEC and free b as well: CTA 102, Mkn 421, S5 0716+71, 3C 279, and PKS 1424−41. The new parameterization contains the LP model (for b = 0), so there is no risk of a worse fit, but nonconvergence can occur when the model is not constrained well enough. The TS threshold for free b is higher in AGN than pulsars, because the curvature is much less in AGN, so Δb is larger. The main objective of fitting PLEC with free b remains to improve the modeling of the surroundings of very bright sources at low energy.

In the DR1 and DR2 catalogs, the sources were represented with a curved spectral model (LP or PLEC) rather than a power law (PL) when (curved spectrum)(PL)) was larger than 9 (3σ). For DR3 we have lowered that threshold to 4 (2σ) for two reasons, illustrated by Figure 5:

• 1.  
  The left-hand plot shows that all sources are curved to some degree, when bright enough that the curvature is significant. This is physically understandable (the PL is an idealized asymptotic case).
• 2.  
  The PL approximation tends to overestimate the source flux outside of the energy range in which the data are well constrained. The right-hand plot quantifies that on the sum of all sources, at three different TS_curv thresholds. Between 200 MeV and 20 GeV, it makes very little difference. High-latitude blazars (blue) are not strongly curved compared to pulsars and other Galactic sources (green and black). But above 20 GeV and below 200 MeV, the TS_curv threshold makes a difference. Most conspicuously, the global spectrum in the plane outside the ridge (green) diverges toward low energies at a threshold of 9 (what was actually used in DR2). The global spectrum is better behaved when the TS_curv threshold is lower.

The implication is that we should choose a TS_curv threshold as low as possible while preserving global convergence of the model. For DR3 we set the threshold to 4 as a compromise because the test with a threshold at 1 required manual intervention, and global convergence could not be guaranteed. This is still a significant improvement as demonstrated in Figure 5 (right).

Zoom In Zoom Out Reset image size

Thanks to the improved statistics and the lowered TS_curv threshold, many more sources are considered significantly curved in DR3 than in DR1. The number of LP spectral models is now 3130 (among which 650 are new sources), while the number of PLEC models is now 258 (among which 2 are new sources). The fraction of sources represented by a curved spectral model increased from 30% to 51% overall.

Changing the spectral models for many sources from PL to curved also has an impact on their energy flux estimate. The PL flux is nearly always larger, due to the integration down to 100 MeV for soft sources (or up to 100 GeV for hard sources Γ < 2). This is quantified in Figure 6. Well-measured pulsars (diamonds at high significance) have on average a 50% larger flux estimate in the PL model than in the PLEC model. The well-measured blazars (plus signs) are in the lower branch at a ratio around 1.1, because their spectra are much less curved than those of pulsars. The diamonds in the lower branch are the very bright blazars fit with PLEC (see above). A ratio less than one occurs when the peak of the curved spectrum is close to either 100 MeV or 100 GeV (the boundaries of the energy interval over which the energy flux is extracted), because the full analysis interval (50 MeV to 1 TeV) extends beyond those boundaries and maintains the PL fit below the peak of the curved model.

Zoom In Zoom Out Reset image size

To provide more information on curved spectra, we now report systematically in the catalog the peak energy in ν F_ν and its uncertainty for all sources (including those not significantly curved) and both models as (Unc_)LP_EPeak and (Unc_)PLEC_EPeak:

The uncertainties are obtained using the covariance matrix. We set those columns to NULL when the peak energy is undefined (when Γ ≥ 2 in PLEC or curvature is upwards) or found too far from the LAT energy range (). E_peak(PLEC) is undefined in only 4 PLEC spectra, but in 516 LP and 1682 PL spectra (mostly because Γ ≥ 2). E_peak(LP) is always defined in PLEC and LP spectra, but remains undefined in 850 PL spectra in which β is found negative or close to 0. The uncertainty on E_peak(LP) is also set to NULL when β is stuck to 1 and has no error (166 occurrences in LP spectra, 6 in PLEC spectra, and 3 in PL spectra). Even when E_peak is formally defined, it should be treated with caution when Γ(PLEC) gets close to 2 or the curvature is close to 0.

In the DR1 and DR2 catalogs, we provided SEDs over N_bands = 7 energy bands. Thanks to the increased statistics after 12 yr (particularly valuable at high energy where the precision is limited by counts rather than background), we can provide them over eight energy bands in DR3, extending to 1 TeV as in the 3FHL catalog (Ajello et al. 2017). We replaced the 30–300 GeV bin by a 30–100 GeV bin and a 100 GeV–1 TeV bin. We computed those two fluxes using unbinned likelihood, as in 4FGL. In the highest-energy band, 172 sources have TS > 25. An example of such a hard source is shown in Figure 7. That source has TS_curv(LP) = 8.70, below the DR2 threshold of 9 but above the DR3 threshold of 4, so it is fit with the LP model in DR3. The significant differences between 1 and 10 GeV are due to the natural variability of that strong source.

Zoom In Zoom Out Reset image size

The systematic uncertainties ⁹² have improved when moving from P8R3_V2 to P8R3_V3 and are now 0.1 in bands 1 and 8 and 0.05 in bands 2–7. We compute the spectral fit quality (SFQ) as before (Equation (5) in Abdollahi et al. 2020). Figure 8 shows how that quantity is distributed. Up to DR2, we flagged (Flag 10 in Table 4) sources with SFQ above the 99.9% quantile (10⁻³ p-value) of the χ² distribution with N_bands −2 degrees of freedom (dof), considering that most sources were fit with PL spectra that have two free parameters. For DR3 that would result in 6 dof. Clearly the data do not follow the χ²(6) distribution, but are much closer to the χ²(4) distribution. This is due to two effects: first, half the sources are now fit with curved spectra (3 free parameters); second, the first (<100 MeV) and last (>100 GeV) bands do not behave as a Gaussian distribution (very large errors in the first band, and very low counts in the last band), so they contribute little to SFQ. For DR3 we flagged sources with SFQ > 18.47, the 10⁻³ p-value of the χ²(4) distribution. The number of sources that trigger this remains small (only 16). The largest SFQ is obtained in PSR J0205+6449 at 49.6. That pulsar indeed has a known PWN (3C 58), which is very significant above the PLEC fit in the two bands above 30 GeV. The other sources flagged with Flag 10 have one spectral band off by more than 3σ.

Zoom In Zoom Out Reset image size

As before, we attempt to flag the sources affected by systematic errors larger than the statistical errors that we report, and which should therefore be treated with more caution. This can impact all source parameters, and even call into question the existence of the faintest flagged sources with significance less than 6σ or 7σ (TS < 50). Because the largest source of systematics is the underlying diffuse emission, the sources close to the Galactic plane, where diffuse emission is both the strongest and most structured, are most affected.

The flags (Flags in the FITS file) are recalled in Table 4 (see the 4FGL paper for the detailed definitions), together with the numbers of sources flagged for each reason and their evolution since DR1. The effect of the underlying IEM was estimated by launching the procedure described in Section 3.3 for a second time using the same seeds but the previous IEM (gll_iem_v06).

We note six changes with respect to the 4FGL procedure:

• 1.  
  For DR1, Flag 1 was applied to the comparison with both pointlike and the alternative IEM. For DR3, the comparison to pointlike was performed with no TS selection in order to apply Flag 3 to sources at TS < 25. Flag 1 was determined only by the comparison with the old IEM, after selecting TS > 25 on both sides.
• 2.  
  Flag 2 mostly relies on the comparison between the new uw1216 pointlike positions of the 4FGL sources (Section 3.2) and the 4FGL positions (uw8606), because no localization was performed over 12 yr with the old IEM. The comparison with the old IEM also flags a few sources (although the seed positions are exactly the same) when among two nearby seeds one survives in DR3 and the other one with the old IEM.
• 3.  
  All sources have been visually screened for potential deficiencies in the underlying diffuse model. The visual screening for diffuse features (Flag 6) depends on the source parameters. It was entirely redone for DR2. To be conservative, we left that flag (and the c suffix) in all the sources that had it in DR1, even if they would not have been flagged in DR2, and added the flag to all the new cases, including the sources already in DR1. For DR3 we checked only the new sources and did not change Flag 6 in the DR1 and DR2 sources. Two sources in the Large Magellanic Cloud (4FGL J0517.9−6930c and J0535.7−6604c) have that flag set as well.
• 4.  
  Flag 12 still applies to PLEC sources with Γ = Γ_S −d/b ≤ 0 (as well as to LP sources with β fixed to 1, as before), even though the new PLEC parameterization (Equation (2)) has no explicit limit at Γ = 0.
• 5.  
  We introduced for DR2 a new Flag (13) to flag explicitly sources from a previous version with TS < 25 in the current version. It now applies to DR2 sources as well (69 have TS < 25 in DR3).
• 6.  
  For DR3, we introduce a new Flag (14) to flag explicitly soft Galactic unassociated (SGU) sources in regions with high source density (Section 6.2). The three strongest SGUs (4FGL J1702.7−5655 at b = −9.2, 4FGL J1823.3−1340 inside HESS J1825−137, and 4FGL J1839.4−0553 inside HESS J1841−055), all at TS > 800, are very clear point sources (good pulsar candidates) and are not flagged.

Overall the fraction of sources with any flag set increased from 23% in DR1 to 32% in DR3. Among the new DR2 and DR3 sources, 75% are flagged in the Galactic plane (∣b∣ < 10°), but only 26% are flagged above the Galactic plane.

In order to verify that the 4FGL-DR3 catalog provides a good representation of the γ-ray sky, we have performed an independent all-sky analysis, named catXcheck, using the p-value statistic (PS) data/model deviation estimator developed by Bruel (2021). catXcheck consists of the analysis of 540 RoIs (12 × 12°, with a pixel size of 0°.1) covering the whole sky. The RoI centers lie on 19 Galactic parallels, whose latitudes range from −90° to 90° with a 9° step. The longitude step is 9° for ∣b∣ ≤ 18°, 10° for ∣b∣ = 27° and 36°, and it is 12°, 15°, 18°, 24°, and 45° for ∣b∣ = 45°, 54°, 63°, 72°, and 81°, respectively. The difference with respect to the catalog main analysis is that we combine all PSF event types together, using all events whose zenith angle is less than 90° to suppress contamination from Earth limb γ rays. The energy range is [100 MeV, 1 TeV], and we use 10 bins per decade.

The source model of each RoI comprises the Galactic diffuse emission, the isotropic template, the Sun and Moon steady emission templates, as well as all point-like and extended sources from DR3, within 5 + 0.015(σ_src − 4) degrees of the RoI border, where σ_src is the source significance (Signif_Avg in the FITS version of the catalog). That formula selects sources outside of the RoI up to 5° for the faintest ones and to 15° for the brightest ones, corresponding respectively to the 68% and 95% containment radii at 100 MeV.

Since we want to check the goodness of the DR3 information, the spectral parameters of the DR3 sources are fixed to their catalog values, and the only free parameters of the binned-likelihood fit are the ones of the Galactic diffuse emission (power-law correction) and isotropic (normalization) templates. After the fit, we compute the 3D map of the number of photons predicted by the source model and compute the PS map in order to detect any significant deviation between the observed and predicted 3D maps.

As explained in Bruel (2021), the PS estimator is computed for each spatial pixel of a RoI, and rather than comparing the observed and predicted count spectra at that pixel, we compare spatially integrated count spectra around that pixel, with an integration radius varying with energy as , where ⊕ is the sum in quadrature. We use the point-source optimized integration parameterization (which is very close to the PSF 68% containment energy dependence): p₀ = 4°, p₁ = 0.9 and p₂ = 0°.1, and we merge the energy bins to have an energy binning of 0.3 in . PS is defined such that , where the p-value is the probability that the statistical fluctuations can reach a level of deviation as large as the one observed in the data, under the assumption that the model represents the data. The PS sign is the sign of the sum of the count spectra residuals in units of σ. The likelihood weights are taken into account in the PS computation. We chose to use the PS map to investigate data/model deviations because, contrary to the TS map, it is sensitive to both positive and negative deviations.

For one pixel, the 3σ, 4σ, and 5σ thresholds correspond to PS = 2.57, 4.20 and 6.24. Since the level of correlation between pixels is negligible with p₂ = 0°.1, we only have to take into account the number of trials when deriving the all-sky significance thresholds. The number of RoIs and their centers have been chosen so that the overlap between neighboring RoIs ensures that, for any given direction in the sky, there is at least one RoI containing this direction and providing a distance to the RoI border larger than ∼1°.5. As a consequence, when combining the results of the 540 PS maps, we loop over the RoIs and the RoI pixels, and we only consider the pixels for which the RoI provides the largest distance to the RoI border. This allows us to avoid double counting the sky directions. The total number of considered pixels is 4,175,468 and adding its logarithm to the one-pixel significance thresholds gives the all-sky 3σ, 4σ, and 5σ thresholds: PS = 9.19, 10.82, and 12.86.

The PS distributions for samples of sources in three different ranges of Galactic latitude are shown in Figure 9. The high-latitude (∣b∣ > 35°) distribution follows the 10^−∣x∣ expectation, whereas some deviations are visible in the low-latitude (∣b∣ < 5°) and mid-latitude (5° < ∣b∣ < 35°) distributions. To show the locations of these deviations, we build an all-sky HEALPix (Górski et al. 2005) map in Galactic coordinates with N_side = 256 from the 540 PS maps. Figure 10 shows the resulting all-sky PS map. Overall the data/model agreement is good over the whole sky. Using a 4σ threshold, we find two positive and five negative deviations. The positive ones are both located inside bright extended sources: 4FGL J1923.2+1408e (W 51C) and 4FGL J0822.1−4253e (Puppis A). They occur because the geometric templates that we use for those sources are far from perfect, and we do not allow additional point sources inside them (Section 3.2). Regarding the negative deviations, the one at a distance of ∼1° from the Galactic center corresponds to only one pixel. The other four are clusters of several pixels with PS < −10.82 and correspond to extended deviations explaining the negative tail in the mid-latitude histogram in Figure 9. They all are close to large molecular clouds: one can be associated with Cepheus, one with Perseus, and two with Orion B.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since these extended negative deviations seem to be related to the Galactic diffuse emission, we also compute the 540 PS maps with the set of PS spatial integration parameters optimized for large extended deviations (p₀ = 5°, p₁ = 0°.8, and p₂ = 1°). With p₂ = 1°, the level of correlation between pixels is expected to increase. Assuming that pixels within 0°.5 of each other are highly correlated, we divide the number of trials by 69, and the 3σ, 4σ, 5σ thresholds become PS = 7.4, 9, 11. The corresponding all-sky PS map, shown in Figure 11, exhibits many deviations above 4σ. A preliminary analysis of the negative ones shows that many of them are correlated with some components of the IEM, especially the local CO template (Galactocentric radius range 7–9 kpc), which may indicate that allowance for variation of the X_CO factor among the local clouds is necessary. The positive deviations are less numerous and less significant. Like the negative ones, some are likely due to imperfections in the modeling of the Galactic diffuse emission, but some could correspond to new extended sources.

Zoom In Zoom Out Reset image size

Although a complete analysis of these extended deviations is outside the scope of this paper, we try to quantify their impact on the DR3 catalog. In order to do so, we look at the 540 PS maps and count the number of pixels with ∣PS∣ > 7.4 and derive the expected number of DR3 sources in these regions from the average source density in the ∣PS∣ < 7.4 region. In the 40 RoIs centered along the Galactic plane, the fraction of negative and positive deviation pixels is 2% and 0.9%, respectively. The numbers of observed and expected DR3 sources are 15 and 24.2 in the negative deviation pixels, and 33 and 10.1 in the positive deviation pixels. In the RoIs farther from the Galactic plane, the fraction of negative and positive deviation pixels is 0.15% and 0.03%, respectively. There are no DR3 sources in the negative deviation pixels, while 8.4 are expected. We conclude that the catalog results are biased in the negative and positive deviation regions, but the effect is small and the number of missed or spurious sources is at most ∼20.

The TS of the 5048 remaining DR1 sources in DR3 increased on average compared to DR1. Among sources that are fit with the PL model in both catalogs (no TS boost from switching to a curved model), TS increased on average by 33% at high latitude (∣b∣ > 10°) and by 22% closer to the Galactic plane. The difference between the two is due to the lower log-likelihood weights, which partly offset the exposure increase close to the Galactic plane where diffuse emission is strong, particularly for soft sources. Even at high latitude, the TS ratio is less than expected from the exposure increase (50%). This is not due to the lower weights (the TS ratio does not depend on index) but to the selection bias (sources tend to be brighter in the interval in which they were defined) and signal-splitting with new sources. The statistical uncertainties on the parameters decreased by 15% at high latitude and 9% close to the plane.

The random character of adding new data and the variable nature of the γ-ray sources inevitably lead to broadening the TS distribution, in addition to shifting it to larger values. The sharp selection threshold at TS = 25 in DR1 becomes a tail of sources at TS < 25 in DR3. Figure 12 shows the TS distribution outside the Galactic plane, where 74 DR1 sources have TS < 25. Over the entire sky, 112 DR1 sources have TS < 25 in DR3, among which 15 have TS < 16, and the smallest TS is 9.8 (4FGL J1816.4−0057).

Zoom In Zoom Out Reset image size

The energy flux decreased on average by 5% (0.2σ) between DR1 and DR3 (considering only PL sources to avoid the effect of switching to a curved model). This is again a selection bias dominated by faint sources: close to the detection threshold, DR1 sources can get either brighter or fainter in DR3, but faint sources in DR3 cannot be fainter in DR1 because they would have TS < 25. Real variability is apparent for bright sources (TS > 1000). The scatter on the energy flux ratio is 5σ on variable DR1 sources, but only 1σ on nonvariable ones. There is no offset and less scatter (only 0.8σ overall) on the difference of PL indices because source fluxes vary strongly, but spectral shapes only slightly.

The 1607 new DR2 and DR3 sources are on average very close to the detection threshold (median TS of 34). This is obvious in Figure 12. There are a few newly active bright blazars at TS > 100, and even three at TS > 800 (4FGL J0154.5−6604, J1544.3−0649, and J2250.0−1250 associated with PKS 0153−663, NVSS J154419−064913, and PKS 2247−131 respectively), but the majority are faint sources that rise above the detection threshold thanks to the additional four years. Outside the Galactic plane, the confusion limit is still distant, because ever-deeper exposures reach to higher energies at which the PSF is better. Among the 1042 faint (TS < 100) new sources at ∣b∣ > 10°, only 7% have the confusion flag (Flag 5 in Table 4) set.

Figure 13 (left) shows that the energy flux distribution of the new sources peaks between 1 and 2 × 10⁻¹² erg cm⁻² s⁻¹. The majority (140) of the 156 new DR2 and DR3 sources considered variable (mostly blazars) are outside the Galactic plane (∣b∣ > 10°). Their average flux is larger than that of the nonvariable sources. This is largely a selection bias (variability cannot be detected in the very faint sources).

Zoom In Zoom Out Reset image size

Figure 13 (right) compares the distribution of PL indices of the new sources with that of the DR1 ones. The variable new sources are soft on average, similar to flat-spectrum radio quasars (FSRQs). The same effect is observed in the 1FLT catalog (Baldini et al. 2021), even more strongly because detections in 1FLT are based on one-month data sets. The nonvariable new sources have a peak on the soft side as well (probably similar to the variable ones, but not bright enough to be significantly variable). But they also have a shoulder on the hard side. This is because adding more exposure is more beneficial at high energies (in the count-dominated regime) than at low energies, and therefore helps with detecting hard sources (like BL Lac objects).

To understand the effects of the analysis changes introduced in DR3 with respect to 4FGL, we have considered their effects on the DR2 data set (10 years). To that end, we started with the same seeds as the DR2 catalog, changed each element in sequence (in the order of the list below), and compared each intermediate result with the previous one.

• 1.  
  We first switched from the Science Tools to the Fermi Tools, and from the P8R3_V2 to the P8R3_V3 data and Instrument Response Functions (changing the isotropic templates accordingly). This led to essentially no change at all, as expected.
• 2.  
  Setting edisp_bins to −2 (Section 3.3) when accounting for energy dispersion had a significant effect. Prior to that, the model counts were underestimated by about 5% in the first energy band (50 to 100 MeV), so the model fit was too high by the same amount at low energy. We found 16 fewer sources with the new settings, due to a slight TS reduction in soft sources. The energy flux decreased by 2% on average. The largest effect was spectral: we found 155 more LP spectra (+10%). Hard PL sources (Γ < 2.1) were unchanged, but Γ decreased more and more for softer sources, reaching ΔΓ = − 0.1 at Γ = 3. Note that this is specific to the catalog setting that uses many components for the summed likelihood in relatively narrow energy bands (Table 1). Analyses using few broader energy bands are not as sensitive to edisp_bins.
• 3.  
  Introducing diffuse priors had very little effect on the sources. We found 7 fewer sources and 11 fewer LP spectra. In the Galactic ridge (where the effect of the priors is strongest), the average energy flux of sources decreased by 2% (0.1σ), and the average PL index decreased (became harder) by 0.005 (0.06σ). This is completely innocuous.
• 4.  
  Lowering the TS_curv threshold from 9 to 4 (Section 3.4) had an important effect on the sources, besides the collective effect illustrated in Figure 5 (right). The number of LP spectra increased by 1000 (+58%). The total number of sources increased by 127 (+2%) because more LP sources allowed more soft neighbors to reach the detection threshold. The energy fluxes of sources that switched from PL to LP decreased by 1.8σ (well above the statistical precision) as illustrated in Figure 6. Their neighbors picked up most of that flux, so their energy fluxes increased. They also became somewhat softer in the process (by 0.02 on average). On the other hand, the neighbors of those neighbors became a little harder, so that the global average PL index did not change.

In conclusion, the change with the greatest consequences was lowering the TS_curv threshold. We are confident that the new setting is closer to the physical reality.

For completeness, we recapitulate below the census of associations of new DR2 sources already given in the DR2 document (Ballet et al. 2020), with some updates on blazar classes:

• 1.  
  two pulsar candidates, PSR J1439−5501 and PSR J1904−11;
• 2.  
  one globular cluster (GLC), NGC 362;
• 3.  
  one high-mass binary (HMB), PSR B1259−63; although detected by the LAT during its 2010 periastron passage (Abdo et al. 2011), this source is included in an FGL catalog for the first time thanks to the gain in significance from its 2017 passage (Johnson et al. 2018);
• 4.  
  three SNRs, 3C 397, SNR G001.4−00.1, and SNR G003.7−00.2; they are considered as reliable associations (as opposed to SPPs, see below) because their radio sizes are quite small (4′–14′) and their positions well match those of the 4FGL sources;
• 5.  
  two star-forming regions, NGC 346, which is the brightest star-forming region in the SMC, and Sh 2–152;
• 6.  
  two galaxies, IC 678 and NGC 5380;
• 7.  
  one starburst galaxy, Arp 299;
• 8.  
  two radio galaxies, NGC 3078 and NGC 4261;
• 9.  
  281 blazars, including 176 blazars of unknown type (BCUs), 66 BL Lacertae objects (BLLs), and 39 FSRQs; ⁹⁵
• 10.  
  16 SPPs (SPP stands for “SNR, Pulsar or PWN” and refers to sources of unknown nature but overlapping with known SNRs or PWNe and thus candidate members of these classes);
• 11.  
  30 sources of unknown nature (UNK, ∣b∣ < 10° sources solely associated with the likelihood-ratio method from large radio and X-ray surveys).

The pulsar PSR J1757−2421, formerly associated with 4FGL J1756.6−2352 in DR1, is now associated with 4FGL J1757.9−2419. The source 4FGL J1756.6−2352 is now a SPP.

Below is a census of the associations for the new DR3 sources:

• 1.  
  three pulsars with detected pulsations, including two young pulsars, PSR J0922+0638 and PSR J1731−4744 (Smith et al. 2019), and one MSP, PSR J1455−3330;
• 2.  
  one MSP candidate, PSR J1957+2516;
• 3.  
  five GLCs, NGC 1851, NGC 5286, NGC 6205, NGC 6712, and M 54;
• 4.  
  two PWNe, PWN G327.1−1.1 and PWN G54.1+0.3;
• 5.  
  one SNR, the Kepler SNR; the DR3 source 4FGL J1730.4−2131 was associated by the Likelihood Ratio method with NVSS J173036−212910; this NVSS source is actually part of the Kepler SNR (extended for the Very Large Array at 3′ diameter), so we report the SNR as the counterpart (see also Acero et al. 2022);
• 6.  
  three novae, V1369 Cen (Cheung et al. 2016), V906 Car (Jean et al. 2018; Aydi et al. 2020), and YZ Ret (MGAB-V207; Li et al. 2020); the position of the DR3 source (4FGL J0358.5−5432) is closer to that of a blazar candidate, PMN J0358−5434, than to that of YZ Ret, but the light curve clearly points to the nova association;
• 7.  
  one high-mass X-ray binary, 1RXS J172006.1−311702 (Esposito et al. 2014); the nature of this source is debated as it was classified as a cataclysmic variable in Fortin et al. (2018);
• 8.  
  one radio galaxy, LEDA 55267, classified as a Fanaroff–Riley type 0 (Baldi et al. 2019) with a LAT detection reported in Paliya (2021), and one AGN, NGC 6454;
• 9.  
  306 blazar candidates including 219 BCUs, 51 BLLs and 36 FSRQs;
• 10.  
  22 SPPs;
• 11.  
  20 sources of unknown nature.

• 1.  
  A specific class (GC) has been created for the Galactic center (4FGL J1745.6−2859), replacing its former SPP class.
• 2.  
  The increase of localization systematic uncertainties for sources lying at ∣b∣ < 5° (Section 3.2) has led to the association of five extra BCUs, five SPPs, and the suppression of the association of 4FGL J1804.9−3001 with the GLC NGC 6528. Four sources classified as unknown in DR1 are now classified as BCUs because they gained Bayesian-based associations. Two other formerly unknown sources are now associated with a pulsar and a binary star (Kleinmann’s star).
• 3.  
  Two associations with pulsar candidates for 4FGL J1618.7−4633 and 4FGL J1720.8−1937 have been discarded on the realization that their estimated γ-ray luminosities would exceed their spin-down powers by factors greater than 30.
• 4.  
  Pulsations have been detected for 22 more DR1 sources, including 10 PSR and 12 MSP (e.g., Clark et al. 2021; Tabassum et al. 2021).
• 5.  
  A total of 29 associations with MSP from the West Virginia University list ⁹⁶ have been added (three replacing associations with BCUs and one with a SPP), along with one association with a young pulsar candidate.
• 6.  
  The pulsar PSR J1909−3744 was mistakenly associated with 4FGL J1912.2−3636. The position of this LAT-detected pulsar was used as a seed, but the resulting 4FGL source ended up too far away from the seed to make the association plausible. The 4FGL source is now unassociated.
• 7.  
  The association of the MSP binary recently discovered by Wang et al. (2020) with 4FGL J0935.3+0901 has been implemented. Its class is binary (BIN).
• 8.  
  Following Bogdanov & Halpern (2015), 4FGL J1544.5−1126 is a candidate transitional MSP binary (1RXS J154439.4−112820). Its class, BCU in DR1, has been changed to “low-mass X-ray binary” (LMB). Another candidate transitional MSP binary, CXOU J110926.4−650224, was found (Coti Zelati et al. 2019, 2021) in a systematic search for X-ray counterparts to unassociated sources in the preliminary LAT 8 yr point-source list (FL8Y). While the position of the X-ray source lies well within the 95% error ellipse of the FL8Y source, it is outside of the ellipse of the corresponding DR3 source 4FGL J1110.3−6501. In consequence, ⁹⁷ we have not retained this association in our tables.
• 9.  
  Another transitional MSP binary, 1SXPS J042749.2−670434, displaying simultaneous optical, X-ray, and γ-ray eclipses has been identified by Kennedy et al. (2020) in 4FGL J0427.8−6704 (previously unassociated).
• 10.  
  Gamma² Velorum (4FGL J0809.5−4714; Martí-Devesa et al. 2020) has been classified as a binary system.
• 11.  
  We have implemented the results of Corbet et al. (2019), who discovered that 4FGL J1405.1−6119 is a HMB.
• 12.  
  The source 4FGL J1832.9−0913, associated with the TeV source HESS J1832−093, has been classified as a HMB (see Martí-Devesa & Reimer 2020).
• 13.  
  Swihart et al. (2021) found a candidate redback MSP binary in a follow-up observation of 4FGL J0940.3−7610. Other redbacks have been discovered in a similar fashion in 4FGL J0407.7−5702 (Miller et al. 2020) and 4FGL J0540.0−7552 (Strader et al. 2021). These associations have been classified as LMBs.
• 14.  
  Similarly, Li et al. (2021) found a black widow binary at the position of 4FGL J0336.0+7502, which has then also been classified as a LMB.
• 15.  
  Three periodic optical and X-ray sources were recently discovered positionally consistent with DR1 sources. These redback MSP candidates have been classified as BIN. They are 4FGL J0212.1+5321 (Li et al. 2016), 4FGL J0523.3−2527 (Strader et al. 2014), and 4FGL J0838.7−2827 (Halpern et al. 2017).
• 16.  
  The tentative association of 4FGL J0647.7−4418 with the HMB RX J0648.0−4418 called out in DR1 has been replaced by association with the BCU SUMSS J064744−441946 following the multiwavelength investigation of Marti et al. (2020).
• 17.  
  Following Järvelä et al. (2020), we have reclassified TXS 2116−077 (4FGL J2118.8−0723) as a Seyfert galaxy instead of a NLSY1.
• 18.  
  The latest version of the Radio Fundamental Catalog ⁹⁸ has enabled the association of six previously unassociated sources and two SPPs with blazar candidates.
• 19.  
  Recent follow-up observations of 4FGL blazars (Peña-Herazo et al. 2017, 2019, 2020, 2021a, 2021b, 2021c; Desai et al. 2019; Rajagopal et al. 2021) have enabled the classification of 240 former BCUs, 2 AGNs, and 2 UNKs into 214 BLLs and 30 FSRQs.
• 20.  
  Two FSRQs (PKS 0736−770, TXS 1530−131) have been reclassified as BCUs and one other (Two Micron All Sky Survey, 2MASS J02212698+2514338) as a BLL.

Note: A typo has been found in § 6.2 of the 4FGL paper concerning the name of the star cluster associated with the extended H ii region encompassing 4FGL J1115.1−6118. The star cluster is NGC 3603 (and not NGC 4603).

Table 5 lists the class tallies of the whole DR3 catalog.

In this section, the spectral properties for the main source classes in the whole 4FGL-DR3 sample are presented. We use these properties to shed light on the nature of the unassociated sources in the next section. Although some classes, most notably pulsars, show significant levels of spectral curvature, the PL photon index, Γ, is a convenient parameter enabling comparison across classes. Figure 14 compares the photon-index distributions for different established classes. Pulsars show remarkably consistent soft spectra, along with GLCs. MSP have somewhat harder spectra on average compared to PSR and GLCs. The observed difference between the MSP and GLC distributions seems to challenge the conventional wisdom that MSP dominate the GLC’s γ-ray emission (e.g., Harding et al. 2005; Venter et al. 2009), and is possibly due to different ambient stellar radiation fields comptonized by the MSP electrons (Bednarek & Sitarek 2007). The spectra of SNRs and PWNe are notably harder than those of pulsars, in agreement with the observation that several SNRs and almost all PWNe are detected in the TeV domain. The distribution for blazars spans a broad range of photon index, the distributions for FSRQs and BLLs exhibiting moderate overlap (see Ajello et al. 2020). Blazars constitute the only abundant class showing a hard-spectrum component.

Zoom In Zoom Out Reset image size

Turning to the spectral curvature, while some blazars show significant curvature, as do most pulsars, the two classes occupy different branches in the LP_TS_curv versus TS plane (Figure 15 top). These branches merge at the low-TS end as curvature cannot be assessed with high significance in that case. Despite the fact that PSR and MSP show slightly different spectra as can be seen from their photon-index distributions, both fall on the same branch in Figure 15 top.

Zoom In Zoom Out Reset image size

Drawing on the class-wise discriminating power of the photon-index distributions illustrated in Figure 14, the distribution for high-latitude unassociated sources is compared to that of BCUs, which are the most relevant blazar set given the low TS, in Figure 16 top. The good agreement between the two distributions supports the conclusion that most high-latitude unassociated sources are indeed blazars. ⁹⁹ This conclusion is slightly qualified when considering the loci of the unassociated sources in the TS versus LP_TS_curv plane (Figure 15 bottom). A significant fraction (≃70%) of the higher-TS (TS > 100) unassociated sources fall onto the pulsar branch, strongly suggesting that they are MSP. This idea is also supported by considering the evolution of the photon-index and variability-index distributions with TS (Figure 17). While for the lower-TS bins, the distributions quite closely match those of BCUs, clear differences are observed for the higher-TS bin (100 < TS < 500), where the distributions are compatible with those expected for MSP. In conclusion, their photon indices, significance of spectral curvature, and lack of variability make many of the high-latitude, 100 < TS < 500 unassociated sources good MSP candidates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since some AGN catalogs suffer from significant incompleteness close to the plane due to extinction, some background blazars detected through the plane could elude association and contribute to the unassociated sample. Based on the tallies observed at high latitude and the difference in flux limits due to the brighter diffuse emission background close to the plane, the number of ∣b∣ < 10° blazars is estimated to be 340 ± 20, while 399 blazars are actually reported as associations in DR3. Comparing the low-latitude photon-index distribution with that obtained for higher latitudes (normalized so that the total yield of Γ < 2.4 sources match), an excess of soft sources (Γ ≃ 2.5) is found in the low-latitude sample relative to the other. From the difference of the spectral-index distributions, this excess, which was less pronounced in 3FGL and is very likely due to a contamination from a nonblazar population, amounts to 75 ± 4 sources. In conclusion, an excess of soft-spectra sources stands out in the low-latitude blazar sample. It can be attributed to a contamination from another population, likely of Galactic origin. After subtracting this excess from the number of reported blazars, the remaining number of detected blazars agrees reasonably well with that extrapolated from the high-latitude population (324 versus 340 ± 20). Few blazars with missed associations are thus expected to contribute to the low-latitude unassociated sample.

Unassociated sources around the Galactic plane (referred to as GUs) account for a growing fraction of the detected sources, from 12% in DR1 to 24% in the new sample. The features presented here were already present in DR1 but have become more salient with the larger population.

The GU latitude distribution, shown in Figure 18, shows a narrow Gaussian-like component (≃210 sources) centered on the plane, and two broad shoulders (≃760 sources) extending out to about ∣b∣ = 10° (these components will be referred to as “spike” and “shoulders” respectively). The spike is narrower ( = 0.005) than the young-pulsar distribution ( =0.012 for sources with ). The loci of the unassociated sources close to the plane exhibit a notable clustering. Some high-density regions (“hot spots”) are readily seen in Figure 19, while other areas are mostly devoid of unassociated sources. The clustering aspect also manifests itself in the large fraction (30%) of GUs with Flag 5 (i.e., close to a brighter neighbor). High-density regions have been determined using a HEALPix tesselation with N_side = 8 produced in Galactic coordinates. The average density of GUs per pixel in the three pixel rows spanning the Galactic Plane is about 7. We define the high-density pixels as those encompassing more than 14 GUs, i.e., twice the average value.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The photon-index distribution is compared to that of PSR in the lower panel of Figure 16. The distribution is remarkably narrow and contrasts sharply with that observed at higher latitude (Figure 16 top). It peaks around Γ = 2.5, with a tail toward lower Γ values. This tail might be due to blazars, SNRs, or PWNe (Figure 14). The high-Γ end of the distribution is peculiar. No established class of Galactic γ-ray emitters exhibits such a soft component. We will refer to sources in this soft component as SGUs.

A dedicated flag (Flag 14) has been devised to indicate SGUs that are in regions of relatively high source density as seen in Figure 19. Flag 14 is set for SGUs lying inside a high-density pixel, having Γ > 2.4, TS < 500, and a closest neighbor (classified either as GU, SPP, or UNK) less than 2° away. Since most SPP or UNK sources probably belong to the same population as GUs as ascertained by their similar photon-index distributions (peaking around Γ = 2.5), they are also flagged if the above conditions are filled. Flag 14 is set for 446 SGUs, 54 SPPs, and 49 UNKs. ¹⁰⁰ SGUs were already significantly present in DR1, but their number in high-density regions has now increased by 50% (see census of Flag 14 in Table 4).

Returning to the spectral properties, the photon-index distribution is essentially constant across the plane, as illustrated in Figure 20 for different (not mutually exclusive) sky regions. In this figure, “clusters” refers to the high-density pixels mentioned above. The unassociated sources preferably lie within the pulsar branch in the LP_TS_curv versus TS plane (Figure 21), supporting their spectral resemblance with pulsars. Their distribution of the LP curvature parameter β is also consistent with that of pulsars. Despite these similarities, other factors distinguish GUs and pulsars. First, the photon-index distributions do not closely match. The significant clustering observed also disfavors the possibility that GUs are mostly pulsars. Moreover, the low-b unassociated sources outnumber the detected PSR by about a factor of 9, making the scenario that so many pulsars have gone undetected so far implausible. The final blow to this possibility comes from their broad latitude distribution, which does not match that of PSR (nor any of the distributions of known γ-ray emitters for that matter). Older pulsars (>10⁶ yr), which could have drifted farther off the plane over time, exhibit a wider Galactic-latitude distribution, more compatible with that observed. However, these pulsars should be undetectable in γ-rays as they would have crossed the “pair-creation death line” (see Abdo et al. 2013).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Most GUs exhibit a low source-to-background ratio (30% of them having the corresponding Flag 4 set). Whether they originate from mismodeled Galactic diffuse emission becomes a valid question. This emission is not fully understood, and its modeling suffers many uncertainties. To produce the Pass 8 model, ¹⁰¹ the contributions of different components, H i, CO, dark neutral medium, and inverse Compton, whose properties are established at larger wavelengths, were fitted to the data in 9 rings of increasing Galactocentric distances. An ad hoc component called the patch encapsulates the contributions of nontemplate components like the Galactic center excess, Loop I, or the Fermi bubbles. Other, yet-to-be-identified large scale structures may be contributing (see Jóhannesson & Porter 2021, for a search performed at high latitude), and a comprehensive analysis of the patch still needs to be carried out. Understanding the patch low-energy behavior would be of particular relevance to this work. The patch covers a limited sky area and was smoothed out not to absorb features with spatial scales less than 4°. Looking for spatial correlations between the SGUs and features present in the different components, the best match was obtained with the patch component at low energy (see Figure 19). The SGU latitude distribution exhibits a north/south symmetry and a relative extension (Figure 18) that find a good correspondence in the patch at ∣b∣ < 10°. Comprehending the photon excess underlying the patch will certainly be key in unveiling the nature of the SGUs.

Below we list some possibilities regarding the SGUs and lay out avenues for further exploring them.

• 1.  
  As mentioned above, SGUs may result from mismodeled Galactic diffuse emission. Allowing more flexibility to the model fitting by leaving free individual diffuse components, like H i and CO, would constrain the robustness of the SGU detections, if the templates used for the model are accurate on scales that are finer than the RoI size, while inaccurate on larger scales.
• 2.  
  A possibility to be considered is that SGUs trace the large-scale regions associated with fresh injection of cosmic-rays, leading to locally harder γ-ray emission than assumed in the model components. A worthwhile test would be to investigate (e.g., by means of simulations) whether the spurious point-like sources spawned from a possible, not-accounted-for diffuse excess would share common spectral features with the observed SGUs.
• 3.  
  Some SGUs could be spatially extended if they are related to structures in the diffuse emission. This would introduce some spectral distortion when analyzing them under the assumption that they are point-like, making them softer than they really are.
• 4.  
  Some correlation with the general structures of interstellar clouds is found by visual inspection, suggesting that the emission ascribed to SGUs could be due to molecular gas unaccounted for on small scales, e.g., because of saturated CO lines. This situation is reminiscent of that encountered in the EGRET era, which led to the discovery of “dark gas” (Grenier et al. 2005) and the suppression of about 30% of unassociated EGRET sources (Casandjian & Grenier 2008) once a modified diffuse emission model was considered. Assessing the effect of saturated CO lines using ¹³CO data would shed light on this possibility.
• 5.  
  Alternative ways could be pursued to build the patch (diffuse emission that is not accounted for by standard components of the interstellar emission model) and distinguish it from point sources (e.g., Selig et al. 2015).
• 6.  
  The alternative to mismodeled diffuse emission, namely the existence of an abundant, entirely new class of sources remains highly speculative. Follow-up multiwavelength studies (as in Bruzewski et al. 2021), e.g., of regions of high SGU densities or surrounding the brightest SGUs, could help elucidate this issue.