Fermi Large Area Telescope Fourth Source Catalog

The LAT detects γ-rays in the energy range from 20 MeV to more than 1 TeV, measuring their arrival times, energies, and directions. The field of view of the LAT is ∼2.7 sr at 1 GeV and above. The per-photon angular resolution (point-spread function, PSF; 68% containment radius) is ∼5° at 100 MeV, improving to 08 at 1 GeV (averaged over the acceptance of the LAT), varying with energy approximately as E^−0.8 and asymptoting at ∼01 above 20 GeV (Figure 1). The tracking section of the LAT has 36 layers of silicon strip detectors interleaved with 16 layers of tungsten foil (12 thin layers, 0.03 radiation length, at the top or Front of the instrument, followed by four thick layers, 0.18 radiation lengths, in the Back section). The silicon strips track charged particles, and the tungsten foils facilitate conversion of γ-rays to positron-electron pairs. Beneath the tracker is a calorimeter composed of an eight-layer array of CsI crystals (∼8.5 total radiation lengths) to determine the γ-ray energy. More information about the LAT is provided in Atwood et al. (2009), and the in-flight calibration of the LAT is described in Abdo et al. (2009d) and Ackermann et al. (2012a, 2012c).

Zoom In Zoom Out Reset image size

The LAT is also an efficient detector of the intense background of charged particles from cosmic rays and trapped radiation at the orbit of the Fermi satellite. A segmented charged-particle anticoincidence detector (plastic scintillators read out by photomultiplier tubes) around the tracker is used to reject charged-particle background events. Accounting for γ-rays lost in filtering charged particles from the data, the effective collecting area at normal incidence (for the P8R3_SOURCE_V2 event selection used here; see below)⁸⁵ exceeds 0.3 m² at 0.1 GeV, 0.8 m² at 1 GeV, and remains nearly constant at ∼0.9 m² from 2 to 500 GeV. The live time is nearly 76%, limited primarily by interruptions of data taking when Fermi is passing through the South Atlantic Anomaly (SAA, ∼15%) and readout dead-time fraction (∼9%).

The data for the 4FGL catalog were taken during the period 2008 August 4 (15:43 UTC) to 2016 August 2 (05:44 UTC) covering eight years. During most of this time, Fermi was operated in sky-scanning survey mode (viewing direction rocking north and south of the zenith on alternate orbits). As in 3FGL, intervals around solar flares and bright GRBs were excised. Overall, about two days were excised due to solar flares, and 39 ks due to 30 GRBs. The precise time intervals corresponding to selected events are recorded in the GTI extension of the FITS file (Appendix A). The maximum exposure (4.5 × 10¹¹ cm² s at 1 GeV) is reached at the North celestial pole. The minimum exposure (2.7 × 10¹¹ cm² s at 1 GeV) is reached at the celestial equator.

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) and a narrower PSF at high energies. Both aspects are very useful for source detection and localization (Ajello et al. 2017). We used the Source class event selection, with the Instrument Response Functions (IRFs) P8R3_SOURCE_V2. Pass 8 introduced a new partition of the events, called PSF event types, based on the quality of the angular reconstruction (Figure 1), with approximately equal effective area in each event type at all energies. The angular resolution is critical to distinguish point sources from the background, so we split the data into those four categories to avoid diluting high-quality events (PSF3) with poorly localized ones (PSF0). We split the data further into 6 energy intervals (also used for the spectral energy distributions in Section 3.5) because the extraction regions must extend further at low energy (broad PSF) than at high energy, but the pixel size can be larger. After applying the zenith angle selection (Section 2.3), we were left with the 15 components described in Table 2. The log-likelihood is computed for each component separately, then they are summed for the SummedLikelihood maximization (Section 3.2).

The lower bound of the energy range was set to 50 MeV, down from 100 MeV in 3FGL, to constrain the spectra better at low energy. It does not help detecting or localizing sources because of the very broad PSF below 100 MeV. The upper bound was raised from 300 GeV in 3FGL to 1 TeV. This is because as the source-to-background ratio decreases, the sensitivity curve (Figure 18 of Abdo et al. 2010e, 1FGL) shifts to higher energies. The 3FHL catalog (Ajello et al. 2017) went up to 2 TeV, but only 566 events exceed 1 TeV over 8 yr (to be compared to 714,000 above 10 GeV).

The zenith angle cut was set such that the contribution of the Earth limb at that zenith angle was less than 10% of the total (Galactic + isotropic) background. Integrated over all zenith angles, the residual Earth limb contamination is less than 1%. We kept PSF3 event types with zenith angles less than 80° between 50 and 100 MeV, PSF2 and PSF3 event types with zenith angles less than 90° between 100 and 300 MeV, and PSF1, PSF2, and PSF3 event types with zenith angles less than 100° between 300 MeV and 1 GeV. Above 1 GeV, we kept all events with zenith angles less than 105° (Table 2).

The resulting integrated exposure over 8 yr is shown in Figure 2. The dependence on decl. is due to the combination of the inclination of the orbit (256), the rocking angle, the zenith angle selection, and the off-axis effective area. The north–south asymmetry is due to the SAA, over which no scientific data is taken. Because of the regular precession of the orbit every 53 days, the dependence on R.A. is small when averaged over long periods of time. The main dependence on energy is due to the increase of the effective area up to 1 GeV, and the addition of new event types at 100 MeV, 300 MeV, and 1 GeV. The off-axis effective area depends somewhat on energy and event type. This, together with the different zenith angle selections, introduces a slight dependence of the shape of the curve on energy.

Zoom In Zoom Out Reset image size

Selecting on zenith angle applies a kind of time selection (which depends on direction in the sky). This means that the effective time selection at low energy is not exactly the same as at high energy. The periods of time during which a source is at zenith angle <105° but (for example) >90° last typically a few minutes every orbit. This is shorter than the main variability timescales of astrophysical sources in 4FGL and is, therefore, not a concern. There remains however the modulation due to the precession of the spacecraft orbit on longer timescales over which blazars can vary. This is not a problem for a catalog (it can at most appear as a spectral effect and should average out when considering statistical properties), but it should be kept in mind when extracting spectral parameters of individual variable sources. We used the same zenith angle cut for all event types in a given energy interval, to reduce systematics due to that time selection.

Because the data are limited by systematics at low energies everywhere in the sky (Appendix B), rejecting half of the events below 300 MeV and 75% of them below 100 MeV does not impact the sensitivity (if we had kept these events, the weights would have been lower).

2.4.1. Diffuse Emission of the Milky Way

2.4.2. Isotropic Background

2.4.3. Solar and Lunar Template

2.4.4. Residual Earth Limb Template

This section describes the generation of a list of candidate sources, with locations and initial spectral fits. This initial stage uses pointlike. Compared with the gtlike-based analysis described in Sections 3.2–3.7, it uses the same time range and IRFs, but the partitioning of the sky, the weights, the computation of the likelihood function, and its optimization are independent. The zenith angle cut is set to 100°. Energy dispersion is neglected for the sources (we show in Section 4.2.2 that it is a small effect). Events below 100 MeV are not useful for source detection and localization, and are ignored at this stage.

3.1.1. Detection Settings

The process started with an initial set of sources, from the 8 yr FL8Y analysis, including the 75 spatially extended sources listed in Section 3.4, and the three-component representation of the Crab (Section 3.3). The same spectral models were considered for each source as in Section 3.3, but the favored model (power law, curved, or pulsar-like) was not necessarily the same. The point-source locations were also re-optimized.

The generation of a candidate list of additional sources, with locations and initial spectral fits, is substantially the same as for 3FGL. The sky was partitioned using HEALPix⁹⁰ (Górski et al. 2005) with N_side = 12, resulting in 1728 tiles of ∼24 deg² area. (Note: references to N_side in the following refer to HEALPix.) The RoIs included events in cones of 5° radius about the center of the tiles. The data were binned according to energy, 16 energy bands from 100 MeV to 1 TeV (up from 14 bands to 316 GeV in 3FGL), Front or Back event types, and angular position using HEALPix, but with N_side varying from 64 to 4096 according to the PSF. Only Front events were used for the two bands below 316 MeV, to avoid the poor PSF and contribution of the Earth limb. Thus, the log-likelihood calculation, for each RoI, is a sum over the contributions of 30 energy and event-type bands.

All point sources within the RoI and those nearby, such that the contribution to the RoI was at least 1% (out to 11° for the lowest energy band), were included. Only the spectral model parameters for sources within the central tile were allowed to vary to optimize the likelihood. To account for correlations with fixed nearby sources, and a factor of three overlap for the data (each photon contributes to ∼3 RoIs), the following iteration process was followed. All 1728 RoIs were optimized independently. Then the process was repeated, until convergence, for all RoIs for which the log-likelihood had changed by more than 10. Their nearest neighbors (presumably affected by the modified sources) were iterated as well.

