{'fraction_non_alphanumeric': 0.035378779822042614, 'fraction_numerical': 0.018481099412716983, 'mean_word_length': 4.3426767676767675, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Young and Millisecond Pulsar GeV Gamma-ray Fluxes from the Galactic Center and Beyond

Gamma-ray observations have shown pulsars to be efficient converters of rotational energy into GeV photons and it is of wide-ranging interest to determine their contribution to the gamma-ray background. We arrive at flux predictions from both the young (<~ Myr) and millisecond (~Gyr) Galactic pulsar populations. We find that unresolved pulsars can yield both a significant fraction of the excess GeV gamma rays near the Galactic Center and an inverse Compton flux in the inner kpc similar to that inferred by Fermi. We compare models of the young pulsar population and millisecond pulsar population to constraints from gamma-ray and radio observations. Overall, we find that the young pulsars should outnumber millisecond pulsars as unassociated gamma-ray point sources in this region. The number of young radio pulsars discovered near the Galactic Center is in agreement with our model of the young pulsar population. Deeper radio observations at higher latitudes can constrain the total gamma-ray emission from both young and millisecond pulsars from the inner galaxy. While this is a step towards better understanding of pulsars, cosmic rays in the Milky Way, and searches for dark matter, we also discuss a few interesting puzzles that arise from the underlying physics of pulsar emission and evolution.

However, these nearby detections imply that there should be far more gamma-ray pulsars throughout the galaxy that are yet below the threshold for individual discovery by Fermi. That these should lead to a diffuse background has long been expected [e.g., [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24], though the flux level depends on a variety of parameters related to pulsars and the Milky Way (MW) that need to be carefully examined.
Here, we elaborate on the model developed in O'Leary et al. [25], where we demonstrated that that young pulsars (YPs) resulting from supernovae (SNe) near the Galactic Center (GC) and Galactic disk lead to a gamma-ray flux and angular distribution closely resembling the GC GeV gamma-ray Excess [26][27][28][29][30][31][32][33][34][35][36][37][38][39] seen in Fermi data. Though in general this follows directly from SN pulsar kicks, e.g., young pulsars escaping the GC should roughly follow a ρ ∼ r −2 number density profile (smoothed by their distribution of velocities) while getting fainter with time so that the flux density further steepens, the details remain to be precisely resolved.
For example, in Fig. 1 we show the all-sky map of the contribution of young pulsars to the gamma-ray sky in our Fiducial Model [25]. Only a small fraction of the visible point sources are bright enough to be blindly detected as pulsars, and a few hundred more are above the Fermi point source detection threshold. Thousands of faint pulsars blend together to produce a diffuse background that peaks at the Galactic Center with emission all along the Galactic plane.
In this work, we present a comprehensive analysis of the contribution of young pulsars and MSPs to the gamma-ray sky. For the young pulsars we follow the approach of O'Leary et al. [25], but explore the expansive parameter space of pulsar properties, with a particular focus on three different models for for gamma-ray emission. In addition, we present a new (forward) model for the population of MSPs within the Galaxy. With these models, we address the young pulsar and MSP contributions of diffuse GeV gamma rays, the number of unassociated Fermi point sources, cosmic-ray electron production within the inner kpc of the MW, as well as the underlying distribution of unresolved gamma-ray point sources. This latter case is particular relevant to the recent evidence that the excess gamma-ray emission in the Galactic Center is not smooth [40,41]. We also estimate the number of radioloud pulsars in the Galactic Center from both YPs and MSPs, and compare our estimates to current and future radio observations.
In § II, we describe how pulsars and their associated gamma-ray properties evolve with age. We also describe the YRL and YRQ populations of pulsars. In § III, we describe our method for modeling the Galactic populations of young and millisecond pulsars, including their birth rate and birth sites. In § IV, we summarize the three young pulsar emission models and the one MSP Model that we focus on throughout the paper. In § V, we present the Galactic gamma-ray emission arising from our four pulsar models, addressing their contribution to unassociated point sources discovered by Fermi ( § V A) and diffuse emission near the GC ( § V C). In § VI, we address contributions to inverse-Compton emission from around the GC. In § VII we predict the distribution of radioloud young and millisecond pulsars in the GC region and compare this to current constraints. In § VIII we estimate the  resolved and unresolved flux distribution of the pulsars and compare our results to flux distribution derived from one-point photon statistics near the GC. Finally, in § IX we summarize our main results and contributions and discuss a variety of additional implications of our models. In the Appendix, we detail many of the underlying uncertainties in our pulsar models and how they impact expected contributions to gamma-ray excesses. Throughout this work, we consider MSPs as pulsars that have been recycled via mass accretion from a binary companion, and all other pulsars as young pulsars, independent of age. Table I collects much of the nomenclature we use in this work for easy reference.

II. PULSAR PROPERTIES AND EVOLUTION
In this section, we give a general overview of pulsar properties and how they evolve in time. We detail the birth properties of pulsars in our Fiducial Model and our empirical model for the gamma-ray properties of pulsars. Finally we describe the properties of both YRQ and YRL pulsars.

A. Pulsar Evolution
The spin-down power of a pulsar is inferred from the observations of the pulsar period, P , and the change of the pulse period,Ṗ , asĖ where I is the neutron star moment of inertia. For a rotating dipole, the spin-down rate is related to the surface field of the pulsar, B, byṖ where R is the neutron star radius and c is the speed of light [42]. While the dipole spin-down power depends on the angle between the spin axis of the pulsar and the magnetic field axis, α, we do not include this effect in our work, and typical of population synthesis models of pulsars [e.g., 42,43], set sin α ≡ 1.
For an isolated pulsar, observations of only P andṖ cannot uniquely determine the surface field, mass, and moment of inertia (or radius). Observations of binary systems have revealed a wide range of pulsar masses, from 1.2 − 2.0 M [44][45][46]. In this work we will use the canonical values of R ≡ 10 6 cm, M ≡ 1.4 M , and I ≡ 1 × 10 45 g cm 2 . The overall uncertainty in I is about a factor of two. Indeed, other works [e.g., 9], have used a higher value of I. A survey of recent neutron star mass-radius studies suggests that one could justify either choice of I [e.g., 47,48]. As such, most population studies set the moment of inertia, radius, and mass of each pulsar to be fixed, with all of the uncertainties incorporated into the underlying distribution of B and the initial spin period P 0 . We do the same in this work. Knowing B and P 0 , we can directly calculate the time evolution of the properties of each pulsar.
The underlying distribution of surface magnetic fields is one of the most important uncertainties in determining the contribution of pulsars to the gamma-ray sky. The B distribution has been modeled by a log-normal distribution in both population synthesis models of radio pulsars and gamma-ray pulsars, with a mean log 10 B = 12.65. However, these two independent population synthesis methods have found different preferred spreads in the distributions, σ log B . In particular, the B distribution used by Watters and Romani [5] (σ log B = 0.3) for modeling early results of the Fermi pulsar catalog is significantly narrower than the spread preferred from radio surveys [43] (σ log B = 0.55). Since Watters and Romani [5], Fermi has discovered many more pulsars, especially at lowerĖ and older ages.
In this work we use σ log B = 0.45, roughly between the preferred value from previous studies in gamma rays and radio. We find that this is most consistent with the entire population of known pulsars, even if it isn't the preferred value for either the radio discovered or gamma-ray discovered populations. In Appendix A 2, we describe the differences between the radio and gamma-ray selected pulsar populations and directly compare our models with σ log B = 0.3 and σ log B = 0.55 in detail.
In order to describe the complete evolution of a pulsar, we also need to know the initial birth period distribution. We follow Watters and Romani [5], and use a truncated normal 10 34 10 35  distribution, with the spread and mean equal to 50 ms. Likewise, we truncate the distribution below 10 ms. Similar to the B distribution, the preferred birth period distribution of the radio pulsars is in disagreement with the gamma-ray population. In Faucher-Giguère and Kaspi [43], the authors found that the preferred mean was closer to 300 ms, significantly longer than used here. Watters and Romani [5] proposed that this was likely biased by the radio luminosity function of Faucher-Giguère and Kaspi [43]. Using the luminosity function of Stollman [49], Watters and Romani [5] found the radio data was consistent with the lower initial spin periods used in their own work. Nevertheless, here we do not include newly born pulsars with ages ≤ 10 4 yr in our analysis, so few pulsars have spin periods comparable to their value at birth since we use such small initial spin periods.
For each pulsar, we use the initial B and P 0 to evolve the observable properties of the pulsar to its present position and age (see § II A for more details). As the pulsar ages, its rotation slows and its period increases as where we take R ≡ 10 6 cm and I ≡ 10 45 g cm 2 , and t is the pulsar's age.

FIG. 3:
Pulsar gamma-ray spectral fits for variousĖ ranges (shown for the log middle energy in each band), which display a hardening at low energies with decreasingĖ.

B. Gamma-ray Properties and Spectra
The explosion in gamma-ray pulsar detections by Fermi has led to a great empirical improvement in understanding of the emission properties. As the data in Fig. 2 suggest, the gammaray luminosity is described fairly well by the relation with C ≈ 1 [3,5], whereĖ is the total spin-down power of the pulsar inferred via Eq. (1). However, the amount of scatter observed in the luminosity function of pulsars exceeds what is expected from uncertainties in the distance to pulsars alone. The detailed light curves measured by Fermi can in most cases be well fit by assuming gamma-ray emission in the outer pulsar magnetosphere; however, no particular model has a clear advantage in all cases [9]. Not having access yet to the true underlying electrodynamical model, we make general empirically-based assumptions and discuss possible deviations. For our Fiducial Model we set C = 1.3, which is closer to reproducing the mean integrated flux of the known pulsar population ( C = 1.6), but doesn't over produce the number of unknown point sources (see § V A).
In the Fermi data, there is a rough correlation of the fit pulsar spectrum power-law index to the spin-down power, a ∝Ė 0.2 (Eq. 5). For a gamma-ray pulsar population dominated by Geminga-like (Ė ≈ 10 34.5 erg s −1 ) pulsars, the average spectrum of the total emission is harder below 1 GeV than from younger pulsars.
We obtain a description of the gamma-ray spectra versuṡ E by taking spectral fits from the Third Fermi Source Catalog (3FGL) [50] and cross correlating with the list of detected pulsars. This is basic, but attempts to take advantage of pulsars  [3]. This method also results in the phase averaged spectrum for each pulsar, which is required to characterize a population of unresolved sources. We exclude pulsars without a counterpart in 3FGL and those whose spectra are fit only with a power law, as these obviously bias the spectral shape. Generally, all pulsars are fit with cutoff spectra, but there is some lag between discovering a pulsar and updating the model in the catalog. There are also objects, like the Crab, that do not get special treatment in the catalog process and have a strong nebular component.
This leaves 73 young, non-recycled pulsars, which we divide into fiveĖ groups. The top (> 10 37 erg s −1 ) and bottom (< 10 34 erg s −1 ) groups each contain only 5 pulsars, while the others contain ∼ 20, so there is likely some small number effects in these bands.
We normalize each spectrum to its total energy flux and take the median spectrum in the band to obtain the spectra in Fig. 3. There is relatively little evolution in the high-energy (>10 GeV) behavior untilĖ drops below 10 34 erg s −1 . This displays evolution at low energies toward harder spectra with lower cutoffs at lowĖ. These are fit by as given in Table II. At energies below the Fermi range, the spectrum is generally expected to steepen for the pulsars of greatest interest here (e.g., [51]). We note that some uncertainty exists in the bottom range set, which contains only 5 pulsars. Our method using this small number yields a = −1.0 and E c = 0.9 GeV (shown as the dashed line in Fig. 3). AtĖ < 10 34 , we have instead used the empirical spectral relations from Pierbattista et al. [9] to obtain a = −1.0 and E c = 1.4 GeV. As we discuss later, this range arguably stands to benefit the most from improved gamma-ray pulsar data.
The > 10 37 erg s −1 range contains the same small number in the sample set. Such very young ( < ∼ 10 kyr) pulsars tend to be abnormal spectrally. The Crab is the canonical example, although it is not alone in this respect (e.g., PSR B1509). Their gamma-ray efficiencies are very low, with an energy spectrum that can peak in the MeV range. Spin-down index estimates are only available from such very young pulsars, with n < 3 displayed by some (unfortunately, measuring the index of the > ∼ 10 4 yr sample is difficult to impossible). However, such rogue very young pulsars would be rare in our simulations and the most subject to Poisson fluctuations, so we do not attempt to refine the spectrum of these sources any further.
Fermi data indicate that the spectra tend to have more emission at higher energies than Eq. (5) alone suggests [52]. For the brightest pulsars with the most spectral information, such as Geminga, residual high-energy emission is often seen beyond a simple exponential cutoff. This work focuses on the region near 2 GeV, the spectral peak for such pulsars and the GC Excess, and so we do not attempt to account for such surplus emission in our analysis.

C. Beaming and Radio Quiet Geminga-like Pulsars
The population of young gamma-ray pulsars can be split into two classes: pulsars that have observed radio emission (radio loud; YRL hereafter) and pulsars that have little or no observed radio emission (radio quiet; YRQ hereafter). Geminga is the prototypical example of a YRQ pulsar -indeed radio-quiet pulsars are often described as "Gemingalike" -with decades of searches finding no clear detection of a pulsed radio signal [cf. 53].
For a given pulsar, the observable signal depends upon the magnetic field geometry and the angle to the observer. YRQ pulsars are caused by the misalignment of the radio signal (which typically has a small opening angle for young pulsars) from that of the gamma-ray emission. The blind detection of pulsed gamma rays is then the best way to determine if a gamma-ray point source is a gamma-ray pulsar. Conversely, the known population of radio pulsars do not all produce gamma-ray emission visible to us. Approximately half of the young gamma-ray pulsar population is YRQ (not correcting for any observational selection effects).
Beaming of both types of emission plays an important role in determining the fractions of YRL and YRQ pulsars. It is expected that the width of the gamma-ray beam decreases as the efficiency of the emission increases. Likewise, as pulsars spin down, the radio beams also shrink. With smaller opening angles for both beams, it is expected that the YRQ pulsar fraction should increase. As pulsars age andĖ decreases, it is also thought that the magnetic field may approach alignment with the rotation axis [54], resulting in a further increase in the YRQ fraction [5]. Overall, recent Fermi results support both scenarios [55].
We consider for beaming the One Pole Caustic model curves in Pierbattista et al. [6,Fig. 14], approximating the beaming fraction of the YRQ, Geminga-like population with an average beaming fraction described as Similarly, for the YRL gamma-ray pulsars For each beaming model, we randomly assign a fraction, f beam , of the pulsars to be beamed towards the observer, and set the flux of all other pulsars to zero. We increase the total luminosity of the visible population to account for the beaming of the pulsars The total luminosity of the entire pulsar population is thus the same as if we included no beaming (cf., Eq. 4). Our Beamed Model, as we outline in § IV uses only the YRQ beaming fraction, with C = 1.0. We defer to Appendix A 5 discussion on relative differences between these beaming models.

III. GALACTIC PULSAR POPULATION
In this section, we discuss in greater detail the birth rate and birth sites of pulsars in the Milky Way incorporated into the population synthesis model developed in [25], which tracks pulsar trajectories through the MW from birth to the present.

A. Birth Rates
The rate of pulsar birth in the MW should generally follow from the rate of formation of stars sufficiently massive to yield a core-collapse (CC) SN, M * > ∼ 8 M . This is not instantaneous though, since even massive stars take ∼ several-tens of Myr to evolve to explosion [58], so that the SN rate is based on the star formation rate (SFR) smoothed over time.
Typically, the stellar evolution time delay will also smooth over SN locations in the Galaxy, so that fine structure becomes less important for the pulsar population in total, e.g., due to particular stellar clusters in spiral arms. An exception to this would arise for star formation occurring in a location fixed in position over a period of tens of Myr. The GC satisfies this condition, which has important implications.
While the prominent Arches and Quintuplet clusters are both within < ∼ 50 pc of the GC, there are young stars observed throughout the inner ∼ 200 pc Central Molecular Zone (CMZ) region [61]. There is also a sizable, sustained SFR within the central few parsecs [62].
For young gamma-ray pulsars, we are interested in the SN rate over the past few Myr. While the Galactic disk rate should remain nearly proportional to the current average SFR, this does not necessarily hold in the GC. A number of indicators suggest that the CMZ SFR may be episodic on ∼ 10 Myr periods [63], and could have been much higher in the past [64][65][66]. One way to circumvent this difficultly is to count the number of massive stars that are close to explosion (effectively integrating over the past SFR convolved with stellar The differential SN rate as a function of distance from the GC in our young pulsar models, where 7% of the rate is from the CMZ (black line). For comparison, we show two distributions used by Fermi [39] to model cosmic-ray injection: "Pulsars" [56] (red dashed) and "OB stars" [57] (blue dotted). All lines have the same total normalization. Here we use two linear scales for the distance, changing at 0.4 kpc, twice the size of the CMZ. Right: The differential SN rate versus distance from the Sun. evolution). Overall, these considerations make estimating the present SN rate in the GC less direct than simply measuring and converting from the present SFR. At high masses, where evolution is fastest, Rosslowe and Crowther [67,68] estimate the number of Wolf-Rayet (W-R) stars (arising from stars with initial mass > ∼ 25 M ) in the CMZ (≈ 250) and the entire Galaxy (≈ 1200 ± 200), finding a ∼ 20% CMZ fraction. For ≈ 250 W-R stars with an estimated W-R lifetime of ∼ 0.25 Myr [67,68], the implied CMZ SN rate from W-R stars alone is ∼ 10 −3 yr −1 , without even considering any M i < ∼ 25 M lower-mass stars or stars that traverse different evolutionary channels.
Figer et al. [69] used top heavy initial mass functions (IMFs) to arrive at SFRs ∼ 1 − 2% of the total Galactic rate. We are interested only in massive stars. They estimate ∼ 150 O-type stars in Arches, about 3 times this in all the clusters, and 6 times this in the entire region, or ∼ 900 O stars. Taking a normal IMF gets us to ∼ (8/20) −1.35 × 900 ∼ 3100 stars above SN threshold in ∼ 3 × 10 6 yr of star formation, or ∼ 10 −3 SN yr −1 . This is consistent with ∼ 10% of all O / W-R stars being born in the GC.
On the present GC SFR, there are young stellar object (YSO) counts both within and outside of clusters throughout the CMZ. These rates depend on the region under consideration and depend on the amount of time objects are in the YSO phase and estimations of contamination from older stars. Galactic SFR estimates from YSOs tend to be on the low side of the various measures [59], despite concerns for contamination. A SFR of 0.06 M yr −1 within b < 1 • and l < 1 • is derived from free-free emission in [70], though this type of measure also tends to be on the low side in the MW [59].
Going from SFR to SNR also contains a dependence on the IMF and binary fraction of massive stars. A number of studies have found evidence for a top-heavy IMF in the central parsec [71]. Previous indications [72] of a top-heavy IMF in the Arches and Quintuplet clusters have since been revised [e.g., 73], suggesting that the IMF in this environment is consistent with the typical Salpeter [74] IMF at the high mass end. In all cases in the GC low mass main sequence stars have yet to be observed, making definitive estimates of the IMF challenging.
One could bypass stellar arguments by counting supernova remnants, with a substantial number observed [75]. Complications in converting to a SN rate include uncertainty in the duration of their visible lifetimes, Poisson variations in the very-recent rate, and clustering of multiple nearby SNe resulting in only a single superbubble. For example, the XMM superbubbles in [75] may have resulted from multiple SNe and there are likely more yet unobserved SN remnants.
We assume for simplicity that the pulsar formation rate is constant in time and that the GC CMZ has a pulsar formation rate 7 % that of the Galactic rate, 1.5 × 10 −3 yr −1 . Considering the mass-dependent delay time between stellar birth and explosion [58], along with a declining IMF, the SN rate typically evolves in a much smoother manner than an episodic SFR. This is roughly the average of the SFR over the prior ∼ 30 Myr, varying at the tens of percent level over Myr timescales. If there were one or more CMZ star formation outbursts over that period, the SN rate today could easily be a factor of a few higher than inferred from the present SFR (depending on the outbursts intensities and durations). Stars formed in clusters may become dispersed after disruption on a ∼ few Myr timescale [63]. We note that rates do not include exotic outcomes such as magnetar or black hole formation.

B. Young Pulsar Birth Sites and Population Synthesis
We generate the Galactic population of young gamma-ray pulsars by using a two component model to account for pulsars born in the Galactic disk and the CMZ, as in [25], with Fig. 4 showing the distribution of Galactic SNe in our model. As detailed in § III A, we assume that pulsars are born in the disk uniformly through time, with a rate of 2.1 × 10 −2 yr −1 [43,60]. The pulsars start on circular orbits within the spiral arms of the MW using the parameters in Faucher-Giguère and Kaspi [43], who successfully modeled the radio population of young pulsars. Pulsars are formed following an exponential disk with scale radius of 3 kpc and scale height of 50 pc. Each pulsar is placed in one of the four spiral arms, modeled as using the parameters derived from [76]. We assume that stars form equally across the four spiral arms. The distribution of known radio pulsars suggest that the sprial arms are not equally important. We note that we use a slightly different coordinate system than Faucher-Giguère and Kaspi [43], rotated by 90 • . For our Fiducial Model, we assume that the birth rate of pulsars in the CMZ is also constant in time, with a rate of 1.5 × 10 −3 yr −1 , which corresponds to approximately 7 % of the Galactic disk rate. This is still significantly smaller than the observed fraction of Wolf-Rayet stars in the GC (≈ 20 %), suggesting that the total birth rate may be substantially higher. We will explore variations of the CMZ SN rate in Appendix A 1. We assume that the CMZ pulsars are born in the Galactic plane in the CMZ between 20 and 200 pc with scale height of 50 pc. We initialize the orbits with random in plane motion, and follow a surface density profile ∝ r −1 .
We assign each pulsar a kick randomly from a threedimensional normal distribution with a ≈ 408 km s −1 mean velocity [77]. The mean kick velocity is significantly larger than the three-dimensional velocity dispersion of GC stars ≈ 95 km s −1 [78], which redistributes the CMZ pulsars throughout the Galactic Bulge. We integrate each pulsar's orbit in the Galactic potential over its lifetime using GALPY [130] with the MWPotential2014 model [79]. We run our simulations for up to 10 8 yr, significantly longer than the lifetime of most gamma-ray pulsars, to ensure that we include the entire population.
In Fig. 5, we show the resulting surface density profile of gamma-ray pulsars in the MW for the four models that we describe in the next section. The individual points for each model show the bright end of the pulsar population, each with a flux greater than 10 −10 erg cm −2 s −1 (we will call pulsars brighter than this flux limit 'bright' throughout the text). The location of the Sun is shown by the red cross. Since the typical young gamma-ray pulsar has a lifetime that is short compared to its orbital time, the pulsars largely map their birthplaces in the spiral arms of the MW.
In our model as in Faucher-Giguère and Kaspi [43], the Sun is located between the Carina-Sagittarius and Perseus spiral arms, 8.5 kpc from the GC. Unlike Watters and Romani [5], we do not include the local Orion spur in our analysis, which may systematically alter the number of bright pulsars in our simulations. As we are mostly interested in the distant and faint pulsars that contribute to the unresolved background, our results remain mostly unaffected.

C. Millisecond Pulsars
The existence of MSPs as strong gamma-ray emitters was one of the surprises of the Fermi mission. It was thought previously that they might emit through a "pair starved" process, tapping most of the volume above the polar cap for acceleration and emission and thus producing gamma-ray beams with shapes similar to radio beams. Instead, it is clear they have outer magnetosphere emission with plenty of pairs to keep the gaps thin -thus sharp caustic-like peaks. It's not yet understood why the pair formation is so efficient -it's possible that the compact MSP magnetosphere allows multipolar fields to contribute, and/or surface X-ray emission.
In stark contrast to the young pulsar population, MSPs are created long after the SN that produced the neutron star via the spinning-up of formerly-young pulsars by stellar material accreted from a binary companion. Due to their weaker magnetic fields MSPs can remain gamma-ray bright for ∼ 10 9 − 10 10 yr. The contribution of unresolved MSPs to the gamma-ray background has been discussed for both the entire sky [24], as well as the GC [8,27,30,[80][81][82]. Here, we forward model the contribution of MSPs to the GC Excess by using a similar process as we have used for young pulsars.
We assume instead that the re-birth places of MSPs follow the present-day MW stellar distribution. In particular, we consider two components. We assume that MSPs form continuously over time in either the Galactic Disk or the Galactic Bulge. To generate the disk population of MSPs, we assume that the MSPs form on circular orbits in a smooth exponential disk. For each MSP, we assign the surface density of MSPs from an exponential distribution with a scale of 2.7 kpc, and a scale height of 50 pc. To generate the Bulge population, we assume that the pulsars formed on isotropic orbits, from a Γ distribution with an inner profile of k = 1.2, and a characteristic scale of 1.9 kpc. In our Fiducial MSP Model, MSPs are randomly assigned to the disk or bulge populations in proportion to their mass ratio, 68 : 5 respectively [79].
Although MSPs may not receive a substantial kick at birth, dynamical processes heat up the population over time through repeated gravitational scatterings. In our model, each MSP receives a random kick using a three-dimensional normal distribution with σ 1D = 80 km s −1 [83]. Some radio observations suggest that the kick velocities of MSPs may be slightly smaller, with σ 1D ≈ 48 [84]. We have modeled the MSPs with this smaller velocity and found consistent results. We model the initial magnetic field distribution similarly to the ordinary pulsars, with a lognormal distribution, with mean of log 10 B = 8.3 where B has units of Gauss. The dispersion is 0.3 dex. The period distribution at re-birth follows from earlier, with a scale of 5 ms and a minimum birth period of 1 ms. This type of distribution, while following a common methodology used for the young population, differs from the canonical population synthesis that has been used in the past for MSPs [84]. Since we wanted to forward model the MSP population, we found this to be a more natural method. Nevertheless, we find this reproduces the observed distribution of the brightest gamma-ray MSPs. In particular it reproduces the observed pulse period distribution and the latitude distribution of the Fermi pulsars.
In contrast to the young pulsars, the re-birth rate of MSPs is known only indirectly. We therefore generate the Galactic MSP population by continuously populating the galaxy until there are five MSPs with a gamma-ray flux greater than 4 × 10 −10 erg cm −2 s −1 . Above this flux, we assume that the MSP population is completely known from gamma-ray observations. We find that the re-birth rate of MSPs in the MW disk is ≈ 8 ± 3 × 10 −7 yr −1 . This is smaller than that derived in [84], 3.0 +1.8 −1.2 × 10 −6 yr −1 , at less than the 2σ level.

IV. OUR FOUR PULSAR EMISSION MODELS
In this work we focus on four different pulsar models with parameters that we outline here. The first three models are for young pulsars, focusing on the range of uncertainties in their emission properties: our Fiducial Model, Scatter Model, and Beamed Model. These three models all have the same birth properties that we describe in § III. The fourth model that we show throughout the work is for MSPs.
Fiducial Model: Our Fiducial Model is the model presented in O'Leary et al. [25]. We assume that there a one-toone relation between L γ andĖ 1/2 as in Eq. (4) with C = 1.3. We impose a firm cutoff in the pulsar luminosity with L γ = 0 forĖ < 10 33.5 erg s −1 . Systematic variations of the model parameters are described in the Appendix.
Scatter Model: The amount of scatter observed in the luminosity distribution of pulsars exceeds what is expected from the uncertainties in the distances to the pulsars alone. For this reason, we assume that the luminosity of the pulsars follows the same relation as our fudicial model, but includes an intrinsic log-normal scatter. This is similar to studies of the radio luminosity functions of pulsars. For each pulsar we randomly select the luminosity from a log-normal distribution with C = 1.0 for the mean of the log luminosity and a scatter, σ γ = 0.32 dex. If the gamma-ray luminosity is greater than the spin-down luminosity, we then set L γ =Ė. The mean luminosity, under these conditions, is then equivalent to C = 1.3. This is still smaller than the mean luminosity of the observed pulsars, which have C = 1.6.
Beamed Model: A number of theoretical models for pulsar gamma-ray emission suggest that as pulsars age, their emission becomes strongly beamed (see § II C). As such, we have included a model to approximate beaming using the relations we presented in § II C, for the YRQ beamed emission. In this model the average gamma-ray luminosity of all pulsars is equal to setting C = 1.0, but a fraction f beam of pulsars are beamed towards the observer, with a luminosity that is f −1 beam greater than the sky average luminosity. MSP ×20 Model: The possible contribution of MSPs as an important constituent of the GC Excess has been discussed extensively (e.g., [27,30,[80][81][82]85]). In this model, we do not include the contribution of young pulsars, even though we argue young pulsars must contribute significantly in this region. As we will soon see, the assumption just discussed in §III C of a direct correspondence between MSPs and bulge-to-disk mass ratio results in a low diffuse MSP GeV flux.
To see if the Excess might be explained by an overabundance of MSPs near the GC, we use the method described in §III C with an enhancement of MSPs in the Galactic Bulge. Such a boost might result if a fraction of the Bulge formed from the inspiral and disruption of globular clusters [e.g., F. Antonini (private comm.); 85]. We find that to explain the full Excess would require a specific MSP formation rate times that expected from point source counts of the local disk MSP population, where M B and M D are the mass of the bulge and disk respectively. The enhancement factor is degenerate with the assumed bulge mass, e.g., M B ≈ 2 × 10 10 M [86] reduces the required boost to a factor of 6. Overall, the Bulge would need be composed of a few hundred massive globular clusters similar to 47 Tuc, a globular cluster known for its rich MSP population. We include this ∼ 20×-enhanced MSP model throughout our analysis, with the MSP luminosity function chosen to be the same as for young pulsars, with C = 1.3 [3]. In Fig. 2, we plot the three different emission models we use for the young pulsars. We also show the observed gamma-ray luminosities of pulsars, L γ , from the 2PC as a function ofĖ. A number of the YRQ pulsars have no known distances, especially those withĖ < ∼ 10 33.5 erg s −1 , and so are not shown. The MSP emission model is identical to our Fiducial Model with C = 1.3, as shown by the solid blue line.

A. Fermi Point Sources
A requirement of any model of the unresolved pulsar population, is that it, at a minimum, not overproduce the number of observed point sources discovered by Fermi. Ideally, it should also reproduce the bright Fermi pulsar luminosity function.
To determine if a pulsar is detectable as a point source by Fermi, we developed an ad hoc model for the detection of point sources in the Fermi 2FGL [87] and 3FGL [50] point source catalogs near the GC [25]. The 3FGL catalog is complete to much fainter point sources, and provides a more robust constraint on the number of point sources overall. Most previous studies of the GC diffuse background have used 2FGL (or an approximation to it) to remove point sources in analyzing residual emission, while the recent study by Fermi uses a unique point source method [39]. In our analysis of the flux from unresolved pulsars we will use the 2FGL completeness limit when comparing to studies assuming 2FGL and 3FGL in comparing to [39].
We estimate the 50% completeness limit of Fermi derived from the brightness distribution of 2FGL point sources within |b| < 20 • and |l| < 20 • [87] as where b is the Galactic latitude of the point source. We find a similar relation for 3FGL, When comparing directly with the 3FGL point source catalog, we also include an additional constraint that the flux must be greater than 4 × 10 −12 erg cm −2 s −1 . This constraint is always in effect for b ≥ 14.7 • . We should note that these equations were empirically derived by the distribution of known point sources in the GC region, rather than by directly injecting and recovering artificial point sources[cf., 3]. It is worth emphasizing that these are the thresholds for being detected as a point source. Much higher fluxes are required to blindly detect pulsations, which is often required to identify which sources are indeed pulsars.
In Fig. 6, we show probability density distributions of pulsar point sources in our four models. These are calculated by examining the number of point sources in 100 random realizations of each model. In the left panels, we first show the probability density distribution of bright pulsars with flux F > 10 −10 erg cm −2 s −1 . For young pulsars, this is sufficiently bright to be discovered directly from the gamma-ray flux alone, without radio observations. Fermi has detected 31 young pulsars (P > 15 ms) and two MSPs (P ≤ 15 ms) above this threshold. The typical number of bright pulsars in our Fiducial Model is ∼ 20, or only two-thirds of the number detected by Fermi. Our MSP model, on the other hand, does a good job at reproducing the entire luminosity function of the known population.
The middle panels of Fig. 6 show the total number of pulsars above the 3FGL 50% completeness limit within the GC region (|l| < 20 • and |b| < 20 • ). There are 100 unassociated point sources and 10 known young gamma-ray pulsars within this region that are above the same threshold. On average the Fiducial Model has 66 pulsars identifiable as point sources within this region. Assuming that 10 of those are identified as pulsars, on average young pulsars in the Fiducial Model constitute 56 % of the unassociated point sources. Examining spectra of all the unassociated point sources in this region, we conclude that approximately half are consistent with those of pulsars.
The right panels of Fig. 6 show the probability density of total number of pulsars above Eq. 12 within the Galactic Plane (|b| < 5 • ) for our four models. This same region contains 222 unassociated Fermi point sources and 57 identified young pulsars that lie above the same threshold. We find that the Fiducial Model, on average, has approximately 176 pulsars above the threshold, constituting ≈ 54 % of the unassociated population.
It has been argued that the observed distribution of Fermi point sources excludes MSPs from contributing more than ∼ 10% of the observed diffuse background [81]. In the bottom panel of Fig. 6, we show the distribution of observable MSP point sources above the 50% detection threshold of Fermi in the MSP ×20 Model. We find that even when the Galactic Bulge MSP population is twenty times larger than the field population, the MSPs contribute only a small fraction of the point source population. Typically, fewer than 20% of the observed point source population would be MSPs. As such, the point source counts near the GC cannot rule out a MSP origin of the Excess.
Overall the Beamed and MSP Models have the least tension with the number of observed point sources. At the same time, indications from the luminosity function show that pulsars likely constitute more than a third, and close to half, of the unidentified point sources near the GC and Galactic plane.
In this region, some of the measured source spectra may be contaminated by the diffuse background or confused sources, making the spectra tougher to identify. The Fermi-LAT Collaboration [50] notes that the Galactic bulge (within 5 • of the GC), along with the Vela (l ∼ 268 • ) and Cygnus (l ∼ 80 • ) regions, are particularly troublesome, with 7 sources within 1 • of the GC and very bright diffuse emission.
Schinzel et al. [88] found that a large number of unassociated Fermi sources within ∼ 10 • of the Galactic plane did not contain radio sources, and that their latitude distribution is similar to the observed gamma-ray pulsar population. They also note that the average spectral index is ∼ −2.34, and that this differs from the young pulsar population. We note that the faintĖ ∼ 6 × 10 33 erg s −1 pulsar J1705-1906 displays a bump in the gamma-ray energy spectrum, although a power law fit yields an index of ∼ −2.3 [89]. So steep power law spectra may not be sufficient to rule out faint young pulsars.

B. The Unresolved Population
We generate maps of the contribution of pulsars to the gamma-ray sky by convolving the gamma-ray flux of each pulsar with the point spread function (PSF) of Fermi. This is done twice, instead of once, to reproduce the reconstructed gamma-ray maps used in the literature [e.g. 35,37], who smooth the observed distribution of gamma rays by the PSF or a Gaussian filter after the image has effectively already been smoothed when measured by Fermi. We calculate the PSF using the energy dependent model of Ackermann et al. [90], who determined that the Fermi PSF is best reproduced by using the sum of two King functions (see their Eq. 5). This model reproduces both the core of the flux distribution as well as the much broader tails. We refer the reader to [90] for details. We focus on the contribution of pulsars at 2 GeV, near the peak of the pulsar spectral contribution.
In Fig. 1, we show a Mollweide projection of gamma-ray flux from young pulsars in the entire MW. Overall, the pulsar emission follows the Galactic plane, with a peak near the GC due to young pulsars formed in the CMZ, which have not been included in previous analyses of gamma-ray pulsars [cf. 8].

C. Diffuse GeV Emission near the Galactic Center
We reconstruct the diffuse emission from the unresolved pulsar population near the GC by first removing all pulsars with fluxes above F ps,2FGL . Although many more point sources are known in the GC region, most current studies mask only those sources discovered in 2FGL [cf. 39]. At present, we do not mask the brightest pulsars (which usually only covers 95 % of the flux), but completely remove them, thus our results should be considered a lower bound on the flux in these regions.
In Fig. 7, we show our reconstructed maps of the GC region for our four models, with the Galactic plane masked (|b| < 2 • ). The Fiducial Model also includes an overlay of the different regions of interest. The High Latitude pulsars are those above the diagonal lines with |b| ≥ |l|. Unlike in the rest of this work, this map shows the effect of the finite number of photons collected. The three young pulsar models show emission in all latitudes, but the flux concentrates around the GC and Galactic plane. The MSP ×20 Model is more spherically symmetric at large angles, as the pulsars are formed using a symmetric Galactic Bulge model.
In Fig. 8, we plot the angular profile of gamma-ray emission from unresolved pulsars near the GC. The Fiducial Model is  as constructed in O'Leary et al. [25]. The lines show the total flux in annular apertures around the GC with width of one degree for each of our four models. The black points show the results from Daylan et al. [35]. The shaded blue regions show the results from Calore et al. [38]. Overall, our the models follow the observed Excess of gamma rays from the GC region out to ≈ 20 • . The three young pulsar models, we note, are not tuned to match amplitude of the GC excess. However, the total flux is uncertain by at least a factor of two (shown in the gray region), due to underlying systematic uncertainties, which we explore further in the Appendix.
In the bottom panel of Fig. 8, we show the unresolved contribution of our Fiducial MSP Model (green solid line), which only reproduces a fraction of the observed diffuse emission near the GC. We estimate approximately 20% of the diffuse emission comes from MSPs in the Galactic disk. Fig. 8 also shows the results of enhancing the bulge MSP population by 5× (dotted) and 20× (dashed) the number of MSPs per stellar mass as the Galactic disk. As we discussed previously, the MSP model requires that the bulge have an enhancement of MSPs relative to the Galactic disk.
Overall, we find that a population of young pulsars born near the GC can roughly reproduce the amplitude, spectrum, and angular distribution of the Excess quite directly. Such a signal has a strong component in the Galactic Plane, consistent with the excess of gamma rays observed in the plane [91], as we discuss further in § A 1. MSPs can also reproduce the spectrum and angular distribution of the signal. The spherical morphology imparted by the model appears to be a better match to the observations by Calore et al. [37]. However, as we described in § III C, this amplitude can only be attained if there is a substantial enhancement of MSPs that extend throughout the bulge, in proportion to the stellar mass. We explore the latitude dependence of the signal in Appendix B. We compare the overall amplitude of the fluxes in the GC region from the Fiducial and ×20 MSP models to Excess spectra from The Fermi-LAT Collaboration [39] in § VI.

VI. INVERSE COMPTON FROM THE INNER KPC
The best fitting models in the detailed recent study for the inner ∼ 1 kpc of the MW by The Fermi-LAT Collaboration [39] include a much higher normalization for the inverse-Compton (IC) component (6 − 30×) than their baseline models, possibly due to a larger cosmic-ray electron population or radiation field intensity. For such an E γ = 1 − 100 GeV IC flux, considering photon backgrounds within 1 kpc, requires e ± of E e ± = 5 − 500 GeV. This energy range corresponds to the positron excess measured by PAMELA [92] and AMS-02 [93,94] most easily ascribed to nearby pulsars (e.g., [95][96][97]).
While the gamma-ray efficiency of young pulsars increases with age, most of the spin-down power is still carried off by a pulsar wind containing energetic electrons and positrons. Aṡ E ∝ t −2 , the youngest pulsars in the population dominate the e ± energetics, so this is not sensitive to any lowĖ cutoff and is the same for all three of our young pulsar models.
We show the MW distribution ofĖ in Fig. 9, which we see is highly concentrated within 1 kpc of the GC and along the Galactic plane (although propagation should stretch this along the local magnetic field direction [98]). As stated in [25], we then expect a substantial contribution to the inner Galactic IC background above that usually assumed (see the left panel of Fig. 4). There could also be an increased contribution from SN remnants. Both of these would be more relevant for the cosmic-ray electrons than for protons, due to the higher rate of e ± energy losses making local production more important, especially at high energies due to a harder spectrum (maybe even more at the GC than indicated by PAMELA/AMS, since at Earth the lower density of sources may result in a substantial variation). This is consistent with the much lower flux attributed to neutral pion decay in the inner ∼ 1 kpc by [39].
We find that a lognormal distribution describes the total summedĖ power within 1 kpc, with log 10 Σ iĖi ≈ 38.2 ± 0.1 erg s −1 for young pulsars, and ≈ 37.5 ± 0.2 erg s −1 for the 20×-enhanced MSP scenario. We note that though we do not include pulsars with age less than 10 4 yr in this calculation, a single Crab can reachĖ > 10 38 erg s −1 (at the far right of Fig. 2), though no such object is necessarily present.
While details of cosmic rays will have to be sorted out with more precise modeling, we can simply scale the e ± luminosity toĖ to estimate their effect. SinceĖ ∝ t −2 , the youngest pulsars (i.e., those near the disk) are most relevant in this regard, and we content ourselves with an equilibrium, one-zone cal- culation [99] using a photon background appropriate for this region [100] and a plausible range of magnetic fields within 1 kpc of B = 10 − 50 µG. For E e > ∼ 10 GeV, the e ± cooling time is then t c < ∼ 10 6 yr, so variations below this scale are smoothed, while this is still a relatively short time for diffusive escape from 1 kpc. Since t c decreases roughly as E −1 e , higher-energy IC emission may be more localized. Fig. 10 shows the IC gamma-ray flux using a young pulsar e ± source spectrum motivated by positron excess measurements, dN/dE e ± ∝ E −1.8 e ± exp[−E e ± /500 GeV], and taking L e ± = 0.5 Σ iĖi , for the range of B. These are shown with our Fiducial composite High-Latitude young pulsar spectrum, scaled by effective area for even comparison with the 15×15 • Fermi NFW component flux [39], and IC spectra for "Pulsar Index-Scaled" and "Pulsar Intensity-Scaled" IEMs [39].
Although young pulsars do not necessarily account for all of the IC emission, and 100% ofĖ need not go into e ± , we see they should help account for the enhanced IC without needing to invoke an enhanced photon background or other ad hoc assumptions. There will also be some GeV-TeV pulsar wind nebula emission (even at old ages as seen around Geminga [101] from such a wind [95]) and variations owing to the local photon field -such as near the GC, Arches, and Quintuplet star clusters -and magnetic field strength [99,102]. Inverse-Compton emission from the disk, on the other hand, should already mostly be accounted for in models assuming a radialdependent source density (e.g., [39]; see again left panel of Fig. 4). While the 20×-enhanced MSP model yields a lower summedĖ within 1 kpc, a less concentrated profile could lead to off-plane IC emission depending on MSP wind properties.

A. Young Pulsars in Radio
Radio surveys of the central 50 pc of the GC probe the underlying distribution of radio pulsars within that region, and hence star formation on a ∼ 10 7 yr timescale. To estimate the total number of radio pulsars in this region we use the same simulation parameters as our Fidicual gamma-ray simulations, except we follow all pulsars down to the radio death line [43,103]. We use the empirical pulsar radio luminosity function from Stollman [49, Eq. 5] to estimate the radio flux from each pulsar. This luminosity function is nearly identical to the one used in Faucher-Giguère and Kaspi [43] forĖ < ∼ 10 33 erg s −1 . More energetic pulsars, however, all follow the luminosity distribution of pulsars withĖ ≈ 10 33 erg s −1 . The radio emission from pulsars is strongly beamed, and so only a fraction of the pulsars are visible to us. Following Tauris and Manchester [104], we assume that only a fraction of pulsars are beamed in our direction f radio = 0.09 [log 10 (P/10 s)] 2 + 0.03.
With these assumptions we find that there should be ≈ 5 − 15 radio pulsars beamed towards us that could have been detected by the GC pulsar survey of Deneva et al. [105], depending on if we include the inner 7 pointings where the sky is brightest. This is consistent with the 5 pulsars observed within the inner ∼ 50 pc, given the uncertainty in the spatial distribution of star formation in our model, and variable sky brightness. Indeed, one could use these radio detections to place a lower limit on the gamma-ray contribution of young pulsars in this region. Searches for compact radio sources can also limit the total number of pulsars in the GC. Lazio and Cordes [106] estimate that no more than 10 4 − 10 5 total pulsars may exist within 2 • of the GC. In our simulations ≈ 1 × 10 4 pulsars are found within the survey area. The pulsed radio searches in this region from Deneva et al. [105] place a much tighter constraint on the overall population.
Future radio searches above the Galactic plane would provide an important constraint on the population of young pulsars. The central ≈ 50 pc are partially obscured by a strong scattering screen that smears out signals on short timescales.
Because of this, many deep surveys have searched at high frequencies where the impact of scattering is less, but the pulsars are intrinsically fainter. The strongest regions of the scattering screen extends only to |b| < ∼ 0.2 • [107]. Above the Galactic plane, the expected amount of scatter is significantly less, although the pulsar density also decreases with latitude. In Figure 11 we show the number density of pulsars in 1 sq. degree regions above and below the Galactic plane with |l| ≤ 0.5 • . We find that near |b| ≈ 1 • , the pulsar density is approximately 45 deg −2 with fluxes above 50 µJy at 1400 MHz. This drops to 25 deg −2 at b = 2 • . The majority of these, ≈ 80 %, originate within the CMZ. The remainder are pulsars from the foreground disk.
The successful detection of pulsars in the CMZ is in stark contrast to the non-detections from deep radio surveys of the central parsec. We estimate, using the same forward modeling and radio luminosity function, that there should be ≈ 10 − 20 detectable pulsars in the deepest surveys of the central parsec for a SNR of 10 −4 yr −1 . This is slightly larger but consistent with the estimates by Dexter and O'Leary [108] who used purely empirical data and relationships to estimate that ≈ 10 radio pulsars should have been detected in this region. Both estimates are at odds with the deepest surveys of the central parsec [109][110][111][112][113], which have found no ordinary, young radio pulsars. This population of pulsars is not included in any of our gamma-ray models.

B. Millisecond Pulsars in Radio
Deep pulsed radio searches may be able to constrain the population of MSPs in the Galactic Bulge. To date, the deepest searches have focused within ≈ 1 pc of the GC, where the sky brightness is high and scattering from the interstellar medium (ISM) may smear out the pulsations. Even these deep radio searches are sensitive to only the brightest few MSPs known in the whole galaxy. Above the plane, the ISM scatters the pulsar signals by orders of magnitudes less. As such, deep pulsar search near ∼ 1 GHz may be able to detect the absolute number of MSPs in the bulge, and constrain the total gammaray background produced by this unresolved population.
MSPs are typically fainter than the YRL pulsar population, even for pulsars with the same spin-down luminosity. For this reason, we use a different radio luminosity distribution for the MSP population. Following Calore et al. [114], we use the best fit model of Bagchi et al. [115] to estimate the flux of MSPs in the Galactic Bulge. We find that the surface density of radio bright (S 1400 > 10 mJy) MSPs is approximately 20 deg −2 at 2 • above the Galactic plane (note that this flux limit is five times lower than quoted above for YRL pulsars). Within 2 • , we do not reproduce the total flux of the gammaray excess. Detailed estimates of the radio detectability of MSPs have recently been carried out by Calore et al. [114] using a complementary approach to model the underlying pulsar population.
As opposed to young pulsars, where there is overwhelming evidence for a substantial radio-quiet population, there is yet no firm evidence for radio-quiet MSPs. While the predicted values for, say, 47 Tuc are larger than the known radio MSPs, this only arises after assuming an efficiency. There are still high-latitude point sources that look like MSPs, though ruling out a radio counterpart is very difficult due to the variety of possible signals.

VIII. (TEN) THOUSAND POINTS OF GAMMA RAYS
Lee et al. [40] use non-Poissonian photon statistics of pixels to derive the flux distribution of unresolved point sources in the GC region on top of a diffuse emission model. They focus on the region above the plane (|b| > 2 • ) and within 10 • of the GC (though practically the large number of point sources excludes most of the < 5 • region), using 1.893 − 11.943 GeV gamma rays. The authors found that the entire Excess could be accounted for with a population of unresolved point sources. In this section, we compare the flux distribution of our four models to see if a population of pulsars is consistent with the unresolved sources inferred in [40].
Only two out of 40 MSPs haveĖ >Ė peak . A much larger fraction of known young gamma-ray pulsars exceedĖ peak , 41 of 77 [3], though only 8 of these are above the plane with |b| > 2 • . Also, the young pulsars have a median characteristic age of only 2 × 10 4 yr. In our models, this population of pulsars is too rare to reproduce the Excess for the GC SFR we have chosen. Within the region of interest of [40], we have only two such pulsars in one realization. Beaming of the gammaray emission at lowĖ ( § II C) may eliminate this discrepancy. However, the discrepancy with the maximum luminosity of the point sources inferred in Bartels et al. [41] is even greater. This is especially the case in the model of Brandt and Kocsis [85], which is dominated by an old (τ > ∼ 10 9 yr) population of pulsars with lowĖ.
In Fig. 12, we directly compare the flux distribution of our four models with the inferred population of unresolved point sources in [40]. The shaded region in Fig. 12 shows the 68% contours for the reconstructed flux distribution of resolved and unresolved point sources in the region of interest. The solid lines show the distribution of point sources in our simulations with a hard cutoff (death-line) atĖ = 10 33.5 erg s −1 , as in our Fiducial Model from [25]. The dashed lines show the same results if we institute a hard cutoff at 10 34 erg s −1 .
In all of our models, the peak of the flux distribution is at lower fluxes, typically corresponding to the flux of pulsars near the gamma-ray death line we imposed. Increasing the flux of all of our pulsars does not cleanly resolve the discrepancy. For the dN/dF distributions shown, increasing the flux of all pulsars shifts the lines to the lower right. This can be seen by comparing the Fiducial Model to the Beamed Model, where the observed lowĖ pulsars have much larger luminosities (see Fig. 2). The amplitude of dN/dF for the beamed model is in better agreement with the reconstructed flux distribution.
It is a challenge to directly compare the results of our analysis with [40], since we are reporting the underlying flux distribution from our models, while the derived flux distributions from [40] are created using templates for the point source distribution. One major assumption of the templates, is that the spectrum of the point sources does not depend on location or brightness. Our pulsar models violate both assumptions. We also note that the dN/dF distribution in [40] does not appear to vary greatly for radial profiles ranging from r −2.2 to r −2.8 , so an even a steeper GC pulsar profile may be consistent with their result, and it may preferentially pick up sources nearer to the plane which are intrinsically brighter. The break seen in the flux distribution may be caused by a combination of two or more unresolved point source distributions (e.g., SN remnants, or both YPs and MSPs). Indeed, we estimate that about half of the resolved 3FGL point sources in this region are from pulsars, leaving the other half to other classes of gamma-ray sources. More detailed modeling may allow us to better recover the properties of the underlying point sources.
The locations of such possible unresolved GC point sources are one place to look for radio pulsars. However, some fraction may be YRQ pulsars and the "hot spots" themselves may not be great proxies for the actual source position. Systematic uncertainties of the background in the GC region may mean that the effective localization of the point source by Fermi will be larger than the beam width of a large radio telescope. A beam width of ∼ few arcmins would require a large area mosaic around such regions, likely requiring a major survey. A better bet could be continuum imaging with one of the new small-dish wide FoV arrays, e.g., ASKAP, to look for correlations between point sources and Fermi hot spots.

IX. DISCUSSION AND CONCLUSIONS
Pulsars are ubiquitous in the Milky Way and studies using Fermi have revealed that pulsed gamma-ray emission is common. However, other than subtracting the nearest and brightest objects as point sources, young pulsars have been largely neglected in diffuse studies near the Galactic Center, despite a large number expected to result from supernovae in the Central Molecular Zone [25]. Furthermore, a large population of millisecond pulsars might also be present throughout the bulge from disrupted globular clusters [85]. In this work we have explored a variety of emission models of pulsars for both young and millisecond pulsars. We have used these models to predict both the expected diffuse emission of gamma rays from the unresolved populations of young and millisecond pulsars, as well as their contribution to the unassociated Fermi point sources in the 3FGL catalog.
We have shown that young pulsars should contribute a significant fraction ( > ∼ 50 %) of the GeV gamma-ray excess observed towards the Galactic Center region using a variety of emission models. Within 5 • of the Galactic Center this emission is effectively spherical, with a radial profile that strongly matches that inferred by a number of groups to within the point spread function of Fermi. Though a typical luminosityweighted average for the young pulsar population is dominated by the softer, brightest pulsars confined to the plane, the summed total young pulsar emission becomes harder at larger angles. Our composite spectra is still somewhat softer than some Excess claims. However, none of our pulsar models are adequately described by NFW profiles, so we do worry that studies that do not include pulsars self-consistently are oversubtracting fluxes at low energies.
Millisecond pulsars may also be an important component in this region, especially at high latitudes, |b| > ∼ 5 • . We forward modeled MSP evolution throughout the Milky Way in a manner similar to the young pulsars. We found that the Galactic bulge must have 6 − 20 times the number of MSPs per stellar mass than the Galactic Disk to explain the GeV excess amplitude outside of 2 • . A few hundred massive, disrupted globular clusters would be sufficient to explain the higher density of pulsars [85]. This model, which largely follows the stellar distribution of the bulge, also matched the morphology of the observed excess outside of 2 • of the Galactic Center. Within 2 • , we require an even higher concentration of millisecond pulsars per stellar mass than explored here.
We compared predictions of our young and millisecond pulsar models to the observed distribution of Fermi unassociated point sources in the Galactic Center region, as well as throughout the Galactic Plane. For the young pulsar population, we expect approximately one-third to one-half of the unassociated point sources near the Galactic Center to be young foreground pulsars depending on the emission model and supernova rate. Contrary to previous statements [e.g., 81,116], we find that millisecond pulsars from the Galactic Bulge do not contribute significantly to the point source population inferred in the Galactic Center region. Old millisecond pulsars are not bright enough to be detected at the distance of the Galactic Center.
Radio observations throughout the Galactic Center region can potentially constrain the contribution of gamma rays from both young and millisecond pulsars. In fact, the young pulsar model discussed here is already consistent with the number of pulsars discovered within the Central Molecular Zone ( § VII). Within the uncertainties of emission models for gamma-raybright pulsars, this suggests young pulsars are a promising source of the gamma-ray emission within this region. Outside of ≈ 0.4 • , no radio surveys have been sufficiently deep to differentiate between a young pulsar origin and a millisecond pulsar origin of the gamma-ray excess. Further observations at higher latitudes, whether deep pulsar searches or surveys for point sources, will further constrain the relative contribution from young pulsars and millisecond pulsars. We found that the young pulsar population can be well characterized by surveys that reach flux limits near 50 µJy at 1.4 GHz. Millisecond pulsars are intrinsically fainter than the young pulsar population, and require surveys reaching ∼ 10 µJy to detect more than a few per square degree.
Overall, we find good agreement between our pulsar models and the observed excess distribution. While neither model exactly reproduces the inferred luminosity function of the unresolved point sources near the GC (see § VIII), we expect that a significant fraction of flux from unresolved point sources come from young and/or millisecond pulsars. One uncertainty is why, in this particular region, the contribution of pulsars is so much more prominent than the rest of the galaxy. Typically, we find young pulsars contribute ≈ 10 % of the to-tal GeV flux of the Galaxy. In the GC region, this is closer to ≈ 30 %. One possibility is that the strong outflows from this region quickly advect cosmic rays away [64,117], which could explain why the CMZ is also fainter than expected in non-thermal radio emission [117].
Our estimates suggest that an important implication of YPs, as opposed to dark matter especially, is a large e ± luminosity above that typically assumed with the inner kpc of the MW (see left panel of Fig. 4). These will lead to gamma-ray emission through inverse-Compton scattering at a level depending on the amount of spin-down energy converted to energetic e ± , likely a large fraction with a hard spectrum, and the ambient magnetic field strength relative to the photon field. We found that an IC flux near to that inferred by Fermi [39] can arise.
There may be other contributions to diffuse emission beyond that associated with gas maps or pulsars (possibly including dark matter). One can also consider primary electron acceleration by supernova remnants, although the escaping flux is uncertain and likely softer. The much lower gamma-ray flux attributed to neutral pion decay in the inner ∼ 1 kpc region by The Fermi-LAT Collaboration [39] does suggest that secondary e ± are far insufficient to account for the IC flux.
We defer more detailed relations to dark matter to [118], where we examine a template for dark matter searches (as well as details of the central degrees of the GC). Obviously we only measure one realization of the pulsar distribution in the MW, which is rather different than for standard assumptions of the smooth dark matter halo (although this too likely depends on peculiarities of Galactic history), although it shows a way forward to extract signals when pulsar backgrounds are relevant.
One could utilize our models with Fermi data, account for point sources in a consistent manner, and cross correlate ob-servations of dwarf galaxies with the GC to find a signal common to both populations. Recent dwarf galaxy limits using Fermi Pass 8 data bite into the parameter space of 30-40 GeV annihilation channel into bb [119,120]. A higher dark matter mass would skirt such present limits (which young pulsar emission automatically satisfy). A higher-energy component should be least noticeable at the GC, becoming gradually more relevant at larger latitudes (depending on the dark matter density profile). Signs of this may be evident [38], although we have also shown that inverse-Compton emission from young pulsars can be at least as luminous than the pulsed emission and likely extends above the disk.
Contrary to a common assumption that the GC GeV excess is already seen with sufficient statistics and that settling the issue will remain entirely in dealing with systematic effects, we suggest that headway can be made with more Fermi data by just finding more sources and characterizing more pulsars. Since these are not all from an invisible population, it is realistic to expect improvement by collecting more photons to better understand the signal and the pulsar physics of those objects that cannot be seen individually at the GC. could themselves inflate bubbles [see, e.g., 64,65]. Indeed, one possible shortcoming of the young pulsar interpretation, is that the gamma rays appear to be suppressed in the GC region compared to the star formation rate.
As we discussed in detail in § III A, the birth rate of neutron stars in the CMZ is fairly uncertain, and likely has changed over time. In our Fiducial Model we assume that pulsars are born at 7 % of the Galactic pulsar formation rate, as this is consistent with the varied estimates of the SNe rate in the CMZ. Nonetheless, the actual rate may be as high as 20 % if we use the WR population as a tracer population for SN.
In Fig. 13, we show the contribution of pulsars to the diffuse background near the GC for three different CMZ SN rates. The solid red line shows our Fiducial Model with a contribution of 7 % of the Galactic rate. The dashed red line shows the contribution from the Galactic Disk only, with no contribution from the CMZ at all. Comparing these two lines, we see that the CMZ contributes only a small fraction of the diffuse flux outside of > ∼ 5 • ; foreground pulsars from the Galactic disk contribute the majority of the flux at higher longitudes. Finally, the dashed line shows the contribution of pulsars if the CMZ SNR is 20 % . We see this is sufficient to reproduce the entire observed Excess. Nevertheless, we expect contributions from other sources in this dense region (e.g., MSPs). In principle, this shows that the CMZ SN rate can not be significantly larger than ∼ 4 × 10 −3 yr −1 , otherwise we expect the unresolved pulsar emission from this region to exceed what has been observed.
Here and in [25], we did not include star formation from the central parsec of the Milky Way, only populating the region within 20 − 200 pc of Sgr A * . Some aspects have been discussed in [99,102] and we will devote even more detailed consideration of this interesting region elsewhere, but do fo-cus on their radio characteristics in § VII. Also, Watters and Romani [5] included a thick disk component in their models with a 10% runaway component and 500 pc full-width-halfmaximum. While we have not included such a process in our results, it would shorten travel times for such YP to reach large angles on average. We take our model to be more conservative in estimating the gamma-ray emission at high latitudes.
Relative to MSPs, we have much better prospects for understanding from first principles the young pulsar luminosity density profile. Further, we also have a much better expectation for the distribution and luminosity of individual pulsars. Since at later ages beaming may become more relevant, the pulsars that are visible to us will be brighter than the average we have assumed above, improving the prospects for detecting this population with more observations.

Magnetic Field Distribution
There is a general disagreement between the intrinsic properties of pulsars that are assumed in population synthesis studies of radio pulsars and gamma-ray pulsars. Although both use the same mean magnetic field strength, radio pulsars are typically modeled with a larger variance σ log B = 0.55 [43] than the gamma-ray population σ log B = 0.3 [5]. We examine here these two populations in more detail, with a focus on the observational selection effects that impact the gamma-ray pulsar population.
When comparing our population synthesis models to the observed population, we focus on the bright pulsars with flux > ∼ 10 −10 erg cm −2 s −1 . This population of pulsars should be complete in the Fermi data. Using the 2PC, we find that the variance of the observed population is σ log B ≈ 0.31. This is the value used by Watters and Romani [5]. However we find that observational selection effects cause the observed variance of the bright pulsars to be smaller than the intrinsic variance, typically by 20 %. This is a natural consequence of pulsars with large magnetic fields evolving more rapidly to lower luminosities. This suggests that the intrinsic value of σ log B is likely closer to 0.4 than 0.3.
In Fig. 14, we compare the P -Ṗ distribution of the pulsars with three different initial magnetic field distributions. With red circles, we show the bright observed gamma-ray pulsars with flux > ∼ 10 −10 erg cm −2 s −1 that are the basis for most of our comparisons. The black points show the entire gamma-ray pulsar population for each of our models, and the blue points show the simulated pulsars with flux > ∼ 10 −10 erg cm −2 s −1 . Each model shown is selected to be a representative sample with a Kolmogorov-Smirnov (KS) p-value near the mean of all the simulations. More generally, for each model we can evaluate the observed distribution of σ log B for the brightest pulsars and determine the fraction of realizations that reproduce the observed distribution. For our Fiducial Model with σ log B = 0.45, we find 15 % of our simulations have a simulated variance less than that observed. With σ log B = 0.40, we find 30 % of our simulations have a inferred dispersion less than the observed dispersion. In contrast, only 2 % of models with σ log B = 0.55 have smaller variance.
In this work we do not directly model the observational selection effects of the radio discovered pulsars, although we expect the bias to be towards lower values of σ log B . Nevertheless, our initial conditions should generally be informed by the properties of the radio population. We estimate σ log B using the Markov chain Monte Carlo ensemble sampler emcee [124], where we iteratively removed 3 − σ outliers between each estimate until we reached a converged value. We use a uniform prior for σ log B . In Table III, we report our estimates of σ log B for both the gamma-ray pulsars and radio pulsars in bins of varyingĖ, as well as the bright population.
Overall, the population of radio pulsars have larger observed σ log B compared to the gamma-ray pulsars, even when removing the outliers. As we incorporate more pulsars into the radio analysis, most of which do not emit in gamma rays, the distribution begins to approach that inferred by Faucher-Giguère and Kaspi [43], but remains smaller. This is likely  15: Intensity of the pulsar background at 2 GeV as a function of angular separation from the Galactic Center, compared to the GC Excess intensity at 2 GeV inferred in Ref. [35] (black points and error bars) and Ref. [38] (shaded blue regions). The solid red line shows the pulsar flux from our Fiducial Model. The other red lines assume a pulsar σ log B = 0.55, the value that best describes the radio population of pulsars, and σ log B = 0.30 as used in [5]. due to similar observational selection effects. One way to alleviate the discrepancy between the gamma-ray pulsar properties and the radio pulsar properties is to allow the initial spin period to depend on the surface magnetic field of the pulsar. Alternatively, the differences observed within the radio population may depend on some true physical effects that we do not model here. For example, some population synthesis models allow the surface field to decay with time. There is still no compelling evidence from radio surveys [43] or in this work that requires us to include surface field decay. In Fig. 15, we show the diffuse gamma-ray flux near the GC for four models with σ log B = 0.3, 0.40, 0.45, and 0.55 dex. We see that the estimated flux from pulsars drops by approximately a factor of two when going from σ log B = 0.55 to 0.30, which correspond to the preferred values from previous analyses of radio and gamma-ray discovered pulsars, respectively. Pulsars with large magnetic fields evolve more quickly and emit gamma rays for a shorter time than pulsars with weaker fields. By having a broad distribution of magnetic fields, such as those implied by radio observations, many more pulsars with weak fields are born. These pulsars can travel farther within their gamma-ray emitting lifetime, where they contribute to the diffuse gamma-ray background.
Overall, we find that σ log B = 0.45 seems to match the properties of the combined population of gamma-ray and radio selected pulsars better than 0.3. It is only slightly larger than most likely intrinsic variance of the gamma-ray population, σ log B = 0.37.
The large systematic difference between the observed gamma-ray pulsars and the radio pulsars in the ATNF catalog [123] really drives our use of a larger dispersion.
The pulsar formation rate and the distribution of magnetic fields are both degenerate with how the pulsars evolve with time. In particular, discrepancies could arise from uncertainties in the moment of inertia in addition to simplifications of the rotating dipole model. The star formation tracers probably have some systematic uncertainties in their normalization as well. One way of getting around this is to populate the galaxy with pulsars until one arrives at the number of objects observed (i.e., obtaining an effective pulsar birth rate that may large compared to the SN rate [6,43]). In order to keep the number of free parameters to a minimum, we choose not to model these directly for the young pulsar population. . The green dashed line shows the contribution of pulsars at high Galactic latitude (|b| > |l| and |b| > 2 • ). These are compared to the GC Excess intensity at 2 GeV inferred in Ref. [35] (black points and error bars) and Ref. [38] (shaded blue regions).

Gamma-ray Death Line
For the heuristic value of C = 1.0, the gamma-ray efficiency of pulsars rises to 100% asĖ approaches 10 33 erg s −1 . For our Fiducial Model C = 1.3, the pulsars reach unity efficiency at a higher energy,Ė ≈ 1.7×10 33 erg s −1 . We counter this effect in our Fiducial Model by setting the gamma-ray flux to zero for pulsars withĖ < ∼ 10 33.5 erg s −1 , which corresponds to a maximum efficiency of 73 %. It is, of course, also implausible that the emission would die just as the efficiency approaches unity. Nevertheless, we first analyze our results looking at fixed cutoffs for a variety ofĖ. In Fig. 16, we show the amount of diffuse flux near the GC for a variety of energy cuts from 10 33 erg s −1 to 10 35 erg s −1 . Later, lacking the true model of pulsars, we take general, empirical approach that can be scaled for a given magnetospheric emission model.
At present, the YRQ gamma-ray pulsar PSR J2139+4716 has lowest inferredĖ∼ 3 × 10 33 erg s −1 [125]. There is indeed much room at the bottom and the luminosity distribution stands to improve [3,89,126]. This population tends to be YRQ, which means these pulsars must be discovered through dedicated blind searches of the gamma rays. Presently, we find that the observed distribution of the brightest pulsars makes it hard to exclude a deathline near 10 33 erg s −1 . That is, when looking at pulsars with luminosities between 10 33 erg s −1 and 10 33.5 erg s −1 , 44 % of our simulations have none with fluxes greater than 10 −10 erg cm −2 s −1 . This is consistent with the observed population.
Allowing the pulsars to approach unity efficiency results in a boost of flux from unresolved point sources near the GC. Primarily, magnetic field curvature and spin-down power determine when gamma-ray emission in young pulsars completely shuts off, i.e., when the pair cascades are not able to sus-tain themselves and can no longer produce enough pairs to short out the accelerating gaps. This is also when the gap size increases (w) so that the emitting surface moves onto lesscurved field lines (i.e, a smaller beam.) Eventually, as the field line curvature and potential become small enough a chargeseparated magnetosphere results (like the classical Goldreich-Julian picture). This is the death line predicted in open gap (OG) models. For example, Wang and Hirotani [127], estimate the death line to be approximately log 10Ṗ = −13.49 + 3.67 log 10 P. (A1) For P = 0.7 s this corresponds toĖ ≈ 1 × 10 33 erg s −1 , the lowest value we consider in Fig. 16. Other estimates of the outer gap death line are at even lower energies [e.g. 128, is closer to ≈ 2 × 10 32 erg s −1 ]. The death line evolves a fair amount depending on the amount of physics included in the calculation. The most sophisticated vacuum models are those based on 2-and 3dimensional solutions to the particle transport and Poisson equation in the gaps. The exact details of gamma-ray emission approaching this regime remain unclear, so showing the results of cutoff at variousĖ values can be used to estimate the effect of a more smooth transition to the death line.

Intrinsic Scatter
We also test whether the pulsar population can be more accurately described by including an intrinsic log-normal scatter in the total flux of the pulsars. This is motivated, in part, by the observed scatter in pulsar fluxes, as well as previous models of the radio luminosity of pulsars [e.g. 49].  In Fig. 17, we show how including the intrinsic scatter of the pulsars impacts the diffuse flux towards the GC. We include our Fiducial Model, and models with σ log L = 0.32, and 0.42. In all cases with scatter we set C = 1.0, yet the scatter effectively increases the mean of the population. For σ log L =0.32, and 0.42, the mean luminosity of the population is equal to 1.3 and 1.6 respectively, when we impose no limits on the luminosity. In practice, when L γ is selected to be above the spin down luminosity of the pulsar, we set L γ =Ė.
Overall, including the scatter only moderately impacts the observed distribution of the diffuse background. With σ log L = 0.32, we get C ≈ 1.3, and do not markedly change the distribution of point sources, except to increase the total number of bright point sources. Finally, this scatter automatically tapers the flux of pulsars withĖ ≈ 10 33 erg s −1 , because we include a hard ceiling on the flux, but impose no floor.
We found, in studying the beaming corrected fluxes of pulsars from Pierbattista et al. [9], that the intrinsic scatter is actually closer to σ log L = 0.42, which would give a mean flux  closer to C = 1.6. Using this high of a scatter requires a larger fraction of the unassociated 3FGL point sources to be pulsars that either σ log L = 0.32 or our Fiducial Model.

Beaming
In § II C, we described the beaming models for each the YRQ population and the YRL population. In Fig. 18, we show how the GC flux depends on the choice of beaming models (currently assuming C=1.0). Overall, the YRQ pulsars tend to result in less excess flux for the same model parameters (e.g., SNR, and σ log B ) as the YRL population or the fidicual model.
One interesting consequence of this empirically motivated model is that the resulting flux density distribution of the pulsar population is in better agreement with some observations. This is discussed in more detail in § VIII. This is because the observable population of YRQ pulsars has a minimum flux near L γ ≈ 10 34 erg s −1 .
One limitation of the YRQ beaming model is the num-ber of bright pulsars withĖ < ∼ 10 33.5 erg s −1 . Within this beaming model we typically expect two pulsars with a spin down luminosity below this threshold to be brighter than 10 −10 erg cm −2 s −1 , where none have been observed. Only 9 % of our simulations are consistent with the observations.

Appendix B: Regions of Interest
A number of studies focused on the emission from different regions near the GC. We show the regions of interest from Calore et al. [38] and The Fermi-LAT Collaboration [39] in Figure 19.
The green dashed line shows the flux of high latitude pulsars with |b| ≥ |l|. This region has much less contribution from the foreground disk in our YP models. However, the Fermi bubbles [129] are in the same region and are not included in our map. These have to be carefully subtracted from diffuse maps of this region. We do not include any systematic residuals that may remain.
In many analyses of the GC region, spherically symmetric templates are used to fit the diffuse emission [c.f. 37]. In order to ascertain the symmetry of the excess towards the GC region, we have divided it into two regions. The main one shown throughout the text is above and below the Galactic Plane with |b| ≥ 2 • . In Figure 20, we also show the average emission in the high latitude region with |b| ≥ |l|. Near the GC, the young pulsars from the CMZ are spherically symmetric, owing to the large kicks they receive at birth. However, outside of 5 • , the foreground YPs begin to dominate the emission. This population follows the Galactic Disk, and so the high latitude emission is reduced.