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

M. Ajello, L. Baldini, J. Ballet, G. Barbiellini, D. Bastieri, R. Bellazzini, E. Bissaldi, R.D. Blandford, E.D. Bloom, E. Bottacini, et al.

► To cite this version:

M. Ajello, L. Baldini, J. Ballet, G. Barbiellini, D. Bastieri, et al. Characterizing the population of pulsars in the inner Galaxy with the Fermi Large Area Telescope. 2020. hal-02024227

HAL Id: hal-02024227

https://hal.science/hal-02024227

Preprint submitted on 5 Aug 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
CHARACTERIZING THE POPULATION OF PULSARS IN THE INNER GALAXY WITH THE *FERMI* LARGE AREA TELESCOPE.

M. AJELLO<fnref n="1" />, L. BALDINI<fnref n="2" />, J. BALLET<fnref n="3" />, G. BARBIELLINI<fnref n="4" />,<fnref n="5" />, D. BASTIERI<fnref n="6" />,<fnref n="7" />, R. BELLAZZINI<fnref n="8" />, E. BISSALDI<fnref n="9" />,<fnref n="10" />, R. D. BLANDFORD<fnref n="11" />, E. D. BLOOM<fnref n="11" />, E. BOTTACINI<fnref n="11" />, J. BREGEON<fnref n="12" />, P. BRUEL<fnref n="13" />, R. BUEHLER<fnref n="14" />, R. A. CAMERON<fnref n="11" />, R. CAPUTO<fnref n="15" />, M. CARAGIULO<fnref n="9" />,<fnref n="10" />, P. A. CARAVEO<fnref n="16" />, E. CAVAZZUTI<fnref n="17" />, C. CECCI<fnref n="18" />,<fnref n="19" />, E. CHARLES<fnref n="11" />,<fnref n="20" />,*, A. CHEKHTMAN<fnref n="21" />, G. CHIARO<fnref n="7" />, S. CIPRINI<fnref n="22" />,<fnref n="18" />, D. COSTANTIN<fnref n="7" />, F. COSTANZA<fnref n="10" />, F. D'AMMANDO<fnref n="23" />,<fnref n="24" />, F. DE PALMA<fnref n="10" />,<fnref n="25" />, R. DESIANTE<fnref n="26" />,<fnref n="27" />, S. W. DIGEL<fnref n="11" />, N. DI LALLA<fnref n="2" />, M. DI MAURO<fnref n="11" />,<fnref n="28" />,*, L. DI VENERE<fnref n="9" />,<fnref n="10" />, C. FAVUZZI<fnref n="9" />,<fnref n="10" />, E. C. FERRARA<fnref n="29" />, A. FRANCKOWIAK<fnref n="14" />, Y. FUKAZAWA<fnref n="30" />, S. FUNK<fnref n="31" />, P. FUSCO<fnref n="9" />,<fnref n="10" />, F. GARGANO<fnref n="10" />, D. GASPARRINI<fnref n="22" />,<fnref n="18" />, N. GIGLIETTO<fnref n="9" />,<fnref n="10" />, F. GIORDANO<fnref n="9" />,<fnref n="10" />, M. GIROLETTI<fnref n="23" />, D. GREEN<fnref n="32" />,<fnref n="29" />, L. GUILLEMOT<fnref n="33" />,<fnref n="34" />, S. GUIRIEC<fnref n="29" />,<fnref n="35" />, A. K. HARDING<fnref n="29" />, D. HORAN<fnref n="13" />, G. JÓHANNESSON<fnref n="36" />,<fnref n="37" />, M. KUSS<fnref n="8" />, G. LA MURA<fnref n="7" />, S. LARSSON³⁸,³⁹, L. LATRONICO<fnref n="26" />, J. LI⁴⁰, F. LONGO<fnref n="4" />,<fnref n="5" />, F. LOPARCO<fnref n="9" />,<fnref n="10" />, M. N. LOVELLETTE⁴¹, P. LUBRANO<fnref n="18" />, S. MALDERA<fnref n="26" />, D. MALYSHEV<fnref n="31" />, L. MARCOTULLI<fnref n="1" />, P. MARTIN⁴², M. N. MAZZIOTTA<fnref n="10" />, M. MEYER⁴³,³⁹, P. F. MICHELSON<fnref n="11" />, N. MIRABAL<fnref n="29" />,<fnref n="35" />, T. MIZUNO⁴⁴, M. E. MONZANI<fnref n="11" />, A. MORSELLI⁴⁵, I. V. MOSKALENKO<fnref n="11" />, E. NUSS<fnref n="12" />, N. OMODEI<fnref n="11" />, M. ORIENTI<fnref n="23" />, E. ORLANDO<fnref n="11" />, M. PALATIELLO<fnref n="4" />,<fnref n="5" />, V. S. PALIYA<fnref n="1" />, D. PANEQUE⁴⁶, J. S. PERKINS<fnref n="29" />, M. PERSIC<fnref n="4" />,<fnref n="7" />, M. PESCE-ROLLINS<fnref n="8" />, F. PIRON<fnref n="12" />, G. PRINCIPE<fnref n="31" />, S. RAINÒ<fnref n="9" />,<fnref n="10" />, R. RANDO<fnref n="6" />,<fnref n="7" />, M. RAZZANO<fnref n="8" />,⁴⁸, A. REIMER⁴⁹,<fnref n="11" />, O. REIMER⁴⁹,<fnref n="11" />, P. M. SAZ PARKINSON⁵⁰,⁵¹,⁵², C. SGRÒ<fnref n="8" />, E. J. SISKIND⁵³, D. A. SMITH⁵⁴, F. SPADA<fnref n="8" />, G. SPANDRE<fnref n="8" />, P. SPINELLI<fnref n="9" />,<fnref n="10" />, H. TAJIMA⁵⁵,<fnref n="11" />, J. B. THAYER<fnref n="11" />, D. J. THOMPSON<fnref n="29" />, L. TIBALDO⁵⁶, D. F. TORRES⁴⁰,⁵⁷, E. TROJA<fnref n="29" />,<fnref n="32" />, G. VIANELLO<fnref n="11" />, K. WOOD⁵⁸, M. WOOD<fnref n="11" />,⁵⁹,*, G. ZAHARIJAS⁶⁰,⁶¹