Another difference from 3FGL was that the diffuse contributions were adjusted globally. We fixed the isotropic diffuse source to be actually constant over the sky but globally refit its spectrum up to 10 GeV, since point-source fits are insensitive to diffuse emission above this energy. The Galactic diffuse emission component also was treated quite differently. Starting with a version of the Galactic diffuse model (Section 2.4.1) without its non-template diffuse γ-ray emission, we derived an alternative adjustment by optimizing the Galactic diffuse normalization for each RoI and the eight bands below 10 GeV. These values were turned into an 8-layer map, which was smoothed and then applied to the PSF-convolved diffuse model predictions for each band. Next, the corrections were remeasured. This process converged after two iterations, such that no further corrections were needed. The advantage of the procedure, compared to fitting the diffuse spectral parameters in each RoI (Section 3.2), is that the effective predictions do not vary abruptly from an RoI to its neighbors and are unique for each point. Also, it does not constrain the spectral adjustment to be a power law.

After a set of iterations had converged, the localization procedure was applied, and source positions were updated for a new set of iterations. At this stage, new sources were occasionally added using the residual TS procedure described in Section 3.1.2. The detection and localization process resulted in 7841 candidate point sources with TS > 10, of which 3179 were new. The fit validation and likelihood weighting were done as in 3FGL, except that, due to the improved representation of the Galactic diffuse, the effect of the weighting factor was less severe.

The pointlike unweighting scheme is slightly different from that described in the 3FGL paper (Section 3.1.2). A measure of the sensitivity to the Galactic diffuse component is the average count density for the RoI divided by the peak value of the PSF, N_diff, which represents a measure of the diffuse background under the point source. For the RoI at the Galactic center, and the lowest energy band, this is $4.15\times {10}^{4}$ counts. We unweight the likelihood for all energy bands by effectively limiting this implied precision to 2%, corresponding to 2500 counts. As before, we divide the log-likelihood contribution from this energy band by $\max (1,{N}_{\mathrm{diff}}/2500)$. For the aforementioned case, this value is 16.6. A consequence is to increase the spectral fit uncertainty for the lowest energy bins for every source in the RoI. The value for this unweighting factor was determined by examining the distribution of the deviations between fluxes fitted in individual energy bins and the global spectral fit (similar to what is done in Section 3.5). The 2% precision was set such that the rms for the distribution of positive deviations in the most sensitive lowest energy band was near the statistical expectation. (Negative deviations are distorted by the positivity constraint, resulting in an asymmetry of the distribution).

An important validation criterion is the all-sky counts residual map. Since the source overlaps and diffuse uncertainties are most severe at the lowest energy, we present, in Figure 3, the distribution of normalized residuals per pixel, binned with N_side = 64, in the 100–177 MeV Front energy band. There are 49,920 such pixels, with data counts varying from 92 to 1.7 × 10⁴. For $| b| \gt 10^\circ$, the agreement with the expected Gaussian distribution is very good, while it is clear that there are issues along the plane. These are of two types. First, around very strong sources, such as Vela, the discrepancies are perhaps a result of inadequacies of the simple spectral models used, but the (small) effect of energy dispersion and the limited accuracy of the IRFs may contribute. Regions along the Galactic ridge are also evident, a result of the difficulty modeling the emission precisely, the reason we unweight contributions to the likelihood.

Zoom In Zoom Out Reset image size

3.1.2. Detection of Additional Sources

As in 3FGL, the same implementation of the likelihood used for optimizing source parameters was used to test for the presence of additional point sources. This is inherently iterative, in that the likelihood is valid to the extent that the model used to calculate it is a fair representation of the data. Thus, the detection of the faintest sources depends on accurate modeling of all nearby brighter sources and the diffuse contributions.

The FL8Y source list from which this started represented several such additions from the 4 yr 3FGL. As before, an iteration starts with choosing a HEALPix N_side = 512 grid, 3.1 M points with average separation  0.15 degrees. But now, instead of testing a single power-law spectrum, we try five spectral shapes; three are power laws with different indices, two with significant curvature. Table 3 lists the spectral shapes used for the templates. They are shown in Figure 4.

Zoom In Zoom Out Reset image size

Table 3.  Spectral Shapes for Source Search

α	β	E0 (GeV)	Template	Generated	Accepted
1.7	0.0	50.00	Hard	471	101
2.2	0.0	1.00	Intermediate	889	177
2.7	0.0	0.25	Soft	476	84
2.0	0.5	2.00	Peaked	686	151
2.0	0.3	1.00	Pulsar-like	476	84

Note. The spectral parameters α, β, and E₀ refer to the LogParabola spectral shape (Equation (2)). The last two columns show the number, for each shape, that were successfully added to the pointlike model, and the number accepted for the final 4FGL list.

Download table as:  ASCIITypeset image

For each trial position, and each of the five templates, the normalizations were optimized, and the resulting TS associated with the pixel. Then, as before, but independently for each template, a cluster analysis selected groups of pixels with TS > 16, as compared to TS > 10 for 3FGL. Each cluster defined a seed, with a position determined by weighting the TS values. Finally, the five sets of potential seeds were compared and, for those within 1°, the seed with the largest TS was selected for inclusion.

Each candidate was added to its respective RoI, then fully optimized, including localization, during a full likelihood optimization including all RoIs. The combined results of two iterations of this procedure, starting from a pointlike model including only sources imported from the FL8Y source list, are summarized in Table 3, which shows the number for each template that was successfully added to the pointlike model, and the number finally included in 4FGL. The reduction is mostly due to the TS > 25 requirement in 4FGL, as applied to the gtlike calculation (Section 3.2), which uses different data and smaller weights. The selection is even stricter (TS > 34, Section 3.3) for sources with curved spectra. Several candidates at high significance were not accepted because they were too close to even brighter sources, or inside extended sources, and thus unlikely to be independent point sources.

3.1.3. Localization

The position of each source was determined by maximizing the likelihood with respect to its position only. That is, all other parameters are kept fixed. The possibility that a shifted position would affect the spectral models or positions of nearby sources is accounted for by iteration. In the ideal limit of large statistics, the log-likelihood is a quadratic form in any pair of orthogonal angular variables, assuming small angular offsets. We define Localization Test Statistic (LTS) to be twice the log of the likelihood ratio of any position with respect to the maximum; the LTS evaluated for a grid of positions is called an LTS map. We fit the distribution of LTS to a quadratic form to determine the uncertainty ellipse (position, major and minor axes, and orientation). The fitting procedure starts with a prediction of the LTS distribution from the current elliptical parameters. From this, it evaluates the LTS for eight positions in a circle of a radius corresponding to twice the geometric mean of the two Gaussian sigmas. We define a measure, the localization quality (LQ), of how well the actual LTS distribution matches this expectation as the sum of squares of differences at those eight positions. The fitting procedure determines a new set of elliptical parameters from the eight values. In the ideal case, this is a linear problem and one iteration is sufficient from any starting point. To account for finite statistics or distortions due to inadequacies of the model, we iterate until changes are small. The procedure effectively minimizes LQ.

We flagged apparently significant sources that do not have good localization fits (LQ $\gt 8$) with Flag 9 (Section 3.7.3), and for them, we estimated the position and uncertainty by performing a moment analysis of an LTS map instead of fitting a quadratic form. Some sources that did not have a well-defined peak in the likelihood were discarded by hand, on the consideration that they were most likely related to residual diffuse emission. Another possibility is that two adjacent sources produce a dumbbell-like shape; for a few of these cases, we added a new source by hand.

As in 3FGL, we checked the sources spatially associated with 984 AGN counterparts, comparing their locations with the well-measured positions of the counterparts. Better statistics allowed examination of the distributions of the differences separately for bright, dim, and moderate-brightness sources. From this, we estimate the absolute precision Δ_abs (at the 95% confidence level) more accurately at ∼00068, up from ∼0005 in 3FGL. The systematic factor f_rel was 1.06, slightly up from 1.05 in 3FGL. Equation (1) shows how the statistical errors Δ_stat are transformed into total errors Δ_tot:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mathrm{tot}}^{2}={\left({f}_{\mathrm{rel}}{{\rm{\Delta }}}_{\mathrm{stat}}\right)}^{2}+{{\rm{\Delta }}}_{\mathrm{abs}}^{2},\end{eqnarray} \tag{ 1 }$$

which is applied to both ellipse axes.

The framework for this stage of the analysis is inherited from the 3FGL catalog. It splits the sky into RoIs, varying typically half a dozen sources near the center of the RoI at the same time. Each source is entered into the fit with the spectral shape and parameters obtained by pointlike (Section 3.1), the brightest sources first. Soft sources from pointlike within 02 of bright ones were intentionally deleted. They appear because the simple spectral models we use are not sufficient to account for the spectra of bright sources, but including them would bias the spectral parameters. There are 1748 RoIs for 4FGL, listed in the ROIs extension of the catalog (Appendix A). The global best fit is reached iteratively, injecting the spectra of sources in the outer parts of the RoI from the previous step or iteration. In this approach, the diffuse emission model (Section 2.4) is taken from the global templates (including the spectrum, unlike what is done with pointlike in Section 3.1), but it is modulated in each RoI by three parameters: normalization (at 1 GeV) and small corrective slope of the Galactic component, and normalization of the isotropic component.

Among the more than 8000 seeds coming from the localization stage, we keep only sources with TS > 25, corresponding to a significance of just over 4σ evaluated from the χ² distribution with four degrees of freedom (position and spectral parameters of a power-law source; Mattox et al. 1996). The model for the current RoI is readjusted after removing each seed below threshold. The low-energy flux of the seeds below threshold (a fraction of which are real sources) can be absorbed by neighboring sources closer than the PSF radius. As in 3FGL, we manually added known LAT pulsars that could not be localized by the automatic procedure without phase selection. However, none of those reached TS > 25 in 4FGL.

