Characterizing the population of pulsars in the inner Galaxy with the Fermi Large Area Telescope

TL;DR

This study identifies approximately 135 potential gamma-ray pulsar candidates in the inner Galaxy using Fermi LAT data, supporting the hypothesis that unresolved pulsars contribute to the Galactic Center gamma-ray excess.

Contribution

The paper develops new selection criteria to identify gamma-ray pulsar candidates and applies them to Fermi LAT data using two different interstellar emission models, revealing a significant candidate population.

Findings

Approximately 135 pulsar candidates identified.

66 candidates are consistent across different emission models.

Supports pulsar population as a source of the Galactic Center gamma-ray excess.

Abstract

An excess of $\gamma$-ray emission from the Galactic Center (GC) region with respect to predictions based on a variety of interstellar emission models and $\gamma$-ray source catalogs has been found by many groups using data from the {\it Fermi} Large Area Telescope (LAT). Several interpretations of this excess have been invoked. In this paper we search for members of an unresolved population of $\gamma$-ray pulsars located in the inner Galaxy that are predicted by the interpretation of the GC excess as being due to a population of such sources. We use cataloged LAT sources to derive criteria that efficiently select pulsars with very small contamination from blazars. We search for point sources in the inner $40^\circ\times40^\circ$ region of the Galaxy, derive a list of approximately 400 sources, and apply pulsar selection criteria to extract pulsar candidates among our source list. We performed the entire data analysis chain with two different interstellar emission models (IEMs), and found a total of 135 pulsar candidates, of which 66 were selected with both IEMs.

Figures

Tables

Table 1. Table 1 Spectral parameters of PSR candidates compared with known PSRs.

IEM	$N_{\mathrm{PSR}}$	$\Gamma$	${\log }_{10}\left(E_{\mathrm{cut}}\left[\mathrm{MeV}\right]\right)$
Off.	86	$1.03\pm 0.52$	$3.28\pm 0.33$
Alt.	115	$1.05\pm 0.50$	$3.27\pm 0.31$
Alt. $\cap$ Off. (Off.)	66	$1.02\pm 0.52$	$3.27\pm 0.32$
Alt. $\cap$ Off. (Alt.)	66	$1.01\pm 0.51$	$3.26\pm 0.30$
Known PSRs (Off.)	172	$1.33\pm 0.54$	$3.43\pm 0.24$
Young PSRs (Off.)	86	$1.46\pm 0.53$	$3.44\pm 0.26$
MSPs(Off.)	86	$1.20\pm 0.50$	$3.42\pm 0.23$

Table 2. Table 2 Contents of the 2FIG source list FITS table.

Contents	Column Name	Units	Uncertainty
Source designation	Source_Name	$\cdots$	$\cdots$
Right ascension	RAJ2000	[deg]	$\cdots$
Declination	DEJ2000	[deg]	$\cdots$
Galactic longitude	GLON	[deg]	$\cdots$
Galactic latitude	GLAT	[deg]	$\cdots$
Containment radius ($68\%$)	pos_68	[deg]	$\cdots$
Containment radius ($95\%$)	pos_95	[deg]	$\cdots$
TS	TS	$\cdots$	$\cdots$
$T{S}_{\mathrm{curv}}^{\mathrm{PLE}}$	TS_curv	$\cdots$	$\cdots$
Integrated photon flux between $E=[0.3,500]$ GeV	Flux300	[ph cm-2 s-1]	Unc_F300
Integrated energy flux between $E=[0.3,500]$ GeV	Energy_Flux300	[MeV cm-2 s-1]	Unc_Energy_Flux300
Functional form of the SED	SpectrumType	$\cdots$	$\cdots$
Spectral index	Spectral_Index	$\cdots$	Unc_Spectral_Index
Cutoff energy (for PLE)	Cutoff	[MeV]	Unc_Cutoff
Curvature parameter, $\beta$ (for LP)	beta	$\cdots$	Unc_beta
IEM with which the source is detected	IEM	$\cdots$	$\cdots$
Associated 3FGL source	3FGL_Name	$\cdots$	$\cdots$
Classification of 3FGL source	3FGL_Class	$\cdots$	$\cdots$
Cluster membership (Off. IEM)	Cluster_Off	$\cdots$	$\cdots$
Cluster membership (Alt. IEM)	Cluster_Alt	$\cdots$	$\cdots$

Full text

1. Introduction

The Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope has been operating since 2008. It has produced the most detailed and precise maps of the $\gamma$-ray sky and collected more than 200 million extraterrestrial $\gamma$ rays in the energy range 0.05–2000 GeV.

The region toward the Galactic Center (GC) is the brightest direction in LAT maps. Along this line of sight (l.o.s.) $\gamma$ rays originate primarily in diffuse processes: interactions of primary cosmic-ray (CR) nuclei with the interstellar gas, bremsstrahlung scattering of CR electrons and positrons with interstellar gas, and inverse Compton scattering of photons from interstellar radiation fields. The LAT also detects individual sources such as pulsars, compact binary systems, supernova remnants, and blazars. In the last seven years many groups analyzing LAT data have reported the detection of an excess of $\gamma$-ray emission at GeV energies with an extent of about $20^{\circ}$ from the GC (we will refer to this as the GC excess).

The GC excess is found with respect to predictions based on a variety of interstellar emission models (IEMs), point source catalogs, and selections of LAT data (e.g., Goodenough & Hooper 2009; Abazajian & Kaplinghat 2012; Hooper & Slatyer 2013; Gordon & Macías 2013; Abazajian et al. 2014; Calore et al. 2015; Daylan et al. 2016). This excess is well modeled with a spherically symmetric generalized Navarro-Frenk-White (NFW) (Navarro et al. 1996; Kravtsov et al. 1998) density profile with index $\alpha=1.25$, $\rho(r)=\rho_{0}/(r/r_{s}(1+r/r_{s}))^{\alpha}$, and its spectral energy distribution (SED) in the inner $10^{\circ}$ from the GC is peaked at a few GeV with an intensity that is approximately one tenth of the total $\gamma$-ray intensity.

The Fermi-LAT Collaboration has performed an analysis using 5.2 years of the Pass 7 reprocessed data in the energy range 1 to 100 GeV for the $15^{\circ}\times 15^{\circ}$ region around the GC. This analysis constructed four dedicated IEMs and produced a point-source catalog (designated 1FIG) which includes 48 sources detected in each of the four IEMs with a Test Statistic (TS) larger than 25 111The TS is defined as twice the difference in maximum log-likelihood between the null hypothesis (i.e., no source present) and the test hypothesis: $TS=2(\log\mathcal{L}_{\rm test}-\log\mathcal{L}_{\rm null})$ (see, e.g., Mattox et al. 1996).  (Ajello et al. 2016).

Recently, the Fermi-LAT Collaboration published an updated analysis (Ackermann et al. 2017) using data from 6.5 years of observation and the new Pass 8 event-level analysis (Atwood et al. 2013). The Pass 8 event-level analysis significantly improves the acceptance, direction and energy reconstruction, and enables sub-selection of events based on the quality of the direction reconstruction. In this updated analysis further investigations of the systematic uncertainties of modeling the diffuse emission region were made using a variety of templates for additional diffuse $\gamma$-ray emission components, such as a data-motivated template for the Fermi bubbles (Su et al. 2010; Ackermann et al. 2014), and with an additional population of electrons used in modeling the central molecular zone, and with three different point source lists.