Draft version October 31, 2017

ABSTRACT

An excess of γ-ray emission from the Galactic Center (GC) region with respect to predictions based on a variety of interstellar emission models and γ-ray source catalogs has been found by many groups using data from the *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 γ-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° × 40° 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.

<footnote n="1">¹ Department of Physics and Astronomy, Clemson University, Kinard Lab of Physics, Clemson, SC 29634-0978, USA</footnote>

<footnote n="2">² Università di Pisa and Istituto Nazionale di Fisica Nucleare, Sezione di Pisa I-56127 Pisa, Italy</footnote>

<footnote n="3">³ Laboratoire AIM, CEA-IRFU/CNRS/Université Paris Diderot, Service d'Astrophysique, CEA Saclay, F-91191 Gif sur Yvette, France</footnote>

<footnote n="4">⁴ Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, I-34127 Trieste, Italy</footnote>

<footnote n="5">⁵ Dipartimento di Fisica, Università di Trieste, I-34127 Trieste, Italy</footnote>

<footnote n="6">⁶ Istituto Nazionale di Fisica Nucleare, Sezione di Padova, I-35131 Padova, Italy</footnote>

<footnote n="7">⁷ Dipartimento di Fisica e Astronomia “G. Galilei”, Università di Padova, I-35131 Padova, Italy</footnote>

<footnote n="8">⁸ Istituto Nazionale di Fisica Nucleare, Sezione di Pisa, I-56127 Pisa, Italy</footnote>

<footnote n="9">⁹ Dipartimento di Fisica “M. Merlin” dell’Università e del Politecnico di Bari, I-70126 Bari, Italy</footnote>

<footnote n="10">¹⁰ Istituto Nazionale di Fisica Nucleare, Sezione di Bari, I-70126 Bari, Italy</footnote>

<footnote n="11">¹¹ W. W. Hansen Experimental Physics Laboratory, Kavli Institute for Particle Astrophysics and Cosmology, Department of Physics and SLAC National Accelerator Laboratory, Stanford University, Stanford, CA 94305, USA</footnote>

<footnote n="12">¹² Laboratoire Univers et Particules de Montpellier, Université Montpellier, CNRS/IN2P3, F-34095 Montpellier, France</footnote>

<footnote n="13">¹³ Laboratoire Leprince-Ringuet, École polytechnique, CNRS/IN2P3, F-91128 Palaiseau, France</footnote>

<footnote n="14">¹⁴ Deutsches Elektronen Synchrotron DESY, D-15738 Zeuthen, Germany</footnote>

<footnote n="15">¹⁵ Center for Research and Exploration in Space Science and Technology (CRESST) and NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="16">¹⁶ INAF-Istituto di Astrofisica Spaziale e Fisica Cosmica Milano, via E. Bassini 15, I-20133 Milano, Italy</footnote>

<footnote n="17">¹⁷ Italian Space Agency, Via del Politecnico snc, 00133 Roma, Italy</footnote>

<footnote n="18">¹⁸ Istituto Nazionale di Fisica Nucleare, Sezione di Perugia, I-06123 Perugia, Italy</footnote>

<footnote n="19">¹⁹ Dipartimento di Fisica, Università degli Studi di Perugia, I-06123 Perugia, Italy</footnote>

<footnote n="20">²⁰ email: echarles@slac.stanford.edu</footnote>

<footnote n="21">²¹ College of Science, George Mason University, Fairfax, VA 22030, resident at Naval Research Laboratory, Washington, DC 20375, USA</footnote>

<footnote n="22">²² ASI Space Science Data Center, Via del Politecnico snc, 00133 Roma, Italy</footnote>

<footnote n="23">²³ INAF Istituto di Radioastronomia, I-40129 Bologna, Italy</footnote>

<footnote n="24">²⁴ Dipartimento di Astronomia, Università di Bologna, I-40127 Bologna, Italy</footnote>

<footnote n="25">²⁵ Università Telematica Pegaso, Piazza Trieste e Trento, 48, I-80132 Napoli, Italy</footnote>

<footnote n="26">²⁶ Istituto Nazionale di Fisica Nucleare, Sezione di Torino, I-10125 Torino, Italy</footnote>

<footnote n="27">²⁷ Università di Udine, I-33100 Udine, Italy</footnote>

<footnote n="28">²⁸ email: mdimauro@slac.stanford.edu</footnote>

<footnote n="29">²⁹ NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="30">³⁰ Department of Physical Sciences, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8526, Japan</footnote>

<footnote n="31">³¹ Erlangen Centre for Astroparticle Physics, D-91058 Erlangen, Germany</footnote>