We introduced a number of improvements with respect to 3FGL (by decreasing order of importance):

• 1.  
  In 3FGL, we had already noted that systematic errors due to an imperfect modeling of diffuse emission were larger than statistical errors in the Galactic plane and were at the same level over the entire sky. With twice as much exposure and an improved effective area at low energy with Pass 8, the effect now dominates. The approach adopted in 3FGL (comparing runs with different diffuse models) allowed us to characterize the effect globally and flag the worst offenders but left purely statistical errors on source parameters. In 4FGL, we introduce weights in the maximum likelihood approach (Appendix B). This allows obtaining directly (although in an approximate way) smaller TS and larger parameter errors, reflecting the level of systematic uncertainties. We estimated the relative spatial and spectral residuals in the Galactic plane where the diffuse emission is strongest. The resulting systematic level  ∼ 3% was used to compute the weights. This is by far the most important improvement, which avoids reporting many dubious soft sources.
• 2.  
  The automatic iteration procedure at the next-to-last step of the process was improved. There are now two iteration levels. In a standard iteration, the sources and source models are fixed and only the parameters are free. An RoI and all its neighbors are run again until $\mathrm{log}{ \mathcal L }$ does not change by more than 10 from the previous iteration. Around that, we introduce another iteration level (superiterations). At the first iteration of a given superiteration, we reenter all seeds and remove (one by one) those with TS < 16. We also systematically check a curved spectral shape versus a power-law fit to each source at this first iteration and keep the curved spectral shape if the fit is significantly better (Section 3.3). At the end of a superiteration, an RoI (and its neighbors) enters the next superiteration until $\mathrm{log}{ \mathcal L }$ does not change by more than 10 from the last iteration of the previous superiteration. This procedure stabilizes the spectral shapes, particularly in the Galactic plane. Seven superiterations were required to reach full convergence.
• 3.  
  The fits are now performed from 50 MeV to 1 TeV, and the overall significances (Signif_Avg) as well as the spectral parameters refer to the full band. The total energy flux, on the other hand, is still reported between 100 MeV and 100 GeV. For hard sources with photon index less than 2, integrating up to 1 TeV would result in much larger uncertainties. The same is true for soft sources with photon indices larger than 2.5 when integrating down to 50 MeV.
• 4.  
  We considered the effect of energy dispersion in the approximate way implemented in the Science Tools. The effect of energy dispersion is calculated globally for each source and applied to the whole 3D model of that source, rather than accounting for energy dispersion separately in each pixel. This approximate rescaling captures the main effect (which is only a small correction, see Section 4.2.2) at a very minor computational cost. In evaluating the likelihood function, the effects of energy dispersion were not applied to the isotropic background and the Sun/Moon components whose spectra were obtained from the data without considering energy dispersion.
• 5.  
  We used smaller RoIs at higher energy because we are interested in the core region only, which contains the sources whose parameters come from that RoI (sources in the outer parts of the RoI are entered only as background). The core region is the same for all energy intervals, and the RoI is obtained by adding a ring to that core region, whose width adapts to the PSF and therefore decreases with energy (Table 2). This does not significantly affect the result, because the outer parts of the RoI would not have been correlated to the inner sources at high energy anyway, but this saves memory and CPU time.
• 6.  
  At the last step of the fitting procedure, we tested all spectral shapes described in Section 3.3 (including log-normal for pulsars and cutoff power law for other sources), readjusting the parameters (but not the spectral shapes) of neighboring sources.

We used only binned likelihood analysis in 4FGL because unbinned mode is much more CPU intensive and does not support weights or energy dispersion. We split the data into fifteen components, selected according to PSF event type and described in Table 2. As explained in Section 2.4.4, at low energy, we kept only the event types with the best PSF. Each event type selection has its own isotropic diffuse template (because it includes residual charged-particle background, which depends on event type). A single component is used above 10 GeV to save memory and CPU time: at high energy, the background under the PSF is small, so keeping the event types separate does not markedly improve significance; it would help for localization, but this is done separately (Section 3.1.3).

A known inconsistency in acceptance exists between Pass 8 PSF event types. It is easy to see on bright sources or the entire RoI spectrum and peaks at the level of 10% between PSF0 (positive residuals, underestimated effective area) and PSF3 (negative residuals, overestimated effective area) at a few GeV. In that range, all event types were considered, so the effect on source spectra average out. Below 1 GeV, the PSF0 event type was discarded but the discrepancy is lower at low energy. We checked by comparing with preliminary corrected IRFs that the energy fluxes indeed tend to be underestimated, but by only 3%. The bias on power-law index is less than 0.01.

The spectral representation of sources largely follows what was done in 3FGL, considering three spectral models (power law, power law with subexponential cutoff, and log-normal). We changed two important aspects of how we parameterize the cutoff power law:

• 1.  
  The cutoff energy was replaced by an exponential factor (a in Equation (4)), which is allowed to be positive. This makes the simple power law a special case of the cutoff power law and allows for fitting of that model to all sources, even those with negligible curvature.
• 2.  
  We set the exponential index (b in Equation (4)) to 2/3 (instead of 1) for all pulsars that are too faint for it to be left free. This recognizes the fact that b < 1 (subexponential) in all six bright pulsars that have b free in 4FGL. Three have b ∼ 0.55 and three have b ∼ 0.75. We chose 2/3 as a simple intermediate value.

For all three spectral representations in 4FGL, the normalization (flux density K) is defined at a reference energy E₀ chosen such that the error on K is minimal. E₀ appears as Pivot_Energy in the FITS table version of the catalog (Appendix A). The 4FGL spectral forms are thus:

• 1.  
  A log-normal representation (LogParabola under SpectrumType in the FITS table) for all significantly curved spectra except pulsars, 3C 454.3 and the Small Magellanic Cloud (SMC):

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

  The parameters K, α (spectral slope at E₀), and the curvature β appear as LP_Flux_Density, LP_Index and LP_beta in the FITS table, respectively. No significantly negative β (spectrum curved upwards) was found. The maximum allowed β was set to 1 as in 3FGL. Those parameters were used for fitting because they allow for the minimization of the correlation between K and the other parameters. However, a more natural representation would use the peak energy E_peak at which the spectrum is maximum (in $\nu {F}_{\nu }$ representation)

  $$\begin{eqnarray}&&{E}_{\mathrm{peak}}={E}_{0}\,\exp \left(\displaystyle \frac{2-\alpha }{2\,\beta }\right).\end{eqnarray} \tag{ 3 }$$