These two analyses confirm the existence of the GC excess. However, the energy spectrum of the excess is found to depend significantly on the choice of IEM and source list (Ajello et al. 2016; Ackermann et al. 2017).

Different interpretations have been proposed to explain the GC excess. Its approximately spherical morphology and energy spectrum are compatible with $\gamma$ rays emitted from a Galactic halo of dark matter (DM). This possibility has been studied in many papers (e.g., Goodenough & Hooper 2009; Abazajian et al. 2014; Calore et al. 2015; Daylan et al. 2016) and the intensity and shape of the GC excess has been found to be compatible with DM particles with mass 40–60 GeV annihilating through the $b\bar{b}$ channel with a thermally averaged cross section close to the canonical prediction for thermal relic DM (roughly $3\times 10^{-26}{\rm cm}^{3}{\rm s}^{-1}$, e.g., Steigman et al. 2012).

However, if DM exists and gives rise to this excess the same particles should also produce measurable emission from dwarf spheroidal satellite galaxies of the Milky Way, which are known to be DM dominated  (Strigari et al. 2008). No evidence of such a flux from dwarf galaxies has been detected so far, and the limits obtained for the annihilation cross section are in tension with the DM interpretation of the GC excess (see Albert et al. 2017, and references therein).

Among alternative interpretations proposed are that the GeV excess is generated by recent outbursts of CR protons interacting with gas via neutral pion production (Carlson & Profumo 2014) or of CR leptons inverse Compton scattering interstellar radiation (Petrović et al. 2014; Cholis et al. 2015a; Gaggero et al. 2015). However, the hadronic scenario predicts a $\gamma$-ray signal that is significantly extended along the Galactic plane and correlated with the distribution of gas, which is highly incompatible with the observed characteristics of the excess (Petrović et al. 2014). The leptonic outburst scenario plausibly leads to a signal that is more smoothly distributed and spherically symmetric; however, it requires at least two outbursts to explain the morphology and the intensity of the excess with the older outbursts injecting more-energetic electrons. An additional population of supernova remnants near the GC that steadily injects CRs is also a viable interpretation for the GC excess (Gaggero et al. 2015; Carlson et al. 2016).

Recently, evidence of the existence of an unresolved population of $\gamma$-ray sources in the inner $20^{\circ}$ of the Galaxy with a total flux and spatial distribution consistent with the GC excess has been published by Lee et al. (2015) and Bartels et al. (2016). These faint sources have been interpreted as belonging to the Galactic bulge PSR population. This interpretation has been investigated, by, e.g., Cholis et al. (2015b) who claim that about 60 Galactic bulge pulsars should have been already present in Fermi-LAT catalogs, though they may not yet be firmly identified as PSRs.

Several authors have examined the properties of detected $\gamma$-ray pulsars (PSRs) and found that the unresolved pulsars in the Galactic bulge could account for a significant fraction of the GC excess (Mirabal 2013; Grégoire & Knödlseder 2013; Yuan & Zhang 2014; Petrović et al. 2015; Brandt & Kocsis 2015; O’Leary et al. 2015, 2016). Throughout this paper we will use “PSR” to refer specifically to detectable $\gamma$-ray pulsars; i.e., pulsars that emit $\gamma$-rays and whose $\gamma$-ray beams cross the Earth. Typically these have spin-down luminosities above the observed “deathline” of $\sim 3\times 10^{33}{\rm\,erg}{\rm s}^{-1}$ (Guillemot et al. 2016). The SEDs of PSRs are compatible with the GC excess spectrum and $\mathcal{O}(1000)$ are required to explain its intensity (Hooper & Goodenough 2011; Abazajian 2011; Calore et al. 2014; Cholis et al. 2015b). This Galactic bulge pulsar population is hypothesized to be distinct from the well-known “disk” population that follows the Galactic spiral arms and from which we detect the mostly local known sample of pulsars in radio and $\gamma$ rays (Manchester et al. 2005; Abdo et al. 2013). See, in particular, Figure 2 of  Calore et al. (2016) for an illustration of the Galactic disk and bulge pulsar populations. Finally, since this PSR population would be distributed in the Galactic bulge its spatial morphology could be consistent with that of the GC excess.

The large amount of data collected by the LAT after 7.5 years of operation and the improvement in energy and spatial resolution brought by Pass 8 enable a deeper search for PSRs in the Galactic bulge. Such a search is highly relevant to testing the potential PSR nature of the GC excess.

Prospects for detecting radio pulsations from the bulge pulsar population were studied by Calore et al. (2016), and the authors found that existing radio pulsar surveys using the Parkes (Keith et al. 2010) and Green Bank (Stovall et al. 2013) telescopes are not quite sensitive enough to detect many pulsars from the bulge population. On the other hand, large area surveys using, e.g., MeerKAT and later SKA (Kramer & Stappers 2015) should detect dozens to hundreds of pulsars from the Galactic bulge.

In this paper we investigate the pulsar interpretation of the GC excess deriving a new catalog of sources detected in the GC region and selecting among them PSR candidates using SED-based criteria. In Section 2 we present the data selection and the background models that we use. In Section 3 we describe our analysis pipeline and derive a list of sources in the GC region. In Section 4 we discuss potential systematic biases between the source lists and PSR candidate lists derived with two different IEMs. In Section 5 we study the SEDs of blazars and PSRs detected by the LAT and introduce a criterion to select PSR candidates from our list of sources. In Section 6 we study the distribution of luminosities of known $\gamma$-ray PSRs. Details of the data analysis pipeline are provided in the appendices.

2. Data selection and background models

The analysis presented in this paper uses 7.5 years of Fermi-LAT data recorded between 2008 August 4 and 2016 February 4 (Fermi mission elapsed time 239557418–476239414 s). We apply the standard data-quality selections222See http://fermi.gsfc.nasa.gov/ssc/data/analysis..

Since we are interested in detecting the emission from point sources, we select events belonging to the “Pass 8 Source” event class and use the corresponding P8R2_SOURCE_V6 instrument response functions. In order to reduce the contamination of $\gamma$ rays generated by cosmic-ray interactions in the upper atmosphere to negligible levels we select events with a maximum zenith angle of $90^{\circ}$. We use events in the energy range $E=[0.3,500]$ GeV.

We select data for a region of interest (ROI) that is a square of side $40^{\circ}$ centered on the GC $(l,b)=(0^{\circ},0^{\circ})$, where $l$ and $b$ are the Galactic longitude and latitude, since the GC excess has an extension of approximately $20^{\circ}$.

We employ two different IEMs to estimate the systematic uncertainties introduced by the choice of IEM. The IEMs are brightest in the Galactic disk where the density of interstellar gas and radiation fields is greatest. Additionally, isotropic emission, mainly due to $\gamma$-ray emission from unresolved sources (see, e.g., Di Mauro & Donato 2015) and residual contamination from interactions of charged particles in the LAT misclassified as $\gamma$ rays, is included in the model (Ackermann et al. 2015).