<footnote n="32">³² Department of Physics and Department of Astronomy, University of Maryland, College Park, MD 20742, USA</footnote>

<footnote n="33">³³ Laboratoire de Physique et Chimie de l'Environnement et de l'Espace – Université d'Orléans / CNRS, F-45071 Orléans Cedex 02, France</footnote>

<footnote n="34">³⁴ Station de radioastronomie de Nançay, Observatoire de Paris, CNRS/INSU, F-18330 Nançay, France</footnote>

<footnote n="35">³⁵ NASA Postdoctoral Program Fellow, USA</footnote>

<footnote n="36">³⁶ Science Institute, University of Iceland, IS-107 Reykjavik, Iceland</footnote>

<footnote n="37">³⁷ Nordita, Roslagstullsbacken 23, 106 91 Stockholm, Sweden</footnote>
# 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 γ-ray sky and collected more than 200 million extraterrestrial γ 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.) γ 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 γ-ray emission at GeV energies with an extent of about 20° 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° from the GC is peaked at a few GeV with an intensity that is approximately one tenth of the total γ-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^1$ (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 γ-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 γ-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} \text{cm}^3\text{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

<footnote n="1">¹ The 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 L_{\text{test}} - \log L_{\text{null}})$ (see, e.g., Mattox et al. 1996).</footnote>

<footnote n="38">³⁸ Department of Physics, KTH Royal Institute of Technology, AlbaNova, SE-106 91 Stockholm, Sweden</footnote>

<footnote n="39">³⁹ The Oskar Klein Centre for Cosmoparticle Physics, AlbaNova, SE-106 91 Stockholm, Sweden</footnote>

<footnote n="40">⁴⁰ Institute of Space Sciences (IEEC-CSIC), Campus UAB, Carrer de Magrans s/n, E-08193 Barcelona, Spain</footnote>

<footnote n="41">⁴¹ Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA</footnote>

<footnote n="42">⁴² CNRS, IRAP, F-31028 Toulouse cedex 4, France</footnote>

<footnote n="43">⁴³ Department of Physics, Stockholm University, AlbaNova, SE-106 91 Stockholm, Sweden</footnote>

<footnote n="44">⁴⁴ Hiroshima Astrophysical Science Center, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8526, Japan</footnote>

<footnote n="45">⁴⁵ Istituto Nazionale di Fisica Nucleare, Sezione di Roma “Tor Vergata”, I-00133 Roma, Italy</footnote>

<footnote n="46">⁴⁶ Max-Planck-Institut für Physik, D-80805 München, Germany</footnote>

<footnote n="47">⁴⁷ Osservatorio Astronomico di Trieste, Istituto Nazionale di Astrofisica, I-34143 Trieste, Italy</footnote>

<footnote n="48">⁴⁸ Funded by contract FIRB-2012-RBFR12PM1F from the Italian Ministry of Education, University and Research (MIUR)</footnote>

<footnote n="49">⁴⁹ Institut für Astro- und Teilchenphysik and Institut für Theoretische Physik, Leopold-Franzens-Universität Innsbruck, A-6020 Innsbruck, Austria</footnote>

<footnote n="50">⁵⁰ Santa Cruz Institute for Particle Physics, Department of Physics and Department of Astronomy and Astrophysics, University of California at Santa Cruz, Santa Cruz, CA 95064, USA</footnote>

<footnote n="51">⁵¹ Department of Physics, The University of Hong Kong, Pokfulam Road, Hong Kong, China</footnote>

<footnote n="52">⁵² Laboratory for Space Research, The University of Hong Kong, Hong Kong, China</footnote>

<footnote n="53">⁵³ NYCB Real-Time Computing Inc., Lattingtown, NY 11560-1025, USA</footnote>

<footnote n="54">⁵⁴ Centre d'Études Nucléaires de Bordeaux Gradignan, IN2P3/CNRS, Université Bordeaux 1, BP120, F-33175 Gradignan Cedex, France</footnote>

<footnote n="55">⁵⁵ Solar-Terrestrial Environment Laboratory, Nagoya University, Nagoya 464-8601, Japan</footnote>

<footnote n="56">⁵⁶ Max-Planck-Institut für Kernphysik, D-69029 Heidelberg, Germany</footnote>

<footnote n="57">⁵⁷ Institució Catalana de Recerca i Estudis Avançats (ICREA), E-08010 Barcelona, Spain</footnote>

<footnote n="58">⁵⁸ Praxis Inc., Alexandria, VA 22303, resident at Naval Research Laboratory, Washington, DC 20375, USA</footnote>

<footnote n="59">⁵⁹ email: mdwood@slac.stanford.edu</footnote>

<footnote n="60">⁶⁰ Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, and Università di Trieste, I-34127 Trieste, Italy</footnote>

<footnote n="61">⁶¹ Laboratory for Astroparticle Physics, University of Nova Gorica, Vipavska 13, SI-5000 Nova Gorica, Slovenia</footnote>