• 2.  
  A subexponentially cutoff power law for all significantly curved pulsars (PLSuperExpCutoff under SpectrumType in the FITS table):

  $$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}=K{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{-{\rm{\Gamma }}}\exp \left(a\,({E}_{0}^{b}-{E}^{b}\right))\end{eqnarray} \tag{ 4 }$$

  where E₀ and E in the exponential are expressed in MeV. The parameters K, Γ (low-energy spectral slope), a (exponential factor in MeV^−b), and b (exponential index) appear as PLEC_Flux_Density, PLEC_Index, PLEC_Expfactor and PLEC_Exp_Index in the FITS table, respectively. Note that in the Science Tools that spectral shape is called PLSuperExpCutoff2 and no ${E}_{0}^{b}$ term appears in the exponential, so the error on K (Unc_PLEC_Flux_Density in the FITS table) was obtained from the covariance matrix. The minimum Γ was set to 0 (in 3FGL it was set to 0.5, but a smaller b results in a smaller Γ). No significantly negative a (spectrum curved upwards) was found.
• 3.  
  A simple power-law form (Equation (4) without the exponential term) for all sources not significantly curved. For those parameters K and Γ appear as PL_Flux_Density and PL_Index in the FITS table.

The power law is a mathematical model that is rarely sustained by astrophysical sources over as broad a band as 50 MeV to 1 TeV. All bright sources in 4FGL are actually significantly curved downwards. Another drawback of the power-law model is that it tends to exceed the data at both ends of the spectrum, where constraints are weak. It is not a worry at high energy, but at low energy (broad PSF), the collection of faint sources modeled as power laws generates an effectively diffuse excess in the model, which will make the curved sources more curved than they should be. Using a LogParabola spectral shape for all sources would be physically reasonable, but the very large correlation between sources at low energy due to the broad PSF makes that unstable.

We use the curved representation in the global model (used to fit neighboring sources) if TS_curv > 9 (3σ significance) where ${\mathrm{TS}}_{\mathrm{curv}}=2\mathrm{log}({ \mathcal L }$(curved spectrum)$/{ \mathcal L }$(power-law)). This is a step down from 3FGL or FL8Y, where the threshold was at 16, or 4σ, while preserving stability. The curvature significance is reported as LP_SigCurv or PLEC_SigCurv, replacing the former unique Signif_Curve column of 3FGL. Both values were derived from TS_curv and corrected for systematic uncertainties on the effective area following Equation (3) of 3FGL. As a result, 51 LogParabola sources (with TS_curv > 9) have LP_SigCurv less than 3.

Sources with curved spectra are considered significant whenever TS > 25 + 9 = 34. This is similar to the 3FGL criterion, which requested TS > 25 in the power-law representation, but accepts a few more strongly curved faint sources (pulsar-like).

One more pulsar (PSR J1057−5226) was fit with a free exponential index, besides the six sources modeled in this way in 3FGL. The Crab was modeled with three spectral components, as in 3FGL, but the inverse-Compton emission of the nebula (now an extended source, Section 3.4) was represented as a log-normal instead of a simple power law. The parameters of that component were fixed to α = 1.75, β = 0.08, K = 5.5 × 10⁻¹³ ph cm⁻² MeV⁻¹ s⁻¹ at 10 GeV, mimicking the broken power law fit by Buehler et al. (2012). They were unstable (too much correlation with the pulsar) without phase selection. Four extended sources had fixed parameters in 3FGL. The parameters in these sources (Vela X, MSH 15−52, γ Cygni, and the Cygnus X cocoon) were freed in 4FGL.

Overall in 4FGL, seven sources (the six brightest pulsars and 3C 454.3) were fit as PLSuperExpCutoff with free b (Equation (4)), 214 pulsars were fit as PLSuperExpCutoff with b = 2/3, the SMC was fit as PLSuperExpCutoff with b = 1, 1302 sources were fit as LogParabola (including the fixed inverse-Compton component of the Crab and 38 other extended sources), and the rest were represented as power laws. The larger fraction of curved spectra compared to 3FGL is due to the lower TS_curv threshold.

The way the parameters are reported has changed as well:

• 1.  
  The spectral shape parameters are now explicitly associated with the spectral model they come from. They are reported as Shape_Param where Shape is one of PL (PowerLaw), PLEC (PLSuperExpCutoff), or LP (LogParabola) and Param is the parameter name. Columns Shape_Index replace Spectral_Index, which was ambiguous.
• 2.  
  All sources were fit with the three spectral shapes, so all fields are filled. The curvature significance is calculated twice by comparing power law with both log-normal and exponentially cutoff power law (although only one is actually used to switch to the curved shape in the global model, depending on whether the source is a pulsar or not). There are also three Shape_Flux_Density columns referring to the same Pivot_Energy. The preferred spectral shape (reported as SpectrumType) remains what is used in the global model, when the source is part of the background (i.e., when fitting the other sources). It is also what is used to derive the fluxes, their uncertainties, and the significance.

This additional information allows us to compare unassociated sources with either pulsars or blazars using the same spectral shape. This is illustrated in Figure 5. Pulsar spectra are more curved than AGNs, and among AGNs, flat-spectrum radio quasars (FSRQs) peak at lower energy than BL Lacs (BLL). It is clear that when the error bars are small (bright sources), any of those plots is very discriminant for classifying sources. They complement the variability versus curvature plot (Figure 8 of the 1FGL paper). We expect most of the (few) bright remaining unassociated sources (black plus signs) to be pulsars, from their location on those plots. The same reasoning implies that most of the unclassified blazars (bcu) should be FSRQs, although the distinction with BL Lacs is less clear-cut than with pulsars. Unfortunately, most unassociated sources are faint (TS < 100) and for those, the same plots are very confused, because the error bars become comparable to the ranges of parameters.

Zoom In Zoom Out Reset image size

As in the 3FGL catalog, we explicitly model as spatially extended those LAT sources that have been shown in dedicated analyses to be spatially resolved by the LAT. The catalog process does not involve looking for new extended sources, testing possible extension of sources detected as point-like, nor refitting the spatial shapes of known extended sources.

Most templates are geometrical, so they are not perfect matches to the data and the source detection often finds residuals on top of extended sources, which are then converted into additional point sources. As in 3FGL, those additional point sources were intentionally deleted from the model, except if they met two of the following criteria: associated with a plausible counterpart known at other wavelengths, much harder than the extended source (Pivot_Energy larger by a factor e or more) or very significant (TS > 100). Contrary to 3FGL, that procedure was applied inside the Cygnus X cocoon as well.

The latest compilation of extended Fermi-LAT sources prior to this work consists of the 55 extended sources entered in the 3FHL catalog of sources above 10 GeV (Ajello et al. 2017). This includes the result of the systematic search for new extended sources in the Galactic plane ($| b| \lt 7^\circ$) above 10 GeV (FGES; Ackermann et al. 2017b). Two of those were not propagated to 4FGL:

• 1.  
  FGES J1800.5−2343 was replaced by the W 28 template from 3FGL, and the nearby excesses (Hanabata et al. 2014) were left to be modeled as point sources.
• 2.  
  FGES J0537.6+2751 was replaced by the radio template of S 147 used in 3FGL, which fits better than the disk used in the FGES paper (S 147 is a soft source, so it was barely detected above 10 GeV).

The supernova remnant (SNR) MSH 15-56 was replaced by two morphologically distinct components, following Devin et al. (2018): one for the SNR (SNR mask in the paper) and the other one for the pulsar wind nebula (PWN) inside it (radio template). We added back the W 30 SNR on top of FGES J1804.7−2144 (coincident with HESS J1804−216). The two overlap but the best localization clearly moves with energy from W 30 to HESS J1804−216.

Eighteen sources were added, resulting in 75 extended sources in 4FGL:

• 1.  
  The Rosette nebula and Monoceros SNR (too soft to be detected above 10 GeV) were characterized by Katagiri et al. (2016a). We used the same templates.
• 2.  
  The systematic search for extended sources outside the Galactic plane above 1 GeV (FHES, Ackermann et al. 2018) found sixteen reliable extended sources. Three of them were already known as extended sources. Two were extensions of the Cen A lobes, which appear larger in γ-rays than the WMAP template that we use following Abdo et al. (2010d). We did not consider them, waiting for a new morphological analysis of the full lobes. We ignored two others: M31 (extension only marginally significant, both in FHES and Ackermann et al. 2017a) and CalTech A (CTA) 1 (SNR G119.5+10.2) around PSR J0007+7303 (not significant without phase gating). We introduced the nine remaining FHES sources, including the inverse-Compton component of the Crab Nebula and the ρ Oph star-forming region (= FHES J1626.9−2431). One of them (FHES J1741.6−3917) was reported by Araya (2018a) as well, with similar extension.
• 3.  
  Four H.E.S.S. sources were found to be extended sources in the Fermi-LAT range as well: HESS J1534−571 (Araya 2017), HESS J1808−204 (Yeung et al. 2016), HESS J1809−193, and HESS J1813−178 (Araya 2018b).
• 4.  
  Three extended sources were discovered in the search for GeV emission from magnetars (Li et al. 2017a). They contain SNRs (Kes 73, Kes 79, and G42.8+0.6) but are much bigger than the radio SNRs. One of them (around Kes 73) was also noted by Yeung et al. (2017).

Table 4 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 distinction being that no position uncertainties are reported and their names end in e (see Appendix A). Unidentified point sources inside extended ones are indicated as "xxx field" in the ASSOC2 column of the catalog.

Thanks to the improved statistics, the source photon fluxes in 4FGL are reported in seven energy bands (1: 50–100 MeV; 2: 100–300 MeV; 3: 300 MeV–1 GeV; 4: 1–3 GeV; 5: 3–10 GeV; 6: 10–30 GeV; 7: 30–300 GeV) extending both below and above the range (100 MeV–100 GeV) covered in 3FGL. Up to 10 GeV, the data files were exactly the same as in the global fit (Table 2). To get the best sensitivity in band 6 (10–30 GeV), we split the data into four components per event type, using pixel size 004 for PSF3, 005 for PSF2, 01 for PSF1 and 02 for PSF0. Above 30 GeV (band 7), we used unbinned likelihood, which is as precise while using much smaller files. It does not allow us to correct for energy dispersion, but this is not an important issue in that band. The fluxes were obtained by freezing the power-law index to that obtained in the fit over the full range and adjusting the normalization in each spectral band. For the curved spectra (Section 3.3), the photon index in a band was set to the local spectral slope at the logarithmic mid-point of the band $\sqrt{{E}_{n}{E}_{n+1}}$, restricted to be in the interval [0,5].

In each band, the analysis was conducted in the same way as for the 3FGL catalog. To adapt more easily to new band definitions, the results (photon fluxes and uncertainties, νF_ν differential fluxes, and significances) are reported in a set of four vector columns (Appendix A: Flux_Band, Unc_Flux_Band, nuFnu_Band, Sqrt_TS_Band) instead of a set of four columns per band as in previous FGL catalogs.

The spectral fit quality is computed in a more precise way than in 3FGL from twice the sum of log-likelihood differences, as we did for the variability index (Section 3.6 of the 2FGL paper). The contribution from each band ${S}_{i}^{2}$ also accounts for systematic uncertainties on effective area via

$$\begin{eqnarray}&&{S}_{i}^{2}=\displaystyle \frac{2{\sigma }_{i}^{2}}{{\sigma }_{i}^{2}+{\left({f}_{i}^{\mathrm{rel}}{F}_{i}^{\mathrm{fit}}\right)}^{2}}\,\mathrm{log}\left[{{ \mathcal L }}_{i}({F}_{i}^{\mathrm{best}})/{{ \mathcal L }}_{i}({F}_{i}^{\mathrm{fit}})\right]\end{eqnarray} \tag{ 5 }$$

where i runs over all bands, ${F}_{i}^{\mathrm{fit}}$ is the flux predicted by the global model, ${F}_{i}^{\mathrm{best}}$ is the flux fitted to band i alone, σ_i is the statistical error (upper error if ${F}_{i}^{\mathrm{best}}\leqslant {F}_{i}^{\mathrm{fit}}$, lower error if ${F}_{i}^{\mathrm{best}}\gt {F}_{i}^{\mathrm{fit}}$) and the spectral fit quality is simply ${\sum }_{i}{S}_{i}^{2}$. The systematic uncertainties⁹¹ ${f}_{i}^{\mathrm{rel}}$ are set to 0.15 in the first band, 0.1 in the second and the last bands, and 0.05 in bands 3–6. The uncertainty is larger in the first band because only PSF3 events are used.

Too large values of spectral fit quality are flagged (Flag 10 in Table 5). Since there are seven bands and (for most sources, which are fit with the power-law model) two free parameters, the flag is set when ${\sum }_{i}{S}_{i}^{2}\gt 20.5$ (probability 10⁻³ for a χ² distribution with five degrees of freedom). Only six sources trigger this. We also set the same flag whenever any individual band is off by more than 3σ (${S}_{i}^{2}\gt 9$). This occurs in 26 sources. Among the 27 sources flagged with Flag 10 (examples in Figure 6), the Vela and Geminga pulsars are very bright sources for which our spectral representation is not good enough. A few show signs of a real second component in the spectrum, such as Cen A (H.E.S.S. Collaboration et al. 2018b). Several would be better fit by a different spectral model: the Large Magellanic Cloud (LMC) probably decreases at high energy as a power law like our own Galaxy, and 4FGL J0336.0+7502 is better fit by a PLSuperExpCutoff model. The latter is an unassociated source at 15° latitude, which has a strongly curved spectrum and is not variable: it is a good candidate for a millisecond pulsar. Other sources show deviations at low energy and are in confused regions or close to a brighter neighbor, such as the Cygnus X cocoon. This extended source contains many point sources inside it, and the PSF below 300 MeV is too broad to provide a reliable separation.

Zoom In Zoom Out Reset image size

Table 5.  Definitions of the Analysis Flags

Flaga	Nsources	Meaning
1	215	Source with TS > 35 which went to TS < 25 when changing the diffuse model (Section 3.7.1) or the analysis method (Section 3.7.2). Sources with TS ≤ 35 are not flagged with this bit because normal statistical fluctuations can push them to TS < 25.
2	215	Moved beyond its 95% error ellipse when changing the diffuse model.
3	342	Flux (>1 GeV) or energy flux (>100 MeV) changed by more than 3σ when changing the diffuse model or the analysis method. Requires also that the flux change by more than 35% (to not flag strong sources).
4	212	Source-to-background ratio less than 10% in highest band in which TS > 25. Background is integrated over $\pi {r}_{68}^{2}$ or 1 square degree, whichever is smaller.
5	398	Closer than ${\theta }_{\mathrm{ref}}$ b from a brighter neighbor.
6	92	On top of an interstellar gas clump or small-scale defect in the model of diffuse emission; equivalent to the c designator in the source name (Section 3.7.1).
7	⋯	Not used.
8	⋯	Not used.
9	136	Localization Quality > 8 in pointlike (Section 3.1) or long axis of 95% ellipse >025.
10	27	${\sum }_{i}{S}_{i}^{2}\gt 20.5$ or ${S}_{i}^{2}\gt 9$ in any band (Equation (5)).
11	⋯	Not used.
12	103	Highly curved spectrum; LP_beta fixed to 1 or PLEC_Index fixed to 0 (see Section 3.3).

Notes.

^aIn the FITS version (see Appendix A), the values are encoded as individual bits in a single column, with Flag n having value ${2}^{(n-1)}$. ^bθ_ref is defined in the highest band in which source TS > 25, or the band with highest TS if all are < 25. θ_ref is set to 377 below 100 MeV, 168 between 100 and 300 MeV (FWHM), 103 between 300 MeV and 1 GeV, 076 between 1 and 3 GeV (in-between FWHM and 2 r₆₈), 049 between 3 and 10 GeV and 025 above 10 GeV (2 r₆₈).

Download table as:  ASCIITypeset image

The fluxes in the 50–100 MeV band are very hard to estimate because of the enormous confusion. The average distance between sources (17) is about equal to the half width at half maximum of PSF3 events in that band, so it is nearly always possible to set a source to 0 and compensate by a suitable combination of flux adjustments in its neighbors. This is why only 34 sources have TS > 25 in that band (all are bright sources with global TS > 700). This is far fewer than the 198 low-energy (30–100 MeV) Fermi-LAT sources reported by Principe et al. (2018, 1FLE). The reason for this is that in 4FGL, we consider that even faint sources in the catalog can have strong low-energy emission, so the total source flux is distributed over 5000 sources, whereas 1FLE focused on finding individual peaks.

At the other extreme, 618 sources have TS > 25 in the 30 to 300 GeV band, which is entirely limited by photon counting (TS > 25 in that band corresponds to about five events). Only 13 of those are not associated with a 3FHL or FHES source. The brightest of them (at TS = 54 in that band) is a hard source associated with 1RXS J224123.5+294244, mostly significant in the last year, after the 3FHL time range.

As in past FGL catalogs, the photon fluxes between 1 and 100 GeV as well as the energy fluxes between 100 MeV and 100 GeV were derived from the full-band analysis assuming the best spectral shape, and their uncertainties come from the covariance matrix. Even though the full analysis is carried out down to 50 MeV and up to 1 TeV in 4FGL, we have not changed the energy range over which we quote fluxes so that they can be easily compared with fluxes in past catalogs. The photon fluxes above 100 GeV are negligible except in the very hardest power-law sources, and the energy fluxes below 100 MeV and above 100 GeV are not precisely measured (even for soft and hard sources, respectively).

3.6.1. One-year Intervals

We started by computing light curves over 1 yr intervals. This is much faster and more stable than fitting smaller time intervals, and provides a good variability assessment already. We used binned likelihood and the same data as in the main run up to 10 GeV (Table 2), but to save disk space and CPU time, we merged event types together. Above 10 GeV, we used unbinned likelihood (more efficient when there are few events). We ignored events above 100 GeV (unimportant for variability).

As in 3FGL, the fluxes in each interval were obtained by freezing the spectral parameters to those obtained in the fit over the full range and adjusting the normalization. As in previous FGL catalogs, the fluxes in each interval are reported as photon fluxes between 0.1 and 100 GeV.

The weights appropriate for one year were computed using the procedure explained in Appendix B, entering the same data cube divided by eight (we use the same weights in each year), and ignoring the last steps specific to splitting event types. The weights are of course much larger than those for 8 yr, but remain a significant correction (the weights are less than 0.2 in the Galactic Ridge up to 300 MeV). We used the same Sun/Moon model for each year. This amounts to neglecting the modulation of their intrinsic flux along the 11 yr solar cycle.

Because of the different weights between the full analysis and that in 1 yr intervals, the average flux from the light curve F_av can differ somewhat from the flux in the total analysis F_glob (low energies are less attenuated in the analysis over 1 yr intervals). This is illustrated in Figure 7. In 4FGL, we compute the variability index TS_var (reported as Variability_Index in the FITS file) as

$$\begin{eqnarray}\begin{array}{rcl}{{TS}}_{\mathrm{var}} & = & 2\sum _{i}\mathrm{log}\left[\displaystyle \frac{{{ \mathcal L }}_{i}({F}_{i})}{{{ \mathcal L }}_{i}({F}_{\mathrm{glob}})}\right]\\ & & \qquad -\max \left({\chi }^{2}({F}_{\mathrm{glob}})-{\chi }^{2}({F}_{\mathrm{av}}),0\right)\end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\chi }^{2}(F)=\sum _{i}\displaystyle \frac{{\left({F}_{i}-F\right)}^{2}}{{\sigma }_{i}^{2}}\end{eqnarray} \tag{ 7 }$$