The first IEM we use is the gll_iem_v06.fits template, released with Pass 8 data (Acero et al. 2016). The corresponding isotropic component is the iso_P8R2_SOURCE_V6_v06.txt template 333For descriptions of these templates see http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html.. These are routinely used for Pass 8 analyses and we refer to this model as the official (Off.) model. The second IEM is the Sample model (Ackermann et al. 2017, Section 2.2), from which we remove the GC excess component and add the Fermi Bubbles template at $|b|<10^{\circ}$ (Ackermann et al. 2017, Section 5.1.3). We refer to that model as the alternate (Alt.) model.

3. Analysis pipeline and source list

We use the Fermipy Python package (version 00-11-00)555See http://fermipy.readthedocs.io/en/latest/. in conjunction with standard LAT ScienceTools666See https://fermi.gsfc.nasa.gov/ssc/. (version 10-03-00) to find and characterize point sources for both IEMs.

To break the analysis into manageable portions we subdivide the $40^{\circ}\times 40^{\circ}$ ROI into 64 smaller $8^{\circ}\times 8^{\circ}$ ROIs with an overlap of $3^{\circ}$ between adjacent ROIs. Sources near the edge of an ROI are thus well contained in an adjacent ROI. Considering the entire $40^{\circ}\times 40^{\circ}$ ROI would imply several hundred free parameters, making the analysis with the LAT ScienceTools prohibitive. In each ROI we bin the data with a pixel size of $0.06^{\circ}$ and 8 energy bins per decade. In general we analyze each ROI separately; however, as discussed below, at certain points in the analysis we merge information from the analyses of the different ROIs.

The first step of the analysis is to find sources in each of the 64 ROIs. For each ROI we construct an initial model consisting of the IEM, the isotropic template and sources detected with $TS>49$ in the Fermi LAT Third Source Catalog (3FGL, Acero et al. 2015). This provides a reasonably good initial representation of the $\gamma$-ray data in each ROI. The procedure selects 116 3FGL sources that we include in the $40^{\circ}\times 40^{\circ}$ ROI. As we will show later in this section we recover the vast majority of the least significant 3FGL sources (i.e., those with $TS$ values ranging from 25 to 49).

We then use Fermipy tools to refine the positions and the SED parameters of 3FGL sources for the larger, Pass 8 data set that we use here, as well as to find new sources in each ROI. The details of this procedure are described in Appendix A. Since the ROIs overlap slightly, as part of this procedure we remove duplicate sources found in more than one ROI.

We detect 374 (385) sources with $TS>25$ when using the Off. (Alt.) IEM model. Combining the list of detected sources with each IEM we detect 469 unique sources of which 290 are found with both models. The positions of these sources are displayed in Figure 1, overlaid on a counts map for the $40^{\circ}\times 40^{\circ}$ ROI. By comparison, the 3FGL catalog contains 202 sources in this region and 189 (182) of them are found with our analysis with the Off. (Alt.) IEM. The 1FIG, which covered only the inner $15^{\circ}\times 15^{\circ}$, contains 48 sources of which we find 38 (41) when we employ in the analysis the Off. (Alt.) IEM. We define associations with 3FGL and 1FIG sources based on the relative positions and the $95\%$ localization uncertainty regions reported in those catalogs and found in our analysis. Specifically, we require that the angular distances of sources in the 3FGL or 1FIG from matching sources in our analysis be smaller than the sum in quadrature of the $95\%$ containment angles in 3FGL or 1FIG and in our analysis. The 3FGL and 1FIG sources that are not present in our lists either have TS near the detection threshold (i.e., $25<TS<36$) or are located within $0\fdg 5$ of the GC.

The GC region is the brightest in the $\gamma$-ray sky and developing a model of the interstellar emission in this region is very challenging (see, e.g., Calore et al. 2015; Ajello et al. 2016; Ackermann et al. 2017). Imperfections of our IEMs could manifest themselves as dense concentrations of sources in regions where the IEMs particularly under-predict the diffuse intensity. To account for this, we employed a source cluster-finding algorithm (described in Appendix B) to identify such regions. We find a total of four clusters of sources with four or more sources within $0\fdg 6$ of at least one other source in the cluster. These clusters are located around the GC, in regions around the W28 and W30 supernova remnants and near 3FGL J1814.1$-$1734c, which is an unassociated source in the 3FGL catalog. (The ‘c’ designation means that it was flagged in that catalog as possibly an artifact.) These clusters are shown in Figure 1.

We removed from further consideration here all sources identified as belonging to clusters.

4. Systematic biases in source-finding and

PSR selection criteria efficiency from uncertainties in the IEMs

In this section we discuss potential systematic biases between the source lists and PSR candidate lists constructed using the “Off.” and “Alt.” IEMs. We find that away from the Galactic plane ($|b|>2.5^{\circ}$) these systematic uncertainties are similar in magnitude to the statistical uncertainties, and that along the plane the systematic differences between the two models are randomly enough distributed that they do not introduce significant biases to the analysis. Finally, we argue that the effect of the uncertainties of the IEMs on the analysis of the Galactic bulge PSR population can be quantified reasonably by evaluating the difference in the results obtained using the two different IEMs.

Following the method we used to associate sources with the 3FGL and 1FIG lists, we match sources between the lists produced with the Off. and Alt. IEMs if the separation between two positions is less than the sum in quadrature of the 95% containment radii of the source localization. In the left panel of Figure 2 we show the significance (i.e., $\sqrt{TS}$) as measured with the two IEMs, considering all source candidates with $TS>16$ in either IEM. We note the high degree of correlation between the significance found with the two different IEMs; the absolute scatter (i.e., the RMS of the difference) of the sources found using both IEMs is $1.7\sigma$ ($4.6\sigma$) for sources with $|b|>2.5^{\circ}$ ($|b|<2.5^{\circ}$). Furthermore, sources were detected with slightly larger significance on average with the Alt. IEM; the difference was $0.2\sigma$ ($1.8\sigma$) for sources away from (along) the plane.

Many of the sources are faint and are detected with TS near the threshold of 25. On the other hand, the requirement $TS_{\rm curv}>9$ in the PSR selection criteria effectively sets a somewhat higher flux threshold, as it is difficult to measure spectral curvature for a faint source near the detection threshold. Therefore it is equally important to consider the effects of fluctuations near threshold on $TS_{\rm curv}$. The $\sqrt{TS_{\rm curv}}$ values found using the two IEMs are shown in the right-hand panel of Figure 2. The absolute scatter of $\sqrt{TS_{\rm curv}}$ is 1.0 (3.0) for sources with $|b|>2.5^{\circ}$ ($|b|<2.5^{\circ}$).

We also explicitly study the differences between the two IEMs and how those correlate with differences in the sources of the respective lists. In Figure 3 we compare locations of sources found with only one of the two IEMs with the estimated statistical significance of the difference between the two. Specifically we used the gtmodel tool to produce the expected counts from diffuse IEMs (including the isotropic component) and computed the difference divided by the square-root of the mean of the two IEMs. The white boxes and circles show sources that are detected with $TS>49$ using one IEM, but that are not detected with the other, and thus primarily attributable to differences in the two IEMs. The large majority of such cases occur near the Galactic plane, which is unsurprising as it is also the region where the difference between the two IEMs is the most significant. On the other hand, the differences between the two IEMs are less statistically significant (generally between $\pm 1\sigma$) away from the plane, and most of the sources found with only one of the two IEMs are either near threshold (i.e, $25\leq TS<49$, green markers) or are seen as sub-threshold source candidates ($16\leq TS<25$, cyan markers) with the other IEM.