<footnote>\* Corresponding authors: M. Di Mauro, mdimauro@slac.stanford.edu; E. Charles, echarles@slac.stanford.edu; M. Wood, md-wood@slac.stanford.edu.</footnote>
interstellar radiation (Petrović et al. 2014; Cholis et al. 2015a; Gaggero et al. 2015). However, the hadronic scenario predicts a γ-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 γ-ray sources in the inner 20° 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 γ-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 γ-ray pulsars; i.e., pulsars that emit γ-rays and whose γ-ray beams cross the Earth. Typically these have spin-down luminosities above the observed “deathline” of ~3 × 10^33 ergs^−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 γ 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 γ-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 selections<fnref n="2" />.

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 γ rays generated by cosmic-ray interactions in the upper atmosphere to negligible levels we select events with a maximum zenith angle of 90°. 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° centered on the GC ($l, b$) = (0°, 0°), where $l$ and $b$ are the Galactic longitude and latitude, since the GC excess has an extension of approximately 20°.

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 γ-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 γ 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 <fnref n="3" />. 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°$ (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)<fnref n="5" /> in conjunction with standard LAT *ScienceTools*<fnref n="6" /> (version 10-03-00) to find and characterize point sources for both IEMs.

To break the analysis into manageable portions we subdivide the 40° × 40° ROI into 64 smaller 8° × 8° ROIs with an overlap of 3° between adjacent ROIs. Sources near the edge of an ROI are thus well contained in an adjacent ROI. Considering the entire 40° × 40° ROI would imply several

<footnote n="2">² See http://fermi.gsfc.nasa.gov/ssc/data/analysis.</footnote>

<footnote n="3">³ For descriptions of these templates see http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html.</footnote>

<footnote n="5">⁵ See http://fermipy.readthedocs.io/en/latest/.</footnote>

<footnote n="6">⁶ See https://fermi.gsfc.nasa.gov/ssc/.</footnote>
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° 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 γ-ray data in each ROI. The procedure selects 116 3FGL sources that we include in the 40° × 40° 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 Fermi 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° × 40° 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° × 15°, 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.5° of the GC.

The GC region is the brightest in the γ-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 clustering 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.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.

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Counts map of the 40° × 40° ROI used in this analysis. The map includes data for the range [0.3, 500] GeV. The map is in Hammer-Aitoff projection, centered on the GC and in Galactic coordinates. The pixel size is 0.1°. The color scale shows the number of photons per pixel. Markers are shown at the positions of sources found in our analysis with the Off. IEM. White markers show sources associated with a 3FGL source and green markers show new sources with no 3FGL counterpart. Stars (squares) indicate sources that are (not) PSR-like and purple markers indicate sources belonging to a cluster, and the clusters are outlined with purple circles (see text for details). Finally, blue stars show PSRs identified as or associated in the 3FGL.</figcaption></figure>

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°$) 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°$ ($|b| < 2.5°$). 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_{curv} > 9$ in the PSR selection criteria effectively sets a somewhat higher flux threshold, as it is difficult to measure
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Left: correlation between the significance ($\sqrt{TS}$) derived using the Off. and Alt. IEMs. The black (red) points are for sources with $|b| > 2.5^\circ$ ($|b| < 2.5^\circ$). Sources that are detected with $TS > 16$ for only one of the two IEMs are assigned a $TS$ value near zero for the other IEM for plotting purposes. Right: same, but for $\sqrt{TS}_{curv}$; again, sources detected with both only one IEM have been a $TS$ value near zero for the other IEM.</figcaption></figure>