where F_i are the individual flux values, ${{ \mathcal L }}_{i}(F)$ is the likelihood in the interval i assuming flux F, and ${\sigma }_{i}$ are the errors on F_i (upper error if F_i ≤ F, lower error if F_i > F). The first term in Equation (6) is the same as Equation (4) of 2FGL. The second term corrects (in the Gaussian limit) for the difference between F_glob and F_av (since the average flux is known only at the very end, it could not be entered when computing ${{ \mathcal L }}_{i}(F)$). We subtract the second term only when it is positive (it is not necessarily positive because the best χ² is reached at the average weighted by ${\sigma }_{i}^{-2}$, not the straight average). On the other hand, we did not correct the variability index for the relative systematic error, which is already accounted for in the weighting procedure.

Zoom In Zoom Out Reset image size

The distribution of observed TS_var is shown in Figure 8. It looks like a composite of a power-law distribution and a χ²(7) distribution with N_int − 1 = 7 degrees of freedom, where N_int is the number of intervals. The left branch corresponds both to constant sources (such as most pulsars) and sources too faint to have measurable variability. There are many blazars among them, which are most likely just as variable as brighter blazars. This contribution of real variability to TS_var is the reason why the histogram is a little offset to the right of the χ²(7) distribution (that offset is absent in the Galactic plane, and stronger off the plane).

Zoom In Zoom Out Reset image size