From these studies we conclude that away from the Galactic plane ($|b|>2.5^{\circ}$) the differences between the two lists are largely attributable to the combination of thresholding and variations between the IEMs that are sub-dominant to the statistical variations of the data. Furthermore, $TS=25$ for four degrees of freedom (source position, spectral index and normalization) corresponds to a detection at the level of 4.1 standard deviations, and the systematic uncertainties in the source significance are smaller than the statistical uncertainties. Therefore, even considering the systematic uncertainties attributable to the IEMs, the pre-trials significance of these source candidates is close to 3 standard deviations.

To estimate how much of the difference between the two lists can be attributed to the combination of small systematic variations with thresholding effects, we simulated the effect of the significance variations attributable to the differences between the two IEMs by adjusting the $\sqrt{TS}$ of each source by a random number drawn from a Gaussian distribution with mean and width given by $\mu=0$, $\sigma=1.7$ ($\sigma=2.6)$ for sources away from (along) the Galactic plane and testing if the adjusted $TS$ value was above the $TS=25$ detection threshold. In this simulation we found that 79% (58%) of the sources found with either IEM would be expected to be found with both IEMs, for sources away from (along) the Galactic plane. In the actual lists, these numbers were 64% and 56%. Since the scatter away from the plane is small compared to the $5.0\sigma$ detection threshold, we believe that ratio between the 79% overlap found in our toy simulation and the 64% overlap found in the two cases sets a lower bound that $64/76=0.81$ of the sources found away from the plane are with either IEM are in fact real. Along the plane, on the other hand, the scatter is comparable in magnitude to the detection threshold. There we believe that the 56% overlap of sources found with both IEMs the lower bound on the fraction of real sources.

We then studied individual ROIs near the Galactic plane to better understand the interplay between the implementation details of our data analysis pipeline and the difference in the IEMs. Specifically, since we fit the normalization of the Galactic diffuse emission in each of the ROIs, and iteratively add point sources in the largest positive residuals, these steps give the resulting models some freedom to compensate for inaccuracies in the IEMs.

Figure 4 is similar to Figure 3, but was produced for a single one of our $8^{\circ}\times 8^{\circ}$ ROIs, and using the fitted values of the IEM normalization.

Figure 5 shows the significance of the residuals for the same ROI as computed with the two different IEMs. Overall, some correlation is evident between the residuals in the two maps, but there are no regions larger than a few degrees across where both models over-predict the data.

To quantify the finding that the two IEMs do not both over-predict the data in the same regions, we have evaluated histograms of the estimated significance of the residuals, $({\rm data}-{\rm model})/\sqrt{\rm model}$, for $0.48^{\circ}\times 0.48^{\circ}$ spatial bins (we simply rebinned the model and data maps combining $8\times 8$ pixels) for both IEMs for all 16 ROIs along the Galactic plane (Fig. 6).

These distributions are fairly similar for the two IEMs. The distribution for the Off. IEM is narrower and more symmetric and is well described by a Gaussian with $\mu=0.1$ and $\sigma=1.6$. On the other hand, for the Alt. IEM, the distribution is wider, and has a non-Gaussian tail extending to around $5\sigma$. This indicates that there are in fact regions where the Alt. IEM over-predicts the data, but that they are moderately-sized. Since we obtain similar results with the Off. and Alt IEMs, we believe that the dominant effects of the uncertainties of the IEMs is to widen the effective threshold. Since there are more sources just below threshold that just above threshold, we expect that this will cause us to detect more sources that we would see without the effect of the IEMs.

From this we conclude that any biases in the source lists due to the IEMs are likely to be slight overestimates of the number of sources. These overestimates would be somewhat larger along the Galactic plane than away from it.

We have also studied the possibility that biases in the IEMs could affect the spectral parameters of the sources, causing either true pulsars to be rejected, or increasing the number of non-pulsars mistakenly selected as PSR candidates. We found that the correlation between the spectral parameters and the background model is modest: the correlation coefficient between the normalization of the IEM and the $\Gamma$ and $E_{\rm cut}$ spectral parameters for the new sources ranges from $-0.25$ to $0.10$. Given the very large statistics along the Galactic plane, small fractional biases in the IEM could cause marked biases in the spectral parameters. Therefore, we have looked at the agreement between the energy fluxes in the low-energy bins and the spectral models and found no evidence of significant biases caused by the IEMs. In Figure 7 we present a histogram of the residual of the single- energy bin fluxes with respect to the spectral models, for two different energy bins. For the lowest energy bin the scaled residuals are very nearly normally distributed ($\mu=0.05$, $\sigma=1.0$); for the energy bin near 2 GeV model slightly overestimates the data ($\mu=-0.2$, $\sigma=1.0$). Overall we believe that these results indicate that the spectra are not significantly biased by any errors in the IEMs.

From these studies we conclude that the systematic variations near the plane appear to be distributed in such a manner that no substantial regions have both IEMs either under- or over-estimate the data. We also conclude that the effect on the measured spectral parameters is small. Therefore, we believe that the systematic uncertainties in the properties of the PSR population are well described by the differences obtained with the two IEMs, or would tend to bias the fitting procedure to find lower numbers of PSRs in the Galactic bulge.

5. SED of pulsars and blazars

In Fermi-LAT catalogs, blazars are the most numerous source population. Blazars are classified as BL Lacertae (BL Lacs) or Flat Spectrum Radio Quasars (FSRQs) depending on the presence of strong emission optical lines. In the 3FGL $95\%$ of BL Lac and $85\%$ of FSRQ spectra are modeled with a power-law (PL) function while blazars with a significant spectral curvature (only about 10% of the entire population) are modeled with a log-parabola (LP)777See http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html for a description of the spectral models implemented in the LAT ScienceTools. (see, e.g., Ajello et al. 2012; Di Mauro et al. 2014; Ackermann et al. 2015; Ghisellini et al. 2017, for a characterization of the Fermi blazar population). On the other hand, a power law with exponential cutoff (PLE) at a few GeV is the preferred model for pulsars (Abdo et al. 2013). Of the 167 PSRs reported in the 3FGL (143 PSRs identified by pulsations and 24 sources spatially associated with radio pulsars) 115 have spectral fits parametrized with a PLE because they have a significant spectral curvature (see, e.g., Lorimer 2004; Calore et al. 2014, for a characterization of the $\gamma$-ray and radio pulsar population). The functional definitions of the PL, LP, and PLE spectra are given in Acero et al. (2015).

As described above, spectral shape is a promising observable to separate PSRs from blazars. We fit the spectrum of each source in the ROI and derive the likelihood values for both the PL ($\mathcal{L}_{\rm{PL}}$) and PLE ($\mathcal{L}_{\rm{PLE}}$) spectra. We introduce for each source in our analysis the TS for a curved spectrum as: $TS^{\rm{PLE}}_{\rm{curv}}=2\cdot(\log{\mathcal{L}_{\rm{PLE}}}-\log{\mathcal{L}_{\rm{PL}}})$. This parameter quantifies the preference to model an SED with a PLE with respect to a PL.