spectral curvature for a faint source near the detection thresh-
old. Therefore it is equally important to consider the effects of
fluctuations near threshold on $TS_{curv}$. The $\sqrt{TS}_{curv}$ values
found using the two IEMs are shown in the right-hand panel
of Figure 2. The absolute scatter of $\sqrt{TS}_{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 pro-
duce 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 us-
ing 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 dif-
ference between the two IEMs is the most significant. On
the other hand, the differences between the two IEMs are less
statistically significant (generally between ±1σ) away from
the plane, and most of the sources found with only one of
the two IEMs are either near threshold (i.e., 25 ≤ TS < 49,
green markers) or are seen as sub-threshold source candidates
(16 ≤ TS < 25, cyan markers) with the other IEM.

From these studies we conclude that away from the Galac-
tic plane (|b| > 2.5°) 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 sta-
tistical 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 uncertain-
ties. 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

<figure><img src="image_3.png" /><figcaption>**Figure 3.** The color scale shows the statistical significance of the difference between the two IEMs integrated over the range $E = [0.3, 500]$ GeV. Sources found with both IEMs are indicated with 'x' markers. Square (circle) markers show sources found only with the Off. (Alt.) IEM. The cyan markers show sources that have a corresponding sub-threshold candidate with $16 \le TS < 25$ using the other IEM. Green markers show low significance sources with $25 \le TS < 49$; white circles and squares show high significance sources with $TS \ge 49$.</figcaption></figure>

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° × 8° ROIs, and using the fitted values of the IEM normalization.

<figure><img src="image_4.png" /><figcaption>**Figure 4.** Same as Figure 3, but for one of the 8° × 8° ROIs used in the analysis pipeline, and for the complete model map at the end of the analysis pipeline, using the fitted values of the IEM normalization.</figcaption></figure>

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, $(\text{data} - \text{model})/\sqrt{\text{model}}$, for 0.48° × 0.48° spatial bins (we simply rebinned the model and data maps combining 8 × 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_{\text{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)<fnref n="7" /> (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 ($L_{\text{PL}}$) and PLE ($L_{\text{PLE}}$) spectra. We introduce for each source in our analysis the TS for a curved spectrum as:

<footnote n="7">⁷ See http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html for a description of the spectral models implemented in the LAT ScienceTools.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Estimated significance of residuals, $(\text{data} - \text{model})/\sqrt{\text{model}}$, for the Off. (left) and Alt. (right). The markers have the same meanings as for Figures 3 and 4. Gaussian smoothing with a width of $\sigma = 0.18^\circ$ has been applied.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 6.** Distribution of estimated significance of residuals for $0.48^\circ \times 0.48^\circ$ spatial bins in the 16 ROIs along the Galactic plane. The solid lines show Gaussian fits to the two distributions. The best-fit mean and width are $\mu = 0.1, \sigma = 1.6$ ($\mu = 0.3, \sigma = 1.8$) for the Off. (Alt.) IEM.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Distribution of the significance of the residual for the energy flux in a single bin with respect to fitted spectral model for the PSR candidates found with the Off. IEM analysis.</figcaption></figure>

This is shown for the lowest-energy bin (black), and for an energy bin near the peak of the statistical power of the analysis (~2 GeV, blue). The solid curves are best-fit Gaussian distributions to the two histograms; the best-fit parameters are given in the text.

$TS_{\text{curv}}^{\text{PLE}} = 2 \cdot (\log \mathcal{L}_{\text{PLE}} - \log \mathcal{L}_{\text{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 γ-ray PSRs with 210 sources<fnref n="8" /> 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_{\text{curv}}^{\text{PLE}}$ is slightly different from the 3FGL *Signif_Curve* parameter ($\sigma_{\text{curv}}$) in that $TS_{\text{curv}}^{\text{PLE}}$ 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 γ-ray PSRs and 218 blazars. From this re-analysis we derived $TS_{\text{curv}}^{\text{PLE}}$, the photon index ($\Gamma$) and the energy cutoff ($E_{\text{cut}}$) for the PLE spectrum.

Of the 210 PSRs, 172 (169) were found to have $TS_{\text{curv}}^{\text{PLE}} > 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_{\text{cut}}[\text{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 $\text{Signif\_Curve} > 3$, 153 have $TS_{\text{curv}}^{\text{PLE}} > 9$. In the left-hand panel of Figure 8 we show $\Gamma$ and $\log_{10}(E_{\text{cut}})$ for PSRs and blazars detected

<footnote n="8">⁸ See https://confluence.slac.stanford.edu/display/GLAMCOG/Public/List+of+LAT-Detected+Gamma-Ray+Pulsars.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Left: photon index Γ and energy cutoff $E_{\text{cut}}$ [MeV] of PSRs and blazars detected in our analysis with $TS_{\text{curv}}^{\text{PLE}} > 9$. MSPs are shown as blue plus signs and young PSRs as red crosses. Here we are showing blazars in the 3FGL (Acero et al. 2015) catalog with curvature significance as in the 3FGL (Signif_Curve) larger than 3 (green circles). Right: same as in the left panel but applied to sources in our 40° × 40° ROI detected with $TS_{\text{curv}}^{\text{PLE}} > 9$ in our analysis with the Off. IEM (black circles) and Alt. IEM (red crosses).</figcaption></figure>

with $TS_{\text{curv}}^{\text{PLE}} > 9$. The two populations are well separated in the plotted SED parameters. Taking $\Gamma < 2.0$ and $E_{\text{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_{\text{cut}})$ for all sources detected with $TS_{\text{curv}}^{\text{PLE}} > 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_{\text{min}}}^{E_{\text{max}}} EdN/dEdE$) over the range from $E_{\text{min}} = 300$ MeV to $E_{\text{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 uncertainties smaller than 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.

<table><caption>**Table 1**

Spectral parameters of PSR candidates compared with known PSRs.</caption><thead><tr><th>IEM</th><th>N<sub>PSR</sub></th><th>Γ</th><th>log<sub>10</sub>(E<sub>cut</sub>[MeV])</th></tr></thead><tbody><tr><td>Off.</td><td>86</td><td>1.03 ± 0.52</td><td>3.28 ± 0.33</td></tr><tr><td>Alt.</td><td>115</td><td>1.05 ± 0.50</td><td>3.27 ± 0.31</td></tr><tr><td>Alt. ∩ Off. (Off.)</td><td>66</td><td>1.02 ± 0.52</td><td>3.27 ± 0.32</td></tr><tr><td>Alt. ∩ Off. (Alt.)</td><td>66</td><td>1.01 ± 0.51</td><td>3.26 ± 0.30</td></tr><tr><td>Known PSRs (Off.)</td><td>172</td><td>1.33 ± 0.54</td><td>3.43 ± 0.24</td></tr><tr><td>Young PSRs (Off.)</td><td>86</td><td>1.46 ± 0.53</td><td>3.44 ± 0.26</td></tr><tr><td>MSPs(Off.)</td><td>86</td><td>1.20 ± 0.50</td><td>3.42 ± 0.23</td></tr></tbody></table>

**Note.** — Mean values and standard deviations of $\Gamma$ and $\log_{10}(E_{\text{cut}}[\text{MeV}])$ for PSR candidates compared with known $\gamma$-ray PSRs. The first two rows are found using different IEMs (Off. first and Alt. second row). The third (fourth) row is for the PSR candidates detected with both IEMs, computed with the parameters derived in the analysis with the Off. (Alt.) IEM. The last three rows list the parameters for all $\gamma$-ray PSRs, young PSRs and MSPs detected with $TS_{\text{curv}}^{\text{PLE}} > 9$.

(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° × 40° 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° × 40° 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} \text{ MeV cm}^{-2} \text{ 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 γ-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 γ-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,<fnref n="9" /> 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 <fnref n="10" />. 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* γ-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⁻¹ and a distance of 1.5 kpc implies energy fluxes in the range $[7 \times 10^{-6}, 3 \times 10^{-3}]$ MeV cm⁻² s⁻¹, 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⁻¹ 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⁻¹ where it is difficult to identify PSRs with γ-ray data.

For reference, the faintest γ-ray pulsar is PSR J0437-4715 with $L = 3.55 \times 10^{31}$ erg s⁻¹, and the faintest γ-ray pulsar found in a blind search is PSR J1741-2054 with $L = 1.51 \times 10^{33}$ erg s⁻¹. 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.

<figure><img src="image_9.png" /><figcaption>**Figure 9.** Observed luminosities for young PSRs (red data), MSPs (blue data) and the whole population of PSRs with $d < 1.5$ kpc (black data). The best fit to the luminosity distribution for $L > 3 \times 10^{33}$ erg s⁻¹ is also reported (black line). The luminosity is integrated over the energy range [0.3, 500] GeV.</figcaption></figure>

Our estimate of $\beta$ for PSRs with $L_{\gamma} > 3 \times 10^{33}$ erg s⁻¹ 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⁻¹ 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⁻¹ 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⁻¹. The least-luminous PSR detected is $3.55 \times 10^{31}$ erg s⁻¹ while the most luminous is $10^{36}$ erg s⁻¹. We therefore simulate luminosities between $10^{31}$ erg s⁻¹ and $10^{36}$ erg s⁻¹ 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⁻¹ in the Galactic disk ($N_{\text{disk}}$) and bulge ($N_{\text{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⁻¹, 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 γ-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

<footnote n="9">⁹ Note that this differs from previous publications, which use 100 MeV as the lower bound of the integration range for the luminosity.</footnote>

<footnote n="10">¹⁰ We 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/).</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 10.** Flux histogram of 3FGL PSRs alone (red triangles) or added to the flux distribution of unassociated 3FGL sources with curvature $Signif\_Curve > 3$ (black points). The cyan band represents the region between the lower limit (already detected PSRs) and upper limit (3FGL PSRs plus unassociated 3FGL sources with detected spectral curvature). Finally the black curve (gray band) represents the benchmark (band between the minimum and maximum) number of disk PSRs. The flux is integrated over the energy range [0.3, 500] GeV.</figcaption></figure>

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_{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_{disk} = [4000, 16000]$. This result is derived with the Galacticcentric 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_{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_{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_{cut} < 10$ GeV for PSRs. Moreover, about 82% of PSRs and only 9% of blazars have $TS_{curv}^{PLE} > 9$. We thus use the selection criteria $TS_{curv}^{PLE} > 9$, $\Gamma < 2$ and $E_{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$^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
<figure><img src="image_11.png" /><figcaption>**Figure 11.** Left panel: latitude profile of the intensity of the GC excess at 2 GeV and longitude 0° for a simulation of a Galactic bulge PSR population that could reproduce the GC excess, i.e., $N_{\text{bulge}} = [800, 3600]$ for $L = [10^{33}, 10^{36}]$ erg s⁻¹ (gray band and black dashed line), and for the GC excess as found by Calore et al. (2015) (cyan band) and Ackermann et al. (2017) (blue and red solid curves). Right panel: energy spectrum from a simulation of a Galactic bulge PSR population that could reproduce the GC excess, i.e., $N_{\text{bulge}} = [800, 3600]$ for $L = [10^{33}, 10^{36}]$ erg s⁻¹ (gray band and black dashed line), and for the GC excess as found by Calore et al. (2015) (cyan band and blue data) and Ackermann et al. (2017) (orange band and red data).</figcaption></figure>

function with slope 1.7 and $L_{\gamma} = [10^{33}, 10^{36}]$ erg s⁻¹. 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).

REFERENCES

Abazajian, K. N. 2011, JCAP, 3, 010

Abazajian, K. N., Canac, N., Horiuchi, S., & Kaplinghat, M. 2014, Phy. Rev. D, 90, 023526

Abazajian, K. N. & Kaplinghat, M. 2012, Phy. Rev. D, 86, 083511

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010a, ApJS, 188, 405

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010b, ApJ, 720, 435

Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17

Acero, F., Ackermann, M., Ajello, M., et al. 2015, ApJS, 218, 23

Acero, F., Ackermann, M., Ajello, M., et al. 2016, ApJS, 223, 26

Ackermann, M., Ajello, M., Albert, A., et al. 2015, ApJ, 799, 86

Ackermann, M., Albert, A., Atwood, W. B., et al. 2014, ApJ, 793, 64

Ackermann, M. et al. 2015, Astrophys. J., 810, 14

Ackermann, M. et al. 2017, Astrophys. J., 840, 43

Ajello, M., Albert, A., Atwood, W. B., et al. 2016, ApJ, 819, 44

Ajello, M., Shaw, M. S., Romani, R. W., et al. 2012, ApJ, 751, 108

Albert, A., Anderson, B., Bechtol, K., et al. 2017, ApJ, 834, 110

Atwood, W., Albert, A., Baldini, L., et al. 2013, ArXiv:1303.3514

Bartels, R., Krishnamurthy, S., & Weniger, C. 2016, Phys. Rev. Lett., 116, 051102

Bartels, R. et al. 2017, ArXiv:1710.XXXXX

Brandt, T. D. & Kocsis, B. 2015, ApJ, 812, 15

Calore, F., Cholis, I., & Weniger, C. 2015, JCAP, 1503, 038

Calore, F., Di Mauro, M., Donato, F., & Donato, F. 2014, ApJ, 796, 1

Calore, F., Di Mauro, M., Donato, F., Hessels, J. W. T., & Weniger, C. 2016, ApJ, 827, 143

Carlson, E., Linden, T., & Profumo, S. 2016, Phys. Rev. Lett., 117, 111101

Carlson, E. & Profumo, S. 2014, Phy. Rev. D, 90, 023015

Cholis, I., Evoli, C., Calore, F., et al. 2015a, JCAP, 1512, 005

Cholis, I., Hooper, D., & Linden, T. 2014, ArXiv:1407.5583

Cholis, I., Hooper, D., & Linden, T. 2015b, JCAP, 1506, 043

Daylan, T., Finkbeiner, D. P., Hooper, D., et al. 2016, Physics of the Dark Universe, 12, 1

Di Mauro, M. & Donato, F. 2015, Phy. Rev. D, 91, 123001

Di Mauro, M., Donato, F., Lamanna, G., Sanchez, D. A., & Serpico, P. D. 2014, Astrophys. J., 786, 129

Dormody, M., Johnson, R. P., Atwood, W. B., et al. 2011, ApJ, 742, 126

Gaggero, D., Taoso, M., Urbano, A., Valli, M., & Ullio, P. 2015, JCAP, 1512, 056

Ghisellini, G., Righi, C., Costamante, L., & Tavecchio, F. 2017, Mon. Not. Roy. Astron. Soc., 469, 255

Goodenough, L. & Hooper, D. 2009, ArXiv:0910.2998

Gordon, C. & Macías, O. 2013, Phy. Rev. D, 88, 083521

Grégoire, T. & Knödlseder, J. 2013, Astron. Astrophys., 554, A62

Guillemot, L., Smith, D. A., Laffon, H., et al. 2016, A&A, 587, A109

Hooper, D. & Goodenough, L. 2011, Phys. Lett., B697, 412

Hooper, D. & Mohlabeng, G. 2016, ApJ, 1603, 049

Hooper, D. & Slatyer, T. R. 2013, Physics of the Dark Universe, 2, 118

Keith, M. J., Jameson, A., van Straten, W., et al. 2010, MNRAS, 409, 619

Kramer, M. & Stappers, B. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 36
<figure><img src="image_12.png" /><figcaption>**Figure 12.** Left panel: Positional uncertainty contours and best-fit position (solid black lines and x marker) from the localization of 3FGL J1709.5–0335. The color scale shows the difference in TS with respect to the best-fit position of the source. The cross is the position of this source in 3FGL. Right panel: TS map of the ROI evaluated with a test source with a point-source morphology and a PL spectrum with photon index of 2. The background model includes the IEM, isotropic template, and 3FGL sources with $TS > 49$.</figcaption></figure>

Kravtsov, A. V., Klypin, A. A., Bullock, J. S., & Primack, J. R. 1998, ApJ, 502, 48

Kruskal, J. B. 1956, Proceedings of the American Mathematical Society, 7, 48

Lee, S. K., Lisanti, M., & Safdi, B. R. 2015, JCAP, 1505, 056

Levin, L., Bailes, M., Barsdell, B. R., et al. 2013, MNRAS, 434, 1387

Lorimer, D. R. 2004, in IAU Symposium, Vol. 218, Young Neutron Stars and Their Environments, ed. F. Camilo & B. M. Gaensler, 105

Lorimer, D. R. 2013, in IAU Symposium, Vol. 291, Neutron Stars and Pulsars: Challenges and Opportunities after 80 years, ed. J. van Leeuwen, 237–242

Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993

Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396

Mirabal, N. 2013, MNRAS, 436, 2461

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

O'Leary, R. M., Kistler, M. D., Kerr, M., & Dexter, J. 2015

O'Leary, R. M., Kistler, M. D., Kerr, M., & Dexter, J. 2016

Petrović, J., Serpico, P. D., & Zaharijaš, G. 2014, JCAP, 1410, 052

Petrović, J., Serpico, P. D., & Zaharijas, G. 2015, JCAP, 1502, 023

Steigman, G., Dasgupta, B., & Beacom, J. F. 2012, Phy. Rev. D, 86, 023506

Stovall, K., Lorimer, D. R., & Lynch, R. S. 2013, Classical and Quantum Gravity, 30, 224003

Strigari, L. E., Bullock, J. S., Kaplinghat, M., et al. 2008, Nature, 454, 1096

Strong, A. W. 2007, Astrophys. Space Sci., 309, 35

Su, M., Slatyer, T. R., & Finkbeiner, D. P. 2010, ApJ, 724, 1044

Winter, M., Zaharijas, G., Bechtol, K., & Vandenbroucke, J. 2016, ApJ, 832, L6

Yuan, Q. & Zhang, B. 2014, Journal of High Energy Astrophysics, 3, 1

# APPENDIX

## A. ANALYSIS PIPELINE AND DESCRIPTION OF FERMIY TOOLS

We analyze each ROI with a pipeline based on the Fermi package and the Fermi Science Tools. In the following description, we denote with italics the Fermi 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° 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° with respect to the 3FGL position and the 3FGL 68% positional uncertainty is 0.064°.

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°. The result is a list of source candidates with $TS > 16$. On average we detected about six sources with $TS > 25$ per ROI.
<figure><img src="image_13.png" /><figcaption>**Figure 13.** SED of 3FGL J1730.6–0357 (black points) and best-fit PLE parameterization with $1\sigma$ uncertainty band (black line and gray band).</figcaption></figure>

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° 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_{\text{curv}}^{\text{PLE}} = 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°. Specifically we keep the version of the source that is closest to the center of the ROI in which it was found.

## 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 \ge 4$ generally avoided creating spurious clusters around chance spatial coincidences. For *nsrc* = 4, we found that using *dist* = 0.6° 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°, resulted in clusters consisting of chains of sources up to ~4° long in the Galactic plane. Using *dist* ≤ 0.3° 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° and *nsrc* = 4 were also found using *dist* = 0.7°, and only the cluster near the GC was found using *dist* = 0.5°. In light of these studies, we adopted the value *dist* = 0.6° and *nsrc* = 4 and obtain the results presented in Section 3.

## 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.

## D. GENERATING SIMULATED DATA WITH FERMIPY

We use the *simulate_roi* method to simulate the binned γ-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 γ-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 γ-ray counts maps that are statistically identical to those produced with *gtobssim*, which simulates individual γ 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.
<table><caption>**Table 2**

Contents of the 2FIG source list *FITS* table.</caption><thead><tr><th>Contents</th><th>Column Name</th><th>Units</th><th>Uncertainty</th></tr></thead><tbody><tr><td>Source designation</td><td>Source_Name</td><td>...</td><td>...</td></tr><tr><td>Right ascension</td><td>RAJ2000</td><td>[deg]</td><td>...</td></tr><tr><td>Declination</td><td>DEJ2000</td><td>[deg]</td><td>...</td></tr><tr><td>Galactic longitude</td><td>GLON</td><td>[deg]</td><td>...</td></tr><tr><td>Galactic latitude</td><td>GLAT</td><td>[deg]</td><td>...</td></tr><tr><td>Containment radius (68%)</td><td>pos_68</td><td>[deg]</td><td>...</td></tr><tr><td>Containment radius (95%)</td><td>pos_95</td><td>[deg]</td><td>...</td></tr><tr><td>TS</td><td>TS</td><td>...</td><td>...</td></tr><tr><td>TS<sup>PLE</sup><sub>curv</sub></td><td>TS_curv</td><td>...</td><td>...</td></tr><tr><td>Integrated photon flux between E = [0.3, 500] GeV</td><td>Flux300</td><td>[ph cm<sup>-2</sup> s<sup>-1</sup>]</td><td>Unc_F300</td></tr><tr><td>Integrated energy flux between E = [0.3, 500] GeV</td><td>Energy_Flux300</td><td>[MeV cm<sup>-2</sup> s<sup>-1</sup>]</td><td>Unc_Energy_Flux300</td></tr><tr><td>Functional form of the SED</td><td>SpectrumType</td><td>...</td><td>...</td></tr><tr><td>Spectral index</td><td>Spectral_Index</td><td>...</td><td>Unc_Spectral_Index</td></tr><tr><td>Cutoff energy (for PLE)</td><td>Cutoff</td><td>[MeV]</td><td>Unc_Cutoff</td></tr><tr><td>Curvature parameter, β (for LP)</td><td>beta</td><td>...</td><td>Unc_beta</td></tr><tr><td>IEM with which the source is detected</td><td>IEM</td><td>...</td><td>...</td></tr><tr><td>Associated 3FGL source</td><td>3FGL_Name</td><td>...</td><td>...</td></tr><tr><td>Classification of 3FGL source</td><td>3FGL_Class</td><td>...</td><td>...</td></tr><tr><td>Cluster membership (Off. IEM)</td><td>Cluster_Off</td><td>...</td><td>...</td></tr><tr><td>Cluster membership (Alt. IEM)</td><td>Cluster_Alt</td><td>...</td><td>...</td></tr></tbody></table>

**Note.** — When a source is detected with both of the IEM models the reported position, SED parameters as well as the photon and energy fluxes are the ones found with the Off. IEM. We report the TS for curvature, and the SED parameters for the PLE only for PSR-like sources as defined in the main text. For sources with a 3FGL association that was modeled with a LP spectrum we also report the curvature parameter, $\beta$. The IEM column has value “Off”, “Alt” or “Off/Alt”. The *Cluster_Off* and *Cluster_Alt* columns give the index of the cluster to which a given source is associated, if any.

<figure><img src="image_14.png" /><figcaption>**Figure 14.** Photon index $\Gamma$ and energy cutoff $E_{cut}$ [MeV] of PSRs (black points) and for a simulated test population of sources isotropically distributed in the sky (red points). See the text for further details on the SED of these sources.</figcaption></figure>

E. TESTING PSR SELECTION CRITERIA WITH SIMULATED DATA

As discussed in Section 5, ~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_{curv}^{PLE}$ and making spectral fits of blazars and PSRs with a PLE model we have shown that the criteria $TS_{curv}^{PLE} > 9$, $\Gamma < 2.0$ and $E_{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_{cut} [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_{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_{cut})$. Then we select the sources detected with $TS > 25$ and $TS_{curv}^{PLE} > 9$. The
result for the values of Γ and log₁₀(Ecut) is displayed in Figure 14. Also the detected sources satisfying TScurv^{PLE} > 9 maintain a very good separation in the Γ − log₁₀(Ecut) plane. Only 6% of the non-PSR sources detected with TS > 25 have measured Γ < 2.0 and Ecut < 10 GeV and TScurv^{PLE} > 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.