Variability is considered probable when TS_var > 18.48, corresponding to 99% confidence in a χ²(7) distribution. We find 1327 variable sources with that criterion. After the χ²-based correction of Equation (6), Vela X remains below that threshold. One extended source still exceeds the variability threshold. This is HESS J1420−607 (Figure 9), confused with its parent pulsar PSR J1420−6048. A similar flux transfer occurred in the third year between the Crab pulsar and the Crab Nebula. This can be understood because the synchrotron emission of the nebula becomes much harder during flares, while our pipeline assumes the soft power-law fit over the full interval applies throughout. None of those variabilities are real.

Zoom In Zoom Out Reset image size

Besides the Crab and the known variable pulsars PSR J1227−4853 (Johnson et al. 2015) and PSR J2021+4026 (Allafort et al. 2013), three other pulsars are above the variability threshold. Two are just above it and could be chance occurrences (there are more than 200 pulsars, so we expect two above the 1% threshold). The last one is PSR J2043+2740 (Figure 10), which looks like a case of genuine variability (secular flux decrease by a factor of three).

Zoom In Zoom Out Reset image size

In 4FGL, we report the fractional variability of the sources in the FITS file as Frac_Variability. It is defined for each source from the excess variance on top of the statistical and systematic fluctuations:

$$\begin{eqnarray}&&\mathrm{Var}=\displaystyle \frac{1}{{N}_{\mathrm{int}}-1}\,\sum _{i}{\left({F}_{i}-{F}_{\mathrm{av}}\right)}^{2}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\delta F=\sqrt{\max \left(\mathrm{Var}-\displaystyle \frac{{\sum }_{i}{\sigma }_{i}^{2}}{{N}_{\mathrm{int}}},0\right)}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\sigma }_{F}}{F}=\max \left(\displaystyle \frac{1}{\sqrt{2({N}_{\mathrm{int}}-1)}}\displaystyle \frac{{V}_{i}}{{F}_{\mathrm{av}}\,\delta F},10\right)\end{eqnarray} \tag{ 10 }$$

where the fractional variability itself is simply δF/F_av. This is similar to Equation (3) of 1FGL, except we omit the systematic error term because it is now incorporated in the ${\sigma }_{i}^{2}$ via the weights. The error σ_F/F is estimated from the expected scatter on the sample variance Var, which is the dominant source of uncertainty. We cap it at 10 to avoid reporting meaningless high uncertainties. Figure 11 can be compared to Figure 8 of Abdo et al. (2009c), which was based on one-week intervals (and contained many fewer sources, of course). The fractional variability is similar in the two figures, going up to 1, reflecting the absence of a preferred variability timescale in blazars. The criterion we use is not sensitive to relative variations smaller than 50% at TS = 100, so only bright sources can populate the lower part of the plot. There is no indication that fainter sources are less variable than brighter ones, but we simply cannot measure their variability.

Zoom In Zoom Out Reset image size

3.6.2. Two-month Intervals

To characterize variability, it is of course useful to have information on shorter timescales than one year. Rather than use monthly bins as in 3FGL (which would have resulted in many upper limits), we have chosen to keep the same number of intervals and build light curves over 48 two-month bins. Because the analysis is not limited by systematics at low energy over two months, we tried to optimize the data selection differently. We used binned likelihood up to 3 GeV and the same zenith angle cuts as in Table 2, but included PSF2 events between 50 and 100 MeV (not only PSF3), and added PSF1 events between 100 and 300 MeV to our standard PSF2+3 selection. This improves the average source significance over one bin, and the Earth limb contamination remains minor. Similarly to the one-year analyses, to save disk space and CPU time we merged event types together in the binned data sets. We used unbinned likelihood above 3 GeV and again ignored events above 100 GeV (unimportant for variability).

The weights appropriate for two months were computed using the same procedure (Appendix B), entering the total data cube divided by 48 (same weights in each interval). The weights are of course larger than those for one year but remain a significant correction in the Galactic plane. Up to 100 MeV the weights range from 0.2 in the Galactic Ridge to 0.85 at high latitude. At 300 MeV, they increase to 0.55 in the Galactic Ridge and 0.99 at high latitude. We used a different Sun/Moon model for each interval (the Sun averages out only over one year) but, again, assuming constant flux.

Variability is considered probable when TS_var > 72.44, corresponding to 99% confidence in a χ² distribution with N_int − 1 = 47 degrees of freedom. We find 1173 variable sources with that criterion, 1057 of which were also considered variable with one-year intervals. Among the 116 sources considered variable only with two-month light curves, 37 (1% of 3738) would be expected by chance, so more than two-thirds must be truly variable. Similarly, 270 sources are considered variable only with one-year intervals (39 expected by chance).

Two extended sources exceed the two-month variability threshold. They are the Monoceros SNR and the Cen A lobes. Both are very extended (several degrees). It is likely that their variability is due to a flaring background source that was missed by the global source detection over eight years. Indeed, the peak in the light curve of the Monoceros SNR is in 2012 June–July, at the time of Nova V959 Mon 2012 (Ackermann et al. 2014a). Another unexpected variable source is the Geminga pulsar. We think that its variability is not real but due to the direct pointings triggered toward the Crab when it was flaring (Geminga is 15° away), combined with details of the effective area or PSF dependence on off-axis angle and azimuthal, that normally average out in scanning mode.

Because the source fluxes are not allowed to be negative, the distribution of fluxes for a given source is truncated at 0. For faint sources, this results in a slight overestimate of the average flux (of no consequence) but also an underestimate of the sample variance (Equation (8)). As a result, the fractional variability (Equation (9)) is underestimated for faint sources and is often zero for weakly variable sources (below threshold). This even happens for two sources considered variable (just above threshold).

More sources are found to be variable using one-year intervals than using two-month intervals. The reason is illustrated in Figure 12, which shows the variability indices divided by N_int − 1 (so that they become directly comparable). If the sources behaved like white noise (as the statistical errors) then the correlation would be expected to follow the diagonal. But blazars behave as red noise (more variability on longer timescales) so the correlation is shifted to the right, and it is more advantageous to use longer intervals to detect variability with that criterion, because statistical errors decrease more than intrinsic variability.

Zoom In Zoom Out Reset image size

Extending this relation to even shorter intervals, the 2FAV catalog of Fermi-LAT flaring sources (Abdollahi et al. 2017), which used one-week intervals, found 518 significantly varying sources. The methodology was completely different (it did not start from a catalog over many years), and the duration was a little shorter (7.4 yr). However, the same trend remains to find fewer variable sources on shorter intervals. Not all sources are dominated by red noise though, and a fraction are above the diagonal in Figure 12. An example is provided in Figure 13 (left panel). In all cases, the variability is of course much better characterized with smaller intervals. An extreme example is provided in Figure 13 (right panel).

Zoom In Zoom Out Reset image size

3.7.1. Diffuse Emission Model

3.7.2. Analysis Method

3.7.3. Analysis Flags

The catalog is available online,⁹² together with associated products. It contains 5064 sources.⁹³ The source designation is 4FGL JHHMM.m+DDMM where the 4 indicates that this is the fourth LAT catalog, and FGL represents Fermi Gamma-ray LAT. Sources confused with interstellar cloud complexes are singled out by a c appended to their names, where the c indicates that caution should be used in interpreting or analyzing these sources. The 75 sources that were modeled as extended for 4FGL (Section 3.4) are singled out by an e appended to their names. The catalog columns are described in Appendix A. Figure 14 illustrates the distribution of the 4FGL sources over the sky, separately for AGN (blue) and other (red) classes.

Zoom In Zoom Out Reset image size

4.2.1. General Comparison

Figure 15 shows the energy flux distribution in 1FGL, 2FGL, 3FGL, and 4FGL outside the Galactic plane. Comparing the current flux threshold with those published in previous LAT Catalog papers, we see that in 4FGL, the threshold is down to ≃2 × 10⁻¹² erg cm⁻² s⁻¹. This is about a factor of two better than 3FGL. In the background-limited regime (up to a few GeV), doubling the exposure time would lead only to a factor $\sqrt{2}$. The remaining factor is due to the increased acceptance, the better PSF, and splitting the data into the PSF event types (Section 2.2). The weights (Appendix B) do not limit the general detection at high latitudes. Above 10⁻¹¹ erg cm⁻² s⁻¹, the 2FGL and 3FGL distributions are entirely compatible with 4FGL. The 1FGL distribution shows a distinct bump between 1 and 2 × 10⁻¹¹ erg cm⁻² s⁻¹. That accumulation of fluxes was clearly incorrect. We attribute it primarily to overestimating significances and fluxes due to the unbinned likelihood bias in the 1FGL analysis, and also to the less accurate procedure then used to extract source flux (see discussion in the 2FGL paper).

Zoom In Zoom Out Reset image size

The threshold at low flux is less sharp in 4FGL than it was in 2FGL or 3FGL. This reflects a larger dependence of the detection threshold on the power-law index (Figure 16). The expected detection threshold is computed from Equation A1 of Abdo et al. (2010e). The systematic limitation (entered in the weighted log-likelihood as described in Appendix B) is accounted for approximately by limiting the integral over angles to θ_max(E) such that g(θ_max, E) = , since $g({\theta }_{\max },E)$ in that equation is exactly the source to background ratio. The detection threshold for soft sources decreases only slowly with exposure due to that. On the other hand, the detection threshold improves in nearly inverse proportion to exposure for hard sources because energies above 10 GeV are still photon-limited (not background-limited).