We perform the same analysis on known PSRs and blazars to study the distribution of spectra of these two populations and develop criteria to select PSR candidates. We use the public list of $\gamma$-ray PSRs with 210 sources888See https://confluence.slac.stanford.edu/display/GLAMCOG/Public+List+of+LAT-Detected+Gamma-Ray+Pulsars. and the sub-sample of sources identified with or associated with blazars in the 3FGL catalog that have significant spectral curvature. Our blazar sub-sample includes all 3FGL blazars that have Signif_Curve greater than 3, where Signif_Curve is the significance in standard deviations of the likelihood improvement between PL and LP spectra. We use this sub-sample to study those blazars most likely to be incorrectly flagged as PSR candidates. This reduced sample of blazars contains 218 objects.

Our definition of $TS^{\rm{PLE}}_{\rm{curv}}$ is slightly different from the 3FGL Signif_Curve parameter ($\sigma_{\rm{curv}}$) in that $TS^{\rm{PLE}}_{\rm{curv}}$ is defined as the likelihood improvement for a PLE spectrum with respect to the PL spectrum. Furthermore, the 3FGL catalog analysis was based on only 4 years of LAT data. Therefore we used Fermipy to re-analyze $10^{\circ}\times 10^{\circ}$ ROIs centered around each source in this sample of 210 $\gamma$-ray PSRs and 218 blazars. From this re-analysis we derived $TS^{\rm{PLE}}_{\rm{curv}}$, the photon index ($\Gamma$) and the energy cutoff ($E_{\rm{cut}}$) for the PLE spectrum.

Of the 210 PSRs, 172 (169) were found to have $TS^{\rm{PLE}}_{\rm{curv}}>9$. The average and standard deviation of their photon indices and energy cutoffs were $\Gamma=1.33\pm 0.54(1.30\pm 0.54)$ and$\log_{10}(E_{\rm{cut}}[\rm{MeV}])=3.43\pm 0.24(3.40\pm 0.24)$ when we employed the Off. (Alt.) IEM.

In Table 1 we report the average photon index and cutoff energy for young PSRs (rotational period $P$ greater than 30 ms) and millisecond PSRs (MSPs). The energy cutoff parameter is consistent between young PSRs and MSPs while the average photon index of MSPs is slightly harder.

In the sample of 218 blazars with Signif_Curve $>3$, 153 have $TS^{\rm{PLE}}_{\rm{curv}}>9$. In the left-hand panel of Figure 8 we show $\Gamma$ and $\log_{10}(E_{\rm{cut}})$ for PSRs and blazars detected with $TS^{\rm{PLE}}_{\rm{curv}}>9$. The two populations are well separated in the plotted SED parameters. Taking $\Gamma<2.0$ and $E_{\rm{cut}}<10$ GeV as selection criteria (shown in cyan in the figure) only 12 blazars, $7\%$ of our blazar sample and less than $1\%$ of the entire 3FGL blazar population, are incorrectly flagged as PSR candidates. Clearly, these selection criteria are effective for distinguishing the PSRs from blazars. Additional studies of the efficiency and false-positive rate of these selection criteria using simulated data are described in Appendix E.

We apply the PSR candidate selection criteria to our source lists. In the list derived with the Off. (Alt.) IEM we find 86 (115) PSR candidates. If we require that the source is selected with both IEMs we find 66 PSR candidates. In the right-hand panel of Figure 8 we show the $\Gamma$ and $\log_{10}(E_{\rm{cut}})$ for all sources detected with $TS^{\rm{PLE}}_{\rm{curv}}>9$ for the analysis with the Off. IEM. The average SED parameters for PLE are shown in Table 1. For PSR-like sources detected with both IEMs the photon index ($1.02\pm 0.52$) is harder with respect to known PSRs (see fifth row in Table 1). This is due to observational biases for the detection of PSRs in direction of the inner part of our Galaxy. We will show this in Section 7.2. We also calculate the integrated energy flux ($S=\int^{E_{\rm{max}}}_{E_{\rm{min}}}EdN/dEdE$) over the range from $E_{\rm{min}}=300$ MeV to $E_{\rm{max}}=500$ GeV.

The full list of sources detected with $TS>25$ is provided as a FITS file and described in detail in Appendix C. We designate the sources with the prefix ‘2FIG’ designation; however we emphasize that many of the fainter sources in the list are detected only with one of the two IEMs we used.

Globular clusters are gravitationally bound concentrations of ten thousand to one million stars and are the most ancient constituents of our Milky Way Galaxy. They are known to contain many pulsars. Among the detected sources we have 11 globular clusters already identified in the 3FGL and among those, 6 satisfy the PSR-like criteria.

6. Luminosity Distribution of PSRs

Of the fraction of the 210 identified $\gamma$-ray PSRs for which we have good distance estimates (roughly 100 have reported fractional uncertaintes smaller that 25%), the large majority are located within 4 kpc of the Solar System (see, e.g., Figure 3 of Abdo et al. 2013). They are thus “local” and belong to the Galactic disk population. A pulsar interpretation of the GC excess requires a Galactic bulge population of PSRs. (Throughout this paper we adopt 8.5 kpc as the distance to the GC.) However the known Galactic disk PSR population is a strong foreground to the putative Galactic bulge population. In this section and the next we describe simulations of the Galactic disk and bulge populations that are based on the morphology and energy spectrum of the GC excess and the characteristics of 3FGL PSRs. We then use these simulations to estimate the number of sources in these populations that would be needed to match both the observed GC excess and the numbers and properties of the detected sources in the $40^{\circ}\times 40^{\circ}$ ROI.

To perform these simulations, we use Fermipy as explained in Appendix D to generate simulated data sets of the individual ROIs and then use the analysis pipeline described in Section 3 to analyze those simulated maps. For these simulations we used the same time and energy ranges and the same ROI as the analysis on the real sky, and we employed the Off. IEM.

For the simulations described in this section we simulate blazars isotropically distributed in the $40^{\circ}\times 40^{\circ}$ ROI, with fluxes taken from the $dN/dS$ derived by Abdo et al. (2010b) using the 1FGL catalog (Abdo et al. 2010a). We simulated individual blazars using a PL SED with $\Gamma$ extracted from a Gaussian distribution with average 2.40 and standard deviation 0.30 as found for the 3FGL (Acero et al. 2015). We simulated blazars with an energy flux integrated between $0.3-500$ GeV of $>9\times 10^{-8}$ MeV cm*-2* s*-1* in order to have a sizable number of simulated sources below the detection threshold.

To model the Galactic disk and bulge PSR populations, ideally we would start with a known luminosity function for PSRs, or derive one starting with the 210 publicly announced $\gamma$-ray PSRs. However, because of complications including the incompleteness of the radio pulsar sample and variations in detection efficiency across the sky, and since the PSR sample covers the entire sky well beyond our ROI, the PSR luminosity function is poorly constrained and difficult to extract. See, e.g., Strong (2007); Cholis et al. (2014); Petrović et al. (2015); Bartels et al. (2016) both for previous estimates of the luminosity function and for discussions of the complications.