Zoom In Zoom Out Reset image size

The power-law index Γ provides a way to compare all sources over all catalog generations, ignoring the complexities of the curved models. Figure 17 shows that the four distributions of the power-law indices of the sources at high Galactic latitude are very similar. Their averages and dispersions are ${{\rm{\Gamma }}}_{1\mathrm{FGL}}=2.22\pm 0.33$ ${{\rm{\Gamma }}}_{2\mathrm{FGL}}=2.17\pm 0.30$, ${{\rm{\Gamma }}}_{3\mathrm{FGL}}=2.22\,\pm 0.31$ and ${{\rm{\Gamma }}}_{4\mathrm{FGL}}=2.23\pm 0.30$.

Zoom In Zoom Out Reset image size

Small differences in the power-law index distributions could be related to slightly different systematic uncertainties in the effective area between the IRFs used, respectively, for 4FGL, 3FGL, 2FGL, and 1FGL (Table 1). There is actually no reason why the distribution should remain the same, since the detection threshold depends on the index and the log N–log S of FSRQs, which are soft Fermi-LAT sources, differs from that of BL Lacs, whose spectra are hard in the LAT band (Ackermann et al. 2015, Figure 7). The apparent constancy may largely be the result of competing effects.

We have compared the distribution of error radii (defined as the geometric mean of the semimajor and semiminor axes of the 95% confidence error ellipse) of the 1FGL, 2FGL, 3FGL, and 4FGL sources at high Galactic latitude. Overall, the source localization improves with time as more photons are added to previously detected sources. We concentrate instead on what happens specifically for faint sources. Figure 18 shows the distribution of 95% confidence error radii for those sources with 25 < TS < 100 in any of the catalogs. The improvement at a given TS level is partly due to the event-level analysis (from Pass 6 to 7 and 8, see Table 1) and partly to the fact that, at a given significance level and for a given spectrum, fainter sources over longer exposures are detected with more photons. This improvement is key to preserving a high rate of source associations (Section 6), even though the source density increases.

Zoom In Zoom Out Reset image size

4.2.2. Step-by-step from 3FGL to 4FGL

6.1.1. Active Galactic Nuclei

6.1.2. Other Galaxies

The Galactic sources include:

• 1.  
  239 pulsars (PSR). The public list of LAT-detected pulsars is regularly updated.⁹⁷ Some 232 pulsars in this list are included in 4FGL (68 would have been missed by the association pipeline using the ATNF catalog), while six are absent because they did not pass the TS > 25 criterion. These pulsars represent by far the largest population of identified sources in 4FGL. Another seven pulsars from the ATNF database are associated with 4FGL sources with high-confidence via the Bayesian method that we consider γ-ray pulsar candidates. This method predicts that about 30 other 4FGL sources are ATNF pulsars. Note that out of the 24 pulsar candidates presented in 3FGL, pulsations have now been detected for 19 of them. The other five are not associated with pulsars in 4FGL.
• 2.  
  40 SNRs. Out of them, 24 are extended and thus firmly identified. The other 16 are not resolved. SNR G150.3+4.5 has a log-normal spectral shape with a very hard photon index Γ of 1.6, which indicates that the emission is most likely leptonic and makes this source an excellent candidate for the Cerenkov Telescope Array (CTA). In contrast, the softer spectrum of the LMC SNR N 132D (photon index = 2.07) makes the hypothesis of a dominant hadronic emission likely. The significant spectral curvature seen in Puppis A is consistent with its non-detection in the TeV domain.
• 3.  
  17 PWNe, 15 of them being extended. New associations are N 157B, PWN G63.7+1.1, HESS J1356−645, FGES J1631.6−4756, FGES J1836.5−0651, FGES J1838.9−0704, HESS J1857+026. The median photon index of the 4FGL PWNe is 2.31. N 157B, located in the LMC, has a photon index of 2.0, hinting at an additional contribution from a (yet-undetected) pulsar at low energy on top of the PWN.
• 4.  
  78 unassociated sources overlapping with known PWNe or SNRs (SPP). Estimation of missed associations of SNR, PWN and SPP sources is made difficult by the intrinsic spatial extension of the sources; no attempts have thus been made along this line.
• 5.  
  30 globular clusters (GLC). Missing relative to 3FGL is 2MS−GC01. The 16 new associations are NGC 362, NGC 1904, NGC 5286, NGC 5904, NGC 6139, NGC 6218, NGC 6304, NGC 6341, NGC 6397, NGC 6402, NGC 6838, NGC 7078, Terzan 1, Terzan 2, GLIMPSE C01, GLIMPSE C02. Only two other 4FGL sources are estimated to be GLCs.
• 6.  
  Six high-mass X-ray binaries (HMB). The three new sources are HESS J0632+057, which has a reported LAT detection after 3FGL (Li et al. 2017b), Cyg X-1, an archetypical black hole binary reported after the 3FGL (Zanin et al. 2016; Zdziarski et al. 2017), and RX J0648.0−4418/HD 49798, which is a peculiar X-ray binary (Mereghetti et al. 2011; Popov et al. 2018). The association probability of RX J0648.0−4418/HD 49798 is just barely larger (0.85 versus 0.80) than that of the blazar candidate SUMSS J064744−441946. Three other 4FGL sources are estimated to be HMBs according to the Bayesian method.
• 7.  
  Three star-forming regions; new since 3FHL is the association of the extended source FHES J1626.9−2431 (Section 3.4) with the ρ Ophiuchi star-forming region. Positional coincidences between 4FGL sources and two of the brightest extended H ii regions present in the catalog of Paladini et al. (2003) have been found. They are reported here as candidate associations: one region corresponds to NGC 6618 in M17, whose extension of 6' at 2.7 GHz encompasses 4FGL J1820.4−1609; the second one corresponds to NGC 4603, which has a similar extension of 6' at 2.7 GHz and encompasses 4FGL J1115.1−6118.
• 8.  
  Two low-mass X-ray binaries (LMB). PSR J1023+0038 is a known binary millisecond pulsar/LMB transition system, with a change in γ-ray flux detected (Stappers et al. 2014) simultaneously with a state change, and was previously detected as 2FGL J1023.6+0040 (but not detected in 3FGL). The LMB 2S 0921−630 (V395 Car) is a well-studied binary involving a neutron star and a K0 III star with an orbital period of 9 days (Shahbaz & Watson 2007) and is a new LAT detection.
• 9.  
  One binary star system (BIN), η Carinae (Abdo et al. 2010f; Reitberger et al. 2015).
• 10.  
  One nova (NOV), V5668 Sagittarii (Cheung et al. 2016). Other novae detected by the LAT are missing. Novae have short durations, and most are below the significance threshold because their signal is diluted over the eight years of 4FGL data. As discussed in Section 3.6.2, Nova V959 Mon 2012 is confused with the SNR Monoceros.

As a new feature relative to previous catalogs, the most probable counterpart to a 4FGL unassociated source is given in a separate column of the FITS table, along with the corresponding association probability (applying a threshold of 0.1 on that probability). This additional information, to be used with care given its low confidence, is meant to foster further investigations regarding the nature of these 4FGL sources and to help clarify why detections claimed in other works are sometimes missing in 4FGL. We report 124 low-confidence (0.1 < P < 0.8) associations with the Bayesian method. Note that the relative distances between γ-ray and counterpart sources remain quite small (53 are within r₉₅ and all within 1.85 r₉₅). This notably small number of low-confidence associations illustrates how quickly the Bayesian association probability drops with increasing relative distance in the case of 4FGL. Except for rare exceptions, the other 1199 4FGL sources (having not even low-confidence associations) will not get associated with any of the tested sources (i.e., belonging to the catalogs listed in Table 6) in a future LAT catalog. We also report 42 matches (classified as UNK) with sources from the Planck surveys (with 0.1 < P ≤ 1) to guide future investigations.

Out of the 1336 sources unassociated in 4FGL, 368 already present in 3FGL had no associations there. Another 27 sources previously associated in 3FGL have now lost their associations because of a shift in their locations relative to 3FGL.

About half of the unassociated sources are located less than 10° away from the Galactic plane. Their wide latitude extension is hard to reconcile with those of known classes of Galactic γ-ray sources. For instance, Figure 20 compares this latitude distribution with that of LAT pulsars. In addition to nearby millisecond pulsars, which have a quasi-isotropic distribution, the LAT detects only young isolated pulsars (age <10⁶ yr), which are by nature clustered close to the plane. Older pulsars, which have had time to drift further off the plane, show a wider Galactic-latitude distribution, more compatible with the observed distribution of the unassociated sources, but these pulsars have crossed the "γ-ray death line" (see Abdo et al. 2013) and are hence undetectable. Attempts to spatially cross correlate the unassociated population with other potential classes, e.g., LMBs (Liu et al. 2007), O stars,⁹⁸ and Be stars⁹⁹ have been unsuccessful. The observed clustering of these unassociated sources in high-density "hot spots" may be a clue that they actually correspond to yet-to-be identified, relatively nearby extended sources. The Galactic latitude distribution near the plane is clearly non-Gaussian as visible in Figure 20, which may indicate the presence of several components.

Zoom In Zoom Out Reset image size