Therefore we adopt a staged approach. We first assume a PL shape $dN/dL\propto L^{-\beta}$ for the luminosity function and estimate the slope ($\beta$) from the data as described below. We then derive the normalization given that slope by using simulations to estimate the number of $\gamma$-ray pulsars that would be required to explain the GC excess and energy flux distribution of 3FGL sources with curved SEDs using simulations. Finally, we reuse those simulations to derive the efficiency for PSRs to pass our selection criteria.

Given for each PSR the energy flux $S$ and distance $d$ we calculate the luminosity: $L=4\pi Sd^{2}$, where $S$ is integrated in the energy range $0.3-500$ GeV, 999Note that this differs from previous publications, which use 100 MeV as the lower bound of the integration range for the luminosity. as derived with the analysis described in Section 5 and the distance ($d$) is taken from the ATNF catalog version 1.54 (Manchester et al. 2005) using the continuously updated web page 101010We always use the Dist_1 parameter, namely the best distance estimate available, when it exists (see the ATNF catalog for more information: http://www.atnf.csiro.au/people/pulsar/psrcat/).. The catalog provides distance measures for 135 out of 210 PSRs and from these we can derive the observed luminosity distribution $dN/dL$. The missing 75 pulsars are mostly Fermi $\gamma$-loud and radio-quiet. We also calculate the $dN/dL$ separately for young PSRs and MSPs without correcting for the detection efficiency.

Since the PSR sample detected by the LAT is known to be incomplete and we do not correct for the detection efficiency, we select sources within a distance of 1.5 kpc from the Earth. Indeed, considering luminosities in the range $[3\times 10^{33},10^{36}]$ erg s*-1* and a distance of 1.5 kpc implies energy fluxes in the range $[7\times 10^{-6},3\times 10^{-3}]$ MeV cm*-2* s*-1*, for which the LAT efficiently detects sources. In Figure 9 we show the luminosity distribution for our sample of PSRs with $d<1.5$ kpc. We then perform a fit to the data starting from $L=3\times 10^{33}$ erg s*-1* to avoid the change in slope at low luminosities due to the incompleteness of the LAT detections at the low-luminosity end. We use a PL shape $dN/dL\propto L^{-\beta}$ and the fit yields $\beta=1.7\pm 0.3$. Our fit differs from the data points only below $10^{33}$ erg s*-1* where it is difficult to identify PSRs with $\gamma$-ray data.

For reference, the faintest $\gamma$-ray pulsar is PSR J0437-4715 with $L=3.55\times 10^{31}$ erg s*-1*, and the faintest $\gamma$-ray pulsar found in a blind search is PSR J1741-2054 with $L=1.51\times 10^{33}$ erg s${{}^{-}1}$. Although the selection we applied in source luminosity and distance reduce significantly the incompleteness of pulsar sample, the efficiency for the detection of PSRs for this distance and luminosity range might be smaller than 100% and thus the intrinsic luminosity distribution even steeper; however we will use the value of $\beta=1.7$ to examine the issue of how many PSRs would be required to produce the GCE.

Our estimate of $\beta$ for PSRs with $L_{\gamma}>3\times 10^{33}$ erg s*-1* is similar to that found for MSPs by Cholis et al. (2014); Hooper & Mohlabeng (2016); Winter et al. (2016). In these papers a break at around $10^{33}$ erg s*-1* or a slightly curved luminosity function is considered. However, since the slope of the luminosity function is 1.7, the integrated luminosity is dominated by the bright sources. Therefore, a change of $dN/dL$ below $10^{33}$ erg s*-1* does not significantly affect our results. We also point out that Winter et al. (2016) have estimated the completeness of the Second Fermi-LAT catalog of pulsars (Abdo et al. 2013) finding that it is almost 100% for pulsars with luminosity greater than $10^{35}$ erg s*-1*. The least-luminous PSR detected is $3.55\cdot 10^{31}$ erg s*-1* while the most luminous is $10^{36}$ erg s*-1*. We therefore simulate luminosities between $10^{31}$ erg s*-1* and $10^{36}$ erg s*-1* to include PSRs below the current LAT detection threshold. Furthermore, throughout this paper, we quote the total number of PSRs with $L=[10^{33},10^{36}]$ erg s*-1* in the Galactic disk ($N_{\rm disk}$) and bulge ($N_{\rm bulge}$) to specify the normalization of $dN/dL$.

7. Simulating the Galactic PSR population

In this section we report our assumptions for the disk and Galactic bulge populations of PSRs and explain how we simulate these two populations.

7.1. Galactic Disk PSRs

For our simulations we use the Galactocentric spatial distribution $\rho(R)$ as modeled by Lorimer (2004): $\rho(R)\propto R^{n}\exp{(-R/\sigma)}$ with $n=2.35$ and $\sigma=1.528$ kpc. The dependence on the distance from the Galactic disk is modeled with an exponential cutoff $\rho(z)\propto\exp{(-|z|/z_{0})}$ with scale height $z_{0}=0.70$ kpc as in Calore et al. (2014). The luminosity function is modeled as a PL with index 1.70 over the range $L=[10^{31},10^{36}]$ erg s*-1*, see Section 6.

Analyses of the Galactic disk pulsar population estimate that it could contain thousands of objects (Levin et al. 2013; Lorimer 2013, 2004). These estimates are derived from radio catalogs of pulsars, correcting their spatial distribution for observational biases and using information for the star formation rate and distribution in the Galaxy.

However, the radio and the $\gamma$-ray emission are only slightly correlated and many of the nearest radio pulsars are not detected by the LAT. The current ATNF catalog lists 714 pulsars within 3 kpc of the Earth that have measured spin down energy loss rates ($\dot{E}$). Of these, 257 have $\dot{E}>10^{33}$ erg s*-1*, the observed minimum for which pulsars emit $\gamma$ rays (Guillemot et al. 2016). The LAT has detected about 30% of these, most likely primarily due to differences in radio and $\gamma$-ray emission beam solid angles of the source and to their distances.

In short, the overall number of $\gamma$-ray PSRs in the disk population is not very well constrained. A lower limit is given by the identified $\gamma$-ray PSRs: in particular, for fluxes integrated between $0.3-500$ GeV larger than $10^{-5}$ MeV cm*-2* s*-1* where the efficiency for the detection of PSRs is almost $100\%$. In Figure 10 we show the energy flux histogram for 3FGL PSRs.

This is, however, only a lower limit because many non-radio PSRs may be present as sources in the 3FGL, but the pulsations have not yet been detected in $\gamma$ rays. Without timing solutions from radio observations, the detection of $\gamma$-ray pulsation is challenging; see e.g., Dormody et al. (2011) for a sensitivity estimate. To obtain an estimate of the upper limit of PSRs in the disk we have selected the 3FGL unassociated sources with curvature significance greater than 3. We added their flux distribution to that of the detected PSRs (Figure 10). We expect that the bright tail ($S>1.8\cdot 10^{-5}$ MeV cm*-2* s*-1*) of the flux distribution for the disk population of PSRs should fall between the already detected PSRs (111) and the sum of this with 3FGL unassociated sources with $\sigma_{\rm{curv}}>3$ (237). This range is represented by the cyan band in Figure 10. From this we estimate that the Galactic disk PSR population consists of between 4000 and 16000 $\gamma$-ray emitting sources with $L=[10^{33},10^{36}]$ erg s*-1* with $\beta=1.7$, i.e., $N_{\rm disk}=[4000,16000]$. This result is derived with the Galactocentric spatial distribution as modeled by Lorimer (2004) and it slightly depends on the assumed radial distribution of PSRs in the disk. The flux distribution of this disk population is displayed with a gray band in Figure 10 and for energy fluxes larger than $10^{-5}$ MeV cm*-2* s*-1* is perfectly consistent with the cyan band.

7.2. Galactic Bulge PSRs

In this section, we consider the properties that would be needed for a Galactic bulge population of pulsars to generate the GC excess, and find good agreement with previous works (Abazajian & Kaplinghat 2012; Brandt & Kocsis 2015; Bartels et al. 2016; Hooper & Mohlabeng 2016).

We model the spatial distribution of the Galactic bulge PSR population as spherically symmetric with respect to the GC with a radial profile $dN/dr\propto r^{-\alpha}$ for $r<3$ kpc and 0 elsewhere and with $\alpha=2.60$ in order to approximately match a generalized NFW with slope of 1.3. This spatial distribution is consistent with the morphology of the GC excess (Calore et al. 2015; Daylan et al. 2016; Ajello et al. 2016; Ackermann et al. 2017) and gives a latitude profile of the $\gamma$-ray intensity from PSRs with a similar shape to the excess (see left-hand panel Figure 11). As with the disk population, we model the luminosity function as a PL with $\beta=1.7$ over the range $L=[10^{31},10^{36}]$ erg s*-1*. For each simulated PSR we draw a location and luminosity from the relevant spatial distribution and $\gamma$-ray PSR luminosity function and sample values from the distributions of $\Gamma$ and $\log_{10}(E_{\rm{cut}})$ given in the fifth row of Table 1. We then derive the SED of each PSR, and simulate PSRs until their total energy spectrum is of the same intensity as the GC excess as reported by Calore et al. (2015) and Ackermann et al. (2017).

In the left panel of Figure 11 we compare the average latitude profile from 20 simulations with the intensity of the GC excess and in the right panel we compare the total SED from simulated PSRs and the GC excess spectrum as derived by Calore et al. (2015) and Ackermann et al. (2017). The gray band in both plots is derived from the possible range of source counts that this component could contain and includes the reported systematic uncertainties of the latitude profile and energy spectrum of the GC excess. This procedure finds that a Galactic bulge PSR population capable of generating the GC excess would need to include 800–3600 sources with $L=[10^{33},10^{36]}$ erg s*-1*, i.e., $N_{\rm bulge}=[800,3600]$.

8. Conclusions

We have analyzed 7.5 years of Pass 8 Fermi-LAT data for the energy range $[0.3,500]$ GeV in a $40^{\circ}\times 40^{\circ}$ around the GC in order to provide a list of PSR candidates and test the pulsar interpretation of the GC excess. Employing two IEMs we detect about 400 sources, a factor of about two more than in the 3FGL catalog (Acero et al. 2015) and five more than in the 1FIG (Ajello et al. 2016) (derived for $[1,100]$ GeV and for $15^{\circ}\times 15^{\circ}$); these latter analyses were based on shorter time intervals of data than we consider here. We then studied the SEDs of $\gamma$-ray PSRs and 3FGL blazars using a PLE shape and found that the distributions of photon index and energy cutoff parameters for these two populations are very well separated, with typical values of $\Gamma<2$ and $E_{\rm{cut}}<10$ GeV for PSRs. Moreover, about 82% of PSRs and only 9% of blazars have $TS^{\rm{PLE}}_{\rm{curv}}>9$. We thus use the selection criteria $TS^{\rm{PLE}}_{\rm{curv}}>9$, $\Gamma<2$ and $E_{\rm{cut}}<10$ GeV to extract PSR candidates from our seed list, finding 66 sources detected with both IEMs.

We took the distribution of spectral parameters from the 210 identified $\gamma$-ray PSRs and the luminosity distribution of PSRs within 1.5 kpc 88footnotemark: 8. We used parameters given by  Lorimer (2004); Calore et al. (2014) to model the spatial distribution of the disk population of PSRs. With this model, we find that given the number and distributions of unassociated 3FGL sources with curved SED we constrain the number of Galactic disk PSRs to be in the range $[4000,16000]$ for a luminosity function with slope 1.7 and $L_{\gamma}=[10^{33},10^{36}]$ erg s*-1*. Similarly, we used the latitude profile and energy spectrum of the GC excess (e.g., Calore et al. 2015; Ackermann et al. 2017) to model the Galactic bulge PSR population and found that it must include 800–3600 sources (most of them unresolved) if it is to explain the GeV excess.

Facilities: Fermi.

The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States; the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France; the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy; the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK), and Japan Aerospace Exploration Agency (JAXA) in Japan; and the K. A. Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Etudes Spatiales in France. MDM and EC acknowledge support by the NASA Fermi Guest Investigator Program 2014 through the Fermi multi-year Large Program N. 81303 (P.I. E. Charles). We would like to thank Richard Bartels, Dan Hooper, Tim Linden, Siddhartha Mishra-Sharma, Nicholas Rodd, Benjamin Safdi and Tracy Slatyer for helping us to identify an error in the maximum likelihood analysis of the Galactic bulge and disk PSR populations that was included in a previous version of this paper, (arXiv:1705.00009, v1). Their work is described in a note will be posted on arXiv at the same time as this draft (Bartels et al. 2017).

Appendix A A. Analysis pipeline and description of Fermipy tools

We analyze each ROI with a pipeline based on the Fermipy package and the Fermi Science Tools. In the following description, we denote with italics the Fermipy methods and configuration parameters used in each step of the pipeline.

We start the analysis of each region with a model that includes 3FGL sources with TS $>49$, the IEM, and the isotropic template and begin by optimizing the spatial and spectral parameters of this model. We first perform a global fit of the spectral parameters for all components in the model. For the global fit we retain the spectral model (PL, LP or PLE) reported by the 3FGL. We then relocalize all 3FGL point sources using the localize method. This method generates a map of the model likelihood versus source position in the vicinity of the nominal 3FGL position and finds the best-fit position and errors by fitting a 2D parabola to the log-likelihood values in the vicinity of the peak. When localizing a source, we free the normalization of the IEM and isotropic template and spectral parameters of sources within 3∘ of the source of interest. After relocalizing 3FGL sources, we repeat the global fit of the spectral parameters of all components.

On average, 3FGL sources move by $0.04^{\circ}$ in the relocalization step. This is of the same order as the $68\%$ location uncertainty radius for most 3FGL sources (Acero et al. 2015). As an example, in the left panel of Figure 12 we show the result of the relocalization for 3FGL J1709.5$-$0335. The new position is offset by $0.087^{\circ}$ with respect to the 3FGL position and the 3FGL $68\%$ positional uncertainty is $0.064^{\circ}$.