The spectral properties of these sources can also provide insight into their nature, as illustrated in Figure 21, which shows the latitude distribution of their spectral indices. The change in spectral hardness with sky location demonstrates the composite nature of the unassociated population. The high-latitude sources have an average photon index compatible with that of blazars of unknown type (Γ = 2.24), a hint that these sources could be predominantly blazars. Unassociated sources lying closer to the Galactic plane have softer spectra, closer to that expected for young pulsars (Γ = 2.42). Another interesting possibility is that some of these unassociated sources actually correspond to WIMP dark matter annihilating in Galactic subhalos (Ackermann et al. 2012d; Coronado-Blázquez et al. 2019). Indeed, ΛCDM cosmology predicts the existence of thousands of subhalos below $\sim {10}^{7}{M}_{\odot }$, i.e., not massive enough to retain gas or stars at all. As a result, they are not expected to emit at other wavelengths, and therefore, they would not possess astrophysical counterparts. Annihilation of particle dark matter may yield a pulsar-like spectrum (Baltz et al. 2007).

Zoom In Zoom Out Reset image size

The correspondence of 4FGL sources with previous Fermi-LAT catalogs (reported in the ASSOC_FGL and ASSOC_FHL columns) was based, as in 3FGL, on error-circle overlap at the 95% confidence level, amounting to

$$\begin{eqnarray}&&{\rm{\Delta }}\leqslant {d}_{x,a}=\sqrt{{\theta }_{x,a}^{2}+{\theta }_{x,4\mathrm{FGL}}^{2}}\end{eqnarray} \tag{ 11 }$$

where Δ is the angular distance between a 4FGL source and a source in catalog a, and the θ_x are derived from the Conf_95_SemiMajor columns in the two catalogs at the x% confidence level (assuming a 2D normal distribution). We also considered that a previous LAT source corresponds to a 4FGL source whenever they have the same association (the associations can have offsets greater than θ₉₅, depending on the density of sources in the catalogs of counterparts at other wavelengths).

We checked all sources that did not have an obvious counterpart in 4FGL inside d₉₅, nor a common association. The procedure is described in detail in Section 4.2.3 of the 3FGL paper. The result is provided in one FITS file per catalog,¹⁰⁰ reporting the same information as Table 11 of the 3FGL paper: counterparts up to 1°, whether they are inside d_99.9 (= $1.52\,{d}_{95}$) or not, and specific conditions (flagged, c source, close to an extended source, split into several sources). The number of missed sources and their nature are provided in Table 10.

We have looked at the most recent catalogs, 3FGL and 3FHL, in more detail. Because the first four years are in common, we expect the 3FGL and 4FGL positions to be correlated. That correlation is however less than one might think because the data have changed (from Pass 7 to Pass 8, Section 2.2). We found that the distribution of ${\rm{\Delta }}/{d}_{\mathrm{95,3}\mathrm{FGL}}$ (when it is less than 1) is narrower by a factor 0.83 than the Rayleigh distribution. This means that, by cutting at ${d}_{\mathrm{95,3}\mathrm{FGL}}$, we expect only 1.3% misses by chance (about 40 sources). With 3FHL, the correlation is larger because it used Pass 8 already, the overlap is 7 yr, and for the hard sources present in 3FHL, the lower-energy photons do not contribute markedly to the localization. The distribution of ${\rm{\Delta }}/{d}_{\mathrm{95,3}\mathrm{FHL}}$ is narrowed by a factor 0.62, and the number of chance misses by cutting at ${d}_{\mathrm{95,3}\mathrm{FHL}}$ should be only 0.04% (about 1 source). The correlation is similarly large with 2FHL (6 yr of Pass 8 data). That correlation effect is less for earlier catalogs, so for them, the fraction of true counterparts that are found outside the combined 95% error circle is closer to 5%. Most of those true sources are expected to have a 4FGL counterpart at the 99.9% level in the FITS files.

Out of 3033 3FGL sources, 469 are missing in 4FGL for various reasons, including the change of diffuse emission model, point sources being absorbed into new extended ones, or variability effects. Most of these missing sources had low significance in 3FGL. Only 72 sources were associated. The majority are blazars (35 BCUs, 17 FSRQs, one BLL, and one SSRQ) plus one AGN. While BLLs are 36% more numerous relative to FSRQs in 3FGL, only one has gone away in 4FGL, an effect possibly related to the larger variability of FSRQs relative to BLLs observed in the LAT energy band (Ackermann et al. 2015). Other missing sources include 11 SPPs, 3 PSRs, one SNR, and one PWN. The nova V407 Cyg is now missing as it no longer fulfills the average-significance criterion.

Two LAT pulsars are considered lost. PSR J1513−5908 (= 3FGL J1513.9−5908) inside the PWN MSH 15−52 is a pulsar peaking at MeV energies (Kuiper et al. 1999), very soft in the LAT band (Pellizzoni et al. 2009; Abdo et al. 2010c), which has gone below threshold after applying the weights. PSR J1112−6103 (= 3FGL J1111.9−6058) was split into two 4FGL sources. One is still associated with the pulsar, but it is not the one closest to the 3FGL position. The third missing pulsar association was between 3FGL J1632.4−4820 and the non-LAT PSR J1632−4818, in a confused region now covered by the extended source 4FGL J1633.0−4746e. Among the five most significant lost 3FGL sources (>20σ), the brightest one (3FGL J1714.5−3832 = CTB 37A) was split into two 4FGL sources, the brighter of which is associated instead with the newly discovered pulsar PSR J1714−3830 (Saz Parkinson et al. 2018) inside the CTB 37A SNR and, hence, was not recognized as a common association. Two others (3FGL J1906.6+0720 and 3FGL J0536.4−3347) were also split, and now both members of each pair are associated. This is definitely an improvement. The last two (3FGL J1745.3−2903c and 3FGL J1747.0−2828) were within 06 of the Galactic center, a region of the sky where changing the diffuse model had a strong impact. They have no 4FGL counterpart at all.

Concerning sources missing from 3FHL, established with Pass 8 data as 4FGL, they amount to 33, with 17 unassociated, nine blazars (four BLLs and five BCUs), one AGN, one SNR, four UNKs, and the transient HMB PSR B1259−63 (diluted over 8 yr). All of these sources had a TS close to the TS = 25 significance threshold.

The synergy between the LAT and the Cerenkov telescopes operating in the TeV energy domain has proven extremely fruitful, in particular by bringing out promising TeV candidates in the LAT catalogs. This approach, further motivated by the upcoming deployment of the CTA, has justified the release of LAT source catalogs above 10 GeV, like the 3FHL (Ajello et al. 2017) based on 7 yr of data. The associations of 4FGL sources with extended sources listed in TeVCat¹⁰¹ are presented in Table 11. Relative to 3FHL, nine new extended TeV sources are associated with 4FGL extended sources (TeV sources: HESS J1534−571, HESS J1808−204, HESS J1809−193, see Section 3.4), or (sometimes multiple) 4FGL point sources (TeV sources: HESS J1718−385, HESS J1729−345, HESS J1848−018, HESS J1858+020, MGRO J1908+06, HESS J1912+101). All TeV blazars have 4FGL counterparts. The median value of Γ for 4FGL point sources associated with TeV point sources is 1.95, indicating hard spectra as expected. In associations with extended TeV sources, the median Γ changes from 2.09 to 2.38 depending on whether the 4FGL sources are extended or not. This fairly large difference favors the interpretation that most associations between extended TeV sources and non-extended 4FGL sources are accidental.

Whenever a high-confidence association with a point-like counterpart is obtained, we provide the most accurate counterpart position available and its uncertainty. In particular, 2775 4FGL AGNs have Very Long Baseline Interferometry (VLBI) counterparts. VLBI, i.e., radio interferometry with baseline lengths of >1000 km, is sensitive to radio emission from compact regions of AGNs that are smaller than 20 milliarcsecond (mas), which corresponds to parsec scales. Such observations allow for the determination of positions of the AGN jet base with mas-level accuracy. We used the RFC catalog based on the dedicated ongoing observing program (Schinzel et al. 2015, 2017) with the Very Long Baseline Array (Napier et al. 1994), as well as VLBI data under other programs. The association between γ-ray source and VLBI counterpart was evaluated along a similar, but distinct, scheme as that presented in Section 5. This scheme (see Petrov et al. 2013, for more details) is based on the strong connection between the γ-ray emission and radio emission at parsec scales and on the sky density of bright compact radio sources being relatively low. The chance of finding a bright background, unrelated compact radio source within the LAT positional error ellipse is low enough to establish an association. The likelihood ratio (with a somewhat different definition from that implemented in the LR-method) was required to be greater than 8 to claim an association, with an estimated false-association fraction of 1%.

For AGNs without VLBI counterparts, the position uncertainties were set to typical values of 20'' for sources associated from the RASS survey and 10'' otherwise. For identified pulsars, the position uncertainties come from the rotation ephemeris used to find γ-ray pulsations, many of which were obtained from radio observations (Smith et al. 2019). If the ephemeris does not include the uncertainties and for pulsar candidates, we use the ATNF psrcat values. If neither of those exist, we use the 01 uncertainties from the list maintained by the West Virginia University Astrophysics group.¹⁰² Ephemeris position uncertainties are often underestimated, so we arbitrarily apply a minimum uncertainty of 1 mas. For GLC, from Harris (1996),¹⁰³ the position uncertainties were assigned a typical value of 2''.