After relocalizing the 3FGL sources, we add new source candidates to the model using the find_sources method. This method iteratively refines the model by identifying peaks in a TS map of the region with $\sqrt{TS}>$ sqrt_ts_threshold and adding a new source at the position of each peak. After each iteration a new TS map is generated with an updated background model that incorporates sources found in the previous iteration. This procedure is repeated until no peaks are found with amplitude larger than sqrt_ts_threshold. To minimize the likelihood of finding multiple peaks associated with a single source, the algorithm restricts the separation between peaks found in an iteration to be greater than min_separation by excluding peaks that are within this distance of a peak with higher TS.

We run find_sources with a point-source test source model with a PL spectrum and a fixed photon index of $\Gamma=2$. We use sqrt_ts_threshold $=4$ and min_separation$=0.4^{\circ}$. The result is a list of source candidates with $TS>16$. On average we detected about six sources with $TS>25$ per ROI.

In the right-hand panel of Figure 12 we display an example TS map derived prior to the source-finding step when only the IEM, isotropic template and 3FGL sources with $TS>49$ are included. From this figure we note that a region with $TS>25$ is located in the center of the ROI, near $(l,b)=(17.5^{\circ},17^{\circ})$ from which our analysis extracted three new sources with $TS>25$.

We derive the SED for each source candidate in our list using the sed method. A likelihood analysis is performed in each energy bin independently, using for the spectrum of the source a PL shape with a fixed photon index of 2 and normalization free to vary. In this procedure we also leave free the normalizations of the IEM, isotropic template, and the normalizations of sources within $3^{\circ}$ of the source of interest. As an example, in Figure 13 the SED of 3FGL J1730.6$-$0357 is reported together with a fit of a PLE spectrum. 3FGL J1730.6$-$0357 has $TS^{\rm{PLE}}_{\rm{curv}}=40$ meaning that it has a significant curvature. The fit with a PLE in fact gives a spectral index of 0.5 and energy cutoff of 1.4 GeV.

To avoid finding duplicate sources in regions where our ROIs overlap, we remove sources that are found in more than one ROI and that have an angular separation smaller than $0.2^{\circ}$. Specifically we keep the version of the source that is closest to the center of the ROI in which it was found.

Appendix B B. Minimum spanning tree source cluster-finding

To avoid placing spurious point sources in regions where the IEM under-predicted the Galactic diffuse emission we applied a source clustering algorithm based on the minimum spanning tree (MST) algorithm. The MST algorithm calculates how to connect points with minimum total length of the connections (Kruskal 1956). The fermipy.cluster_sources module applies the MST algorithm to identify clusters of sources by joining sets of sources within a maximum connection distance (the $dist$ parameter) and retaining those clusters with at least a specified minimum number of sources (the $nsrc$ parameter). We tested different values for $dist$ and $nsrc$. We found that using $nsrc\geq 4$ generally avoided creating spurious clusters around chance spatial coincidences. For $nsrc=4$, we found that using $dist=0.6^{\circ}$ selected a few clusters that were spatially associated with known regions of complicated Galactic diffuse emission, such as supernova remnants and the GC. Using larger values, such as $dist=1.0^{\circ}$, resulted in clusters consisting of chains of sources up to $\sim 4^{\circ}$ long in the Galactic plane. Using $dist\leq 0.3^{\circ}$ and $nsrc=4$, resulted in no clusters being found with either IEM. The same four clusters found using the Off. IEM and the values $dist=0.6^{\circ}$ and $nsrc=4$ were also found using $dist=0.7^{\circ}$, and only the cluster near the GC was found using $dist=0.5^{\circ}$. In light of these studies, we adopted the value $dist=0.6^{\circ}$ and $nsrc=4$ and obtain the results presented in Section 3.

Appendix C C. Source list and contents of FITS file

Together with this paper, we are releasing the list of sources detected in our analysis as a FITS file. The file contains a single binary table with the source data. The list includes 469 sources detected with $TS>25$ in a region with $|b|<20^{\circ}$ and $|l|<20^{\circ}$. The table has one row per source; the column names and contents are described in Table 2. When applicable the units of the columns are given by the header keywords following the FITS standard. All of the spectral parameters are taken from the ROI optimization procedure described in Appendix A.

Appendix D D. Generating simulated data with Fermipy

We use the simulate_roi method to simulate the binned $\gamma$-ray counts data in each ROI using the maximum-likelihood model of the ROI. Specifically, the method generates “model cubes” of the expected number of $\gamma$-ray counts in each pixel and energy bin in the ROI for the time interval of our analysis. The method then generates Poisson-distributed random numbers with expectation values drawn from the model cube for each pixel and energy bin and produces a simulated binned counts maps for each ROI. This procedure results in simulated $\gamma$-ray counts maps that are statistically identical to those produced with gtobssim, which simulates individual $\gamma$ rays and convolves them with the instrument response model. The simulate_roi method is many times faster than gtobssim, making the extensive simulations we have performed much more tractable.

Appendix E E. Testing PSR selection criteria with simulated data

As discussed in Section 5, $\sim$90% of blazars in the 3FGL catalog have an SED modeled with PL shape while the remaining 10% are modeled with a LP. On the other hand, about 82% of PSRs have energy spectra consistent with a PLE. Employing the parameter $TS^{\rm{PLE}}_{\rm{curv}}$ and making spectral fits of blazars and PSRs with a PLE model we have shown that the criteria $TS^{\rm{PLE}}_{\rm{curv}}>9$, $\Gamma<2.0$ and $E_{\rm{cut}}<10$ GeV work very well to separate the PSR and blazar populations.

In this Appendix we investigate how these criteria work for a simulated test population of sources with curved spectra compatible with a PLE, but with a slightly larger energy cutoff and softer photon index with respect to PSRs. We want to test if this additional population would severely contaminate our PSR candidates. For this we simulate a bulge population of PSRs as explained in Section 7.2 and 1500 sources with a photon index of $2.3\pm 0.2$ and energy cutoff of $\log_{10}{(E_{\rm{cut}}[\rm{MeV}])}=4.48\pm 0.25$ uniformly distributed in the GC region. We choose these distributions for the energy spectrum parameters to demonstrate that a putative population of sources with a curved SED and with a distribution of $\Gamma$ and $E_{\rm{cut}}$ that is fairly well separated from the PSR-like criteria is not going to contaminate significantly our selection of PSR-like sources because of mis-estimation of the spectral parameters.

We simulate fluxes of sources for this population from the source count distribution of blazars derived by Abdo et al. (2010b). We use the same analysis as used for the derivation of the source list in the real sky. For the SED we consider a PLE shape and evaluate the best-fit parameters for $\Gamma$ and $\log_{10}(E_{\rm{cut}})$. Then we select the sources detected with $TS>25$ and $TS^{\rm{PLE}}_{\rm{curv}}>9$. The result for the values of $\Gamma$ and $\log_{10}(E_{\rm{cut}})$ is displayed in Figure 14. Also the detected sources satisfying $TS^{\rm{PLE}}_{\rm{curv}}>9$ maintain a very good separation in the $\Gamma-\log_{10}(E_{\rm{cut}})$ plane. Only $6\%$ of the non-PSR sources detected with $TS>25$ have measured $\Gamma<2.0$ and $E_{\rm{cut}}<10$ GeV and $TS^{\rm{PLE}}_{\rm{curv}}>9$. This result means that the presence of a putative source population with an SED modeled with a PLE but with a softer photon index and higher-energy cutoff would produce a contamination that is small with respect to our PSR candidates.