# Millisecond Pulsars from Accretion Induced Collapse naturally explain the
Galactic Center Gamma-ray Excess

Anuj Gautam1 Roland M. Crocker1 corresponding author:
<EMAIL_ADDRESS>ORCID: orcid.org/0000-0002-2036-2426 Lilia
Ferrario2 Ashley J. Ruiter3 Harrison Ploeg4 Chris Gordon4 Oscar Macias5,6
1Research School of Astronomy and Astrophysics, Australian National
University, Canberra, Australia  2Mathematical Sciences Institute, Australian
National University, Canberra, Australia  3School of Science, University of
New South Wales Canberra, Australian Defence Force Academy, Canberra,
Australia  4School of Physical and Chemical Sciences, University of
Canterbury, Christchurch, New Zealand  5Kavli Institute for the Physics and
Mathematics of the Universe (WPI), University of Tokyo, Kashiwa, Chiba
277-8583, Japan  6GRAPPA Institute, Institute of Physics, University of
Amsterdam, 1098 XH Amsterdam, The Netherlands

###### Abstract

Gamma-ray data from the Fermi-Large Area Telescope reveal an unexplained,
apparently diffuse, signal from the Galactic bulge. The origin of this
“Galactic Center Excess” (GCE) has been debated with proposed sources
prominently including self-annihilating dark matter and a hitherto undetected
population of millisecond pulsars (MSPs). We use a binary population synthesis
forward model to demonstrate that an MSP population arising from the accretion
induced collapse of O-Ne white dwarfs in Galactic bulge binaries can naturally
explain the GCE. Synchrotron emission from MSP-launched cosmic ray electrons
and positrons seems also to explain the mysterious “haze” of hard-spectrum,
non-thermal microwave emission from the inner Galaxy detected in WMAP and
Planck data.

##  I Introduction

Data collected by the Large Area Telescope (LAT) on board of NASA’s Fermi
Gamma Ray Space Telescope have revealed an excess of extended, apparently
diffuse $\gamma$-ray emission from the Galactic Center and Galactic bulge
(Goodenough and Hooper, 2009; Hooper and Goodenough, 2011; Abazajian  _et al._
, 2014; Gordon and Macías, 2013). This $\gamma$-ray feature, peaking near
${\sim}2$ GeV with an approximately spherical Daylan  _et al._ (2016)
intensity profile $\propto r^{-2.4}$ (Abazajian  _et al._ , 2014; Calore  _et
al._ , 2015), is known as the _Galactic Center Excess_ (GCE). The GCE extends
to angular radial scales of at least ${\sim}10^{\circ}$ but possibly to
${\gtrsim}20^{\circ}$ (Hooper and Slatyer, 2013; Ackermann and al. (2017),
Fermi Collaboration).

While the morphology, spectrum, and luminosity of the GCE are broadly
consistent with weakly interacting massive particle dark matter annihilation
(Goodenough and Hooper, 2009), another explanation for the signal is that it
is caused by an otherwise-undetected population of $\gamma$-ray millisecond
pulsars (MSPs) in the Galactic bulge, such a population being anticipated by
(Wang  _et al._ , 2005). Indeed, the spectrum of the GCE resembles not only
that of individual $\gamma$-ray MSPs but also that of globular clusters whose
integrated $\gamma$-ray emission was predicted, and later confirmed, to arise
from large MSP populations (Chen, 1991; Abdo and al. (2010), Fermi
Collaboration; Abazajian, 2011; Gordon and Macías, 2013).

One potential discriminant between WIMP DM and MSPs as the origin of the GCE
is that, while the former should be characterised by “true” diffuse emission,
the latter is, at least partially, expected to generate the GCE via the
combined emission from many discrete, but dim and individually-unresolvable
MSP sources. In some circumstances, statistical methods of point source
inference Collin  _et al._ (2021) may distinguish between these two types of
emission (Bartels  _et al._ , 2016; Lee  _et al._ , 2016). Efforts in this
direction using a number of different techniques are on-going but have not
been uncontested and there is on-going controversy about the reliability of
these studies in the presence of large systematic uncertainties due to
foreground model uncertainties (e.g., Leane and Slatyer, 2020a; Buschmann  _et
al._ , 2020; Chang  _et al._ , 2020).

Finally, four studies (Macias  _et al._ , 2018; Bartels  _et al._ , 2018;
Macias  _et al._ , 2019; Calore  _et al._ , 2021) have now shown that the
distribution of the GCE signal traces the stars of the inner Galaxy better
than it traces a radially symmetric distribution (as would be expected for
dark matter). This constitutes a strong argument for a stellar origin of the
GCE that is completely independent of claims and counter-claims concerning the
presence of a sub-threshold source contribution.

An argument raised against an MSP origin for the GCE is that the implied
$\gamma$-ray luminosity distribution is not consistent with the local,
resolved MSP disk population since it predicts too many bright MSPs (Hooper
and Mohlabeng, 2016). This has been disputed Ploeg  _et al._ (2020) but, in
any case, the bulge MSP population need not display the same characteristics
as that of the local Galactic disk MSP population: star formation ceased $\sim
8-10$ Gyrs ago in the bulge (Nataf, 2016; Crocker  _et al._ , 2017) while it
is ongoing in the disk. Because MSPs spin down over time and become less
$\gamma$-ray luminous – with the initially brightest MSPs evolving fastest –
there is a systematic age difference between bulge and disk MSP populations,
so that we expect a relative deficit of bright MSPs in the bulge.

To date, much of the discussion surrounding the tenability of the bulge MSP
explanation for the GCE has implicitly adopted the canonical “recycling”
scenario for the origin of the MSPs whereby an old and slowly rotating neutron
star (NS) formed in a core-collapse supernova accretes material and angular
momentum from a companion during Low-Mass X-Ray Binary (LMXB) phases
Radhakrishnan and Srinivasan (1982). Accretion spins up the NS to millisecond
periods thus enabling it to emit detectable $\gamma$ radiation (Bhattacharya
and van den Heuvel, 1991). On this basis, Cholis  _et al._ (2015) and Haggard
_et al._ (2017) raised another argument against an MSP origin of the GCE: the
Galactic bulge should host ${\sim}10^{3}$ LMXBs, while observations by the
INTEGRAL X-ray telescope have uncovered only 42 LMXBs with an additional 46
candidates within a $10^{\circ}$ radius of the Galactic Centre (Haggard  _et
al._ , 2017).

This argument, however, neglects the existence of a separate channel which may
produce $\gtrsim$ 50% of all observed MSPs (Ferrario and Wickramasinghe, 2007;
Hurley  _et al._ , 2010; Ruiter  _et al._ , 2019), namely, the accretion
induced collapse (AIC) of white dwarfs into neutron stars (NS). Here we
concentrate on AIC events that occur when ultra-massive oxygen-neon (O-Ne)
white dwarfs, accreting matter from their binary companions, approach the
Chandrasekhar mass limit (${\sim}1.4M_{\odot}$), lose electron degeneracy
pressure support due to electron capture on 24Mg and 24Na nuclei (see Jones
_et al._ , 2019, and references therein), and collapse into NSs (Miyaji  _et
al._ , 1980; Takahashi  _et al._ , 2013; Schwab  _et al._ , 2015).
Conservation of angular momentum in these events yields spin periods for the
nascent NSs (whose radii are three orders of magnitude smaller than the parent
WDs) of just a ${\sim}$few ms. Similarly, flux conservation during AIC means
that progenitor O-Ne WD magnetic fields $\sim 10^{3}-10^{5}\ \rm G$ result in
field strengths of ${\sim}10^{7}-10^{9}$ G for the newly born NSs, directly
matching observations without invoking hypothetical mechanisms for magnetic
field decay (crucial in the recycling scenario). Furthermore, because very
little mass is lost in the collapse, the binary is not disrupted, nor does the
system receive a significant natal kick (Fryer  _et al._ , 1999; Kitaura  _et
al._ , 2006). This means that NSs born via AIC have much smaller space
velocities than those imparted to NSs during the generally asymmetrical
explosion of a core collapse supernova (Lyne and Lorimer, 1994;
Wongwathanarat, A.  _et al._ , 2013; Bear and Soker, 2018). Thus, the MSPs
born via AIC are more likely to remain trapped in the gravitational potential
of the bulge allowing for a large population to build up over its history and
proffering an explanation of the match Macias  _et al._ (2018); Bartels  _et
al._ (2018); Macias  _et al._ (2019); Calore  _et al._ (2021) between the
projected bulge stellar distribution and the detailed GCE morphology.

In this paper we present the results of a binary population synthesis (BPS)
forward model for the formation and evolution of the Galactic bulge MSP
population born from the AIC of accreting O-Ne WDs to ascertain whether the
$\gamma$-ray (and X-ray) phenomenology of this population allows it to explain
the GCE.

##  II AIC Evolution scenario

O-Ne WDs are generated when stars develop carbon-oxygen (CO) cores in the
range of ${\sim}1.6-2.25\ \rm M_{\odot}$ that ignite carbon under semi-
degenerate conditions. A major aspect of binary evolution, crucial to the
production of compact binaries (Paczynski, 1976), is common envelope (CE)
evolution which occurs when a star (generally a giant) overfills its Roche
lobe engulfing both stars. A giant star with a convective envelope tends to
further expand in response to mass loss, thus triggering unstable mass
transfer at rates that are so high that the transferred matter cannot all be
captured by the companion. Instead, the matter accumulates in a common
envelope surrounding both the giant core and the companion star (de Kool,
1992). Since the envelope rotates more slowly than the stellar orbit, friction
causes the stars to spiral together and transfer orbital energy to the
envelope. This process may release sufficient energy to shed the entire
envelope, leaving either a close binary containing a WD and its companion
star, or a coalescence of the giant’s core with the companion star (see
Appendix S2).

In an AIC, while the rotational angular momentum of the collapsing O-Ne WD is
largely retained, there is some angular momentum loss to gravitational wave
(GW) radiation induced by the rapidly changing quadrupole moment
characterising the matter distribution. Although an AIC event imparts a
negligible kick (e.g. Fryer  _et al._ , 1999), some baryonic mass is suddenly
lost, causing a drop in the effective gravitational mass and an expansion of
the binary. Depending on system parameters, the stars may come into contact
again and mass transfer resume. We find that donor stars can transfer more
than $0.1M_{\odot}$ of material to the NS and thus post-AIC accretion is a
very important determinant of the spin evolution of the AIC MSP population.
Given the $\gamma$-ray luminosity of a NS is powered by its rotational kinetic
energy, our modelling includes all accretion torques acting on the NS
following its formation (see Appendix S4).

##  III Generation of reference data set

To evolve a model population of binaries up to AIC, we have used the state-of-
the-art binary-star evolution (BSE) code Hurley  _et al._ (2000, 2002, 2013).
BSE takes as input stellar zero age main sequence (ZAMS) masses, $M_{\rm
1,ZAMS}$ and $M_{\rm 2,ZAMS}$ for the primary and secondary respectively
(where $M_{\rm 1,ZAMS}$ $\geq$ $M_{\rm 2,ZAMS}$), and the initial binary
separation, $a$. We adopt the following empirically-motivated, “off-the-shelf”
parameter choices (see Appendix S1.1 for detailed justifications and
motivations): i) We assume a binary star fraction of 70% Ruiter  _et al._
(2019); ii) We draw $M_{1,ZAMS}$ from the initial mass function of Kroupa  _et
al._ (1993) iii) Secondary masses are drawn from a flat mass ratio
distribution Ferrario (2012); Hurley  _et al._ (2002); iv) Initial
eccentricities are set to zero; v) Orbital periods at zero age are drawn from
a log-uniform distribution and range from $10$ days to $10^{4}$ days Duchêne
and Kraus (2013); vi) The common envelope efficiency factor (see Appendix S2)
is set to $\alpha_{\rm CE}=1$ Ruiter  _et al._ (2019); vii) The binding energy
parameter $\lambda_{\rm b}$ is calculated using the Cambridge STARS code
(Eggleton, 1971; Pols  _et al._ , 1995, 1998); and viii) Magnetic fields are
assigned to each nascent AIC NS via random sampling from the log-normal
magnetic field distribution of Ploeg  _et al._ (2020). Our synthetic
population consists of ${\sim}6.3\times 10^{7}$ binary systems which, given
the parameters above, would be hosted by an overall stellar population of
$2.0\times 10^{9}M_{\odot}$ ZAMS mass.

We numerically evolve the binaries up to 14 Gyrs generating a total of
${\sim}9.1\times 10^{3}$ AIC events, which start when the population is only
${\sim}0.2$ Gyr old (Figure 1). In the present study we have neglected Merger-
Induced Collapse (MIC) events that occur when two massive WDs merge due to
gravitational wave radiation (Ivanova  _et al._ , 2008). This is because they
occur relatively early post star-formation (Ruiter  _et al._ , 2019; Lyutikov
and Toonen, 2019), supply less than 10% of the MSP population (Ruiter  _et
al._ , 2019), and, in any case, produce nascent NSs with super-strong magnetic
fields (e.g. Wickramasinghe  _et al._ , 2014) that quickly spin down and no
longer qualify as MSPs today.

##  IV Application to the Galactic bulge

In order to apply our model results to the Galactic bulge, we need to scale up
from the $2.0\times 10^{9}\ M_{\odot}$ ZAMS equivalent stellar mass (whose
binary population we have modeled) to the bulge ZAMS mass. The latter quantity
can be inferred from its current measured (Portail  _et al._ , 2015) stellar
mass $(1.55\pm 0.15)\times 10^{10}\rm M_{\odot}$ (including stellar remnants),
its star formation history (with the overwhelming majority of stars formed
more than 8 Gyr ago Nataf (2016); Crocker  _et al._ (2017)), and a model for
the mass loss by stars in winds and ejecta over the course of their evolution
(Maraston, 1998). The current bulge stellar mass needs to scaled upwards by
$1.81\pm 0.02$ to obtain the ZAMS mass 111This error range purely reflects the
uncertainty arising from the star formation history., leading to an overall
rescaling from the ZAMS mass implied by the binary population we have modeled
to the bulge of $14.0\pm 4.1$. This implies $(1.3\pm 0.4)\times 10^{5}$ AIC
events have occurred in the bulge up to the present (cf. Figure 1) and, in the
present epoch, $(1.1\pm 0.3)\times 10^{5}$ MSPs in the bulge, and an AIC
expected every few Myr.

Figure 1: Number of model bulge MSPs (periods as labeled) and cumulative AIC
events over cosmological time (note that some AIC NSs spin down to $P>40$ ms
in which case they cease being classed as MSPs by us). The $\pm 1\sigma$ error
band on the red curve for all AIC events reflects the uncertainties stemming
from the bulge stellar mass determination and in the bulge binarity fraction.
(For clarity, equivalent error bands on the other curves are not shown.)

###  IV.1 Predicted Bulge MSP $\gamma$-ray emission

Empirically, ${\sim}10\%$ of the kinetic energy of rotation of an MSP is lost
to dipole radiation Abdo and al. (2013) (Fermi Collaboration) which emerges at
$\gamma$-ray wavelengths due to “prompt” emission by $e^{+}-e^{-}$ pairs
within the magnetosphere. A further, poorly-determined fraction (1 – 90%:
(Sudoh  _et al._ , 2020; Song  _et al._ , 2021)), is carried away in
relativistic pairs that escape the magnetosphere and, after possible re-
acceleration at pulsar wind termination or intra-binary shocks, generates
additional synchrotron and $\gamma$-ray inverse Compton (IC) signals on
interstellar medium (ISM) fields. Ultimately, both prompt and delayed signals
are powered by the liberation, via magnetic dipole braking, of the rotational
kinetic energy of the MSP population. In Figure 2 we show the evolution of
this spin down power – which thus sets an upper limit on the total (i.e.,
prompt + IC) $\gamma$-ray luminosity of the model MSP population – together
with a datum indicating the measured Macias  _et al._ (2018) GCE luminosity.
Evidently, there is sufficient power available to explain the signal at 13.8
Gyr cosmological (current) time.

Figure 2: Spin-down power liberated by magnetic dipole braking over
cosmological time (band indicates $\pm 1\sigma$ error). The time series has
been smoothed over 0.5 Gyr. The luminosity datum at 13.8 Gyr, $2.2\pm
0.4\times 10^{37}$ erg/s, is inferred from the Boxy Bulge spectrum measured by
Macias  _et al._ (2018) together with a weighted mean distance (as obtained
from Monte Carlo sampling of the spatial model (Ploeg  _et al._ , 2020)) to
each MSP of $\langle d_{\rm MSP}\rangle=7.89$ kpc).

We next wish to determine the spectrum at Earth of the $\gamma$-ray emission
due to our model bulge MSP population. Given that our BPS gives us the
instantaneous period and period derivative of each MSP, we use an empirical
prompt $\gamma$-ray luminosity model of the form Gonthier  _et al._ (2018):

$L_{\rm\gamma,prompt,i}=f_{\rm\gamma,prompt}\left(P_{i}/{\rm
s}\right)^{\alpha_{\gamma}}\dot{P}_{\rm MB,i}^{\beta_{\gamma}}\ \rm erg\ \
s^{-1},$ (1)

where $P_{i}$ is the current period of the $i$th MSP and $\dot{P}_{\rm MB,i}$
is the time derivative of its period due to magnetic dipole braking. This
luminosity is distributed over a spectrum (Ploeg  _et al._ , 2020):

$\frac{dN_{\rm\gamma,prompt,i}}{dE_{\gamma}}\propto
E_{\gamma}^{-\Gamma_{i}}e^{(-E_{\gamma}/E_{\rm cut,i})^{2/3}}\,.$ (2)

Then the prompt flux received from each model MSP is determined from

$F_{\rm\gamma,prompt,i}\equiv\frac{L_{\rm\gamma,prompt,i}}{4\pi d_{i}^{2}}\,.$
(3)

Finally, after assigning such distances and the spectral and other parameters
entering into the equations above – which we do by sampling from the spatial
distribution determined in Ploeg  _et al._ (2020); see Appendix S6 – we can
generate an initial estimate of the prompt spectrum detected from the entire
bulge MSP population. This is shown as the dot-dash blue curve in Figure 3.

Figure 3: Fitted $\gamma$-ray spectrum at Earth of the model MSP population at
13.8 Gyr together with the measured GCE spectral data points Macias  _et al._
(2019). The dot-dash blue line shows the central value of the spectrum
directly-predicted by our BPS model; all solid lines have been fitted
following the procedure outlined in Appendix S6 (bands indicate $\pm 1\sigma$
error). The spectral index of the IC-emitting CR $e^{\pm}$ population is
$\gamma_{e}=-2.1^{+0.1}_{-0.2}$ and its cutoff energy is $(9^{+10}_{-2})\times
10^{11}$ eV.

From Figure 3 it is evident that the predicted prompt spectrum at Earth is
promisingly similar to that measured and, in particular, does a good job of
reproducing the ${\sim}2$ GeV bump). On the other hand, i) the amplitude of
the central value of the predicted signal is somewhat low and ii) the spectral
shape is not a good fit to the tail at high energies ($E_{\gamma}\gtrsim 10$
GeV)222The overall $\chi^{2}$ obtained from the fit of the predicted prompt
spectrum to the GCE data points is 346 for 15 data points (and no tunable
parameters), accounting for the errors in both the predicted prompt spectrum
and the spectral measurements.. However, there is a potential further source
of $\gamma$-ray emission powered by the MSP spin down that could contribute to
the observed GCE signal, in particular, the high-energy tail: IC off photons
of the interstellar radiation field by $e^{\pm}$ that escape MSP
magnetospheres (e.g. Yuan and Ioka (2015); Petrović  _et al._ (2015); Calore
_et al._ (2015); Horiuchi  _et al._ (2016); Linden  _et al._ (2016); Song  _et
al._ (2019, 2021); Macias  _et al._ (2021)). This process can supply
spectrally hard emission up to $\gtrsim 30$ GeV, as required (e.g., Song  _et
al._ (2021)). Moreover, our $\gamma$-ray template analysis reveals that the
GCE signal is smeared out over a vertical angular scale of $\sim 5^{\circ}$,
an observation that strongly suggests that IC emission significantly
contributes to it (see Appendix S7.2).

Now, while the total prompt emission is prescribed (up to statistical
fluctuations) once we have the $\\{P_{i},\dot{P}_{\rm MB,i}\\}$ from our model
MSP population, this is not the case for the delayed IC signal which is
subject to a number of factors not determined within out model (see Appendix
S7). On the other hand, the IC and IC + prompt signals are subject to two
significant constraints: i) synchrotron radiation off the bulge magnetic field
by the population of CR $e^{\pm}$ escaping MSP magnetospheres (and also
supplying the IC signal) should not exceed the measured, microwave frequency
flux density emitted from the inner Galaxy Planck Collaboration (2013) and ii)
the overall prompt + IC + synchrotron luminosity may not exceed the total
power liberated by magnetic braking of the MSP population.

Using these constraints, we perform a simultaneous fitting of the IC and
prompt signals to the GCE data points allowing the amplitude of the prompt
signal to float within its uncertainty range; the fitted spectra and their
uncertainty bands are shown in Figure 3. Evidently, an old population of MSPs
can reproduce very well the measured $\gamma$-ray emission from the GCE.
Turning to Figure 4 which shows the broad-band, non-thermal emission from our
fit, note that we also saturate (by construction) the microwave data points
Planck Collaboration (2013) that describe the spectrum of the “microwave
haze”, an enduringly mysterious, non-thermal (total) emission feature of the
inner Galaxy first discovered in WMAP data Finkbeiner (2004) and for which a
number of previous works have considered a potential origin in inner Galaxy
pulsars or MSPs Zhang  _et al._ (2009); Kaplinghat  _et al._ (2009); Malyshev
_et al._ (2010); Harding and Abazajian (2010); McQuinn and Zaldarriaga (2011).
Also note that our fit here requires that a significant fraction, close to
100%, of the spin-down power liberated by magnetic braking ends up in
relativistic pairs (cf.Sudoh  _et al._ (2020)).

Figure 4: Fitted $\gamma$-ray (colors as labeled in Figure 3) and
radio/microwave synchrotron (orange) spectra at Earth of the model MSP
population at 13.8 Gyr together with the measured GCE Macias  _et al._ (2019)
and microwave haze Planck Collaboration (2013) spectral data points. The color
band indicates the $\pm 1\sigma$ error on the fitted synchrotron spectrum
(error bands on the $\gamma$-ray curves are suppressed for clarity, but are as
shown in Figure 3); all solid lines have been fitted following the procedure
outlined in Appendix S6. Best fit parameters are as given in the caption of
Figure 3 plus a magnetic field $(2.3^{+0.6}_{-0.4})\times 10^{-5}\ \mu$G.

###  IV.2 Predicted population of bright, $\gamma-$ray resolvable Bulge MSPs

In Figure 5 (top) we show the absolute (energy) flux $F_{E}$ distribution of
our model bulge AIC MSP population at 13.8 Gyr. The middle panel shows a zoom-
in, close to the detection threshold (corresponding to 0 on the x-axis), of
the histogrammed distribution of the ratios of each MSP’s prompt flux to the
flux required for likely detection, either as a $\gamma-$ray point source
($F_{\rm sens,point}$) or (more stringently) as an individual MSP ($F_{\rm
sens,MSP}$), at the MSP’s modelled position; see Appendix S10. Clearly very
few MSPs are likely to be individually detected. Note that most of the IC-
emitting CR $e^{\pm}$ launched by each MSP will be transported outside the
aperture 333The Fermi-LAT has a 68% containment radius of $\sim 0.1^{\circ}$
at $\sim$ 10 GeV Atwood and al. (2009) (Fermi Collaboration) which corresponds
to a lineal size of $\sim$ 14 pc in the bulge. defined by the Fermi-LAT point
spread function before IC radiating (cf. Figure S19), so the IC flux
contributed by each MSP does not contribute to its point source flux.

Figure 5: Top: Total and prompt energy flux distributions of the modeled MSP
population at $t=13.8$ Gyrs cosmological time. The vertical dashed lines show
the mean and $\pm 1\sigma$ point-source detection thresholds Ploeg  _et al._
(2020). Center: zoom-in (close to the flux detection thresholds) of the MSP
population prompt energy fluxes normalised to the Fermi-LAT detection
threshold for i) MSPs and ii) $\gamma$-ray point sources Bottom: the number
flux distribution of the model bulge MSP population over the solid angle
defined by $|b|>2^{\circ}$ and $r<25^{\circ}$; the vertical, dashed, blue line
is the number flux equivalent to an expected 1 photon per source; the
vertical, dashed, gray lines subtend the (variable) approximate point source
flux sensitivity threshold Chang  _et al._ (2020). Note that the population
contains no MSPs with $F_{E}\gtrsim 10^{-11}$ erg/s/cm2 or $F\gtrsim 10^{-10}$
ph/s/cm2

Considering only prompt emission, we determine a mean number of Bulge MSPs
expected to be detected by Fermi of $0.15_{0.0}^{0.37}$444Including IC
emission into the point source flux – certainly an overestimate – we determine
$7.3\pm 3.3$ MSPs are expected to be resolvable. Ploeg  _et al._ (2020)
identified three resolved MSPs that have non-negligible probabilities of
belonging to the bulge.. The photon number flux distribution of the population
is shown in the bottom panel of Figure 5. We discuss our predictions for the
flux distribution of the putative MSP population in the context of results
obtained using point source inference techniques in Appendix S11.

##  V Conclusions

We have evolved a model population of Galactic bulge binaries to the present
day; these generate ${\sim}1.3\times 10^{5}$ accretion induced collapse (AIC)
neutron stars whose present-day spin periods ($<40$ ms) qualify them as MSPs.
Individually, these MSPs cannot be distinguished from those that form via
“recycling” of old NSs born from core collapse supernovae. At a population
level, bulge MSPs are older and thus dimmer than typical local MSPs belonging
to the Galactic disk population and characterised by a luminosity distribution
$10^{31\pm 3}$ erg/s including both prompt and IC components (and taking the
mean and standard deviation in the log). The model bulge AIC MSP population
can fully explain both the spectrum and overall luminosity of the GCE and, as
a bonus, the microwave haze from the inner Galaxy. We find that the synthetic
MSP population also reproduces the correct spectral shape of the GCE yet does
not predict too many bright and resolvable $\gamma$-ray MSPs; we show in
Appendix S12 that, in addition, our model does not predict too many bright,
individual X-ray sources nor an excessive integrated X-ray luminosity.

The upcoming Cherenkov Telescope Array (CTA) will look at the sky in photons
with energies ranging from $20\ \rm GeV$ to $300\ \rm TeV$ and is expected to
spend at least 500 hours observing the inner Galaxy Cherenkov Telescope Array
Consortium (2019) in which case we expect its point source sensitivity at 0.1
TeV to approach few $\times 10^{-13}$ erg cm-2 s-1 Maier (2019). Thus, CTA
should resolve at least a handful of GCE MSPs according to the best-fit
spectral parameters reported above (cf. Macias  _et al._ (2021)). In addition,
we expect point source inference techniques directed at inferring the flux
distribution of the sub-threshold population in Fermi-LAT data (e.g. List  _et
al._ , 2020) will, as they reach maturity and become robust to systematics,
start to test our MSP population predictions.

## Acknowledgements

RMC thanks Mark Krumholz, Ivo Seitenzahl, Gavin Rowell, Celine Boehm, and
Matthew Baring for enlightening conversations. Nick Rodd and Shunsaku Horiuchi
are thanked for a close reading of the paper in draft. RMC acknowledges
support from the Australian Government through the Australian Research
Council, award DP190101258 (shared with Prof. Mark Krumholz). AJR has been
supported by the Australian Research Council through grant number FT170100243.
Parts of this research were undertaken with the assistance of resources and
services from the National Computational Infrastructure (NCI), which is
supported by the Australian Government, through the National Computational
Merit Allocation Scheme and the UNSW HPC Resource Allocation Scheme. The work
of O.M. was supported by JSPS KAKENHI Grant Number JP20K14463 and by the World
Premier International Research Center Initiative (WPI Initiative), MEXT,
Japan.

## References

  * Goodenough and Hooper (2009) L. Goodenough and D. Hooper, Possible evidence for dark matter annihilation in the inner milky way from the fermi gamma ray space telescope (2009), arXiv:0910.2998 [hep-ph] .
  * Hooper and Goodenough (2011) D. Hooper and L. Goodenough, Phys. Lett. B 697, 412 (2011).
  * Abazajian  _et al._ (2014) K. N. Abazajian, N. Canac, S. Horiuchi, and M. Kaplinghat, Physical Review D 90, 10.1103/physrevd.90.023526 (2014).
  * Gordon and Macías (2013) C. Gordon and O. Macías, Physical Review D 88, 10.1103/physrevd.88.083521 (2013).
  * Daylan  _et al._ (2016) T. Daylan, D. P. Finkbeiner, D. Hooper, T. Linden, S. K. N. Portillo, N. L. Rodd, and T. R. Slatyer, Physics of the Dark Universe 12, 1 (2016), arXiv:1402.6703 [astro-ph.HE] .
  * Calore  _et al._ (2015) F. Calore, I. Cholis, C. McCabe, and C. Weniger, Physical Review D 91, 10.1103/physrevd.91.063003 (2015).
  * Hooper and Slatyer (2013) D. Hooper and T. R. Slatyer, Physics of the Dark Universe 2, 118 (2013), arXiv:1302.6589 [astro-ph.HE] .
  * Ackermann and al. (2017) (Fermi Collaboration) M. Ackermann and al. (Fermi Collaboration), Astrophys. J. 840, 43 (2017).
  * Wang  _et al._ (2005) W. Wang, Z. J. Jiang, and K. S. Cheng, Mon. Not. Roy. Astron. Soc. 358, 263 (2005), arXiv:astro-ph/0501245 [astro-ph] .
  * Chen (1991) K. Chen, Nature (London) 352, 695 (1991).
  * Abdo and al. (2010) (Fermi Collaboration) A. A. Abdo and al. (Fermi Collaboration), Astron. and Astrophys. 524, A75 (2010), arXiv:1003.3588 [astro-ph.GA] .
  * Abazajian (2011) K. N. Abazajian, J. Cosmol. Astropart. Phys. 2011, 010–010 (2011).
  * Collin  _et al._ (2021) G. H. Collin, N. L. Rodd, T. Erjavec, and K. Perez, arXiv e-prints , arXiv:2104.04529 (2021), arXiv:2104.04529 [astro-ph.IM] .
  * Bartels  _et al._ (2016) R. Bartels, S. Krishnamurthy, and C. Weniger, Phys. Rev. Lett. 116, 10.1103/physrevlett.116.051102 (2016).
  * Lee  _et al._ (2016) S. K. Lee, M. Lisanti, B. R. Safdi, T. R. Slatyer, and W. Xue, Phys. Rev. Lett. 116, 051103 (2016).
  * Leane and Slatyer (2020a) R. K. Leane and T. R. Slatyer, Phys. Rev. Lett. 125, 121105 (2020a), arXiv:2002.12370 [astro-ph.HE] .
  * Buschmann  _et al._ (2020) M. Buschmann, N. L. Rodd, B. R. Safdi, L. J. Chang, S. Mishra-Sharma, M. Lisanti, and O. Macias, Phys. Rev. D 102, 023023 (2020), arXiv:2002.12373 [astro-ph.HE] .
  * Chang  _et al._ (2020) L. J. Chang, S. Mishra-Sharma, M. Lisanti, M. Buschmann, N. L. Rodd, and B. R. Safdi, Phys. Rev. D 101, 023014 (2020), arXiv:1908.10874 [astro-ph.CO] .
  * Macias  _et al._ (2018) O. Macias, C. Gordon, R. M. Crocker, B. Coleman, D. Paterson, S. Horiuchi, and M. Pohl, Nature Astronomy 2, 387 (2018).
  * Bartels  _et al._ (2018) R. Bartels, E. Storm, C. Weniger, and F. Calore, Nature Astronomy 2, 819–828 (2018).
  * Macias  _et al._ (2019) O. Macias, S. Horiuchi, M. Kaplinghat, C. Gordon, R. M. Crocker, and D. M. Nataf, J. Cosmol. Astropart. Phys. 2019, 042–042 (2019).
  * Calore  _et al._ (2021) F. Calore, F. Donato, and S. Manconi, arXiv e-prints , arXiv:2102.12497 (2021), arXiv:2102.12497 [astro-ph.HE] .
  * Hooper and Mohlabeng (2016) D. Hooper and G. Mohlabeng, J. Cosmol. Astropart. Phys. 2016, 049–049 (2016).
  * Ploeg  _et al._ (2020) H. Ploeg, C. Gordon, R. Crocker, and O. Macias, J. Cosmo. Astropart. Phys. 2020, 035 (2020), arXiv:2008.10821 [astro-ph.HE] .
  * Nataf (2016) D. M. Nataf, Publications of the Astronomical Society of Australia 33, e023 (2016).
  * Crocker  _et al._ (2017) R. M. Crocker, A. J. Ruiter, I. R. Seitenzahl, F. H. Panther, S. Sim, H. Baumgardt, A. Möller, D. M. Nataf, L. Ferrario, J. J. Eldridge, M. White, B. E. Tucker, and F. Aharonian, Nature Astronomy 1, 0135 (2017).
  * Radhakrishnan and Srinivasan (1982) V. Radhakrishnan and G. Srinivasan, Current Science 51, 1096 (1982).
  * Bhattacharya and van den Heuvel (1991) D. Bhattacharya and E. P. J. van den Heuvel, Physics Reports 203, 1 (1991).
  * Cholis  _et al._ (2015) I. Cholis, D. Hooper, and T. Linden, Physical Review D 91, 10.1103/physrevd.91.083507 (2015).
  * Haggard  _et al._ (2017) D. Haggard, C. Heinke, D. Hooper, and T. Linden, J. Cosmol. Astropart. Phys. 2017, 056–056 (2017).
  * Ferrario and Wickramasinghe (2007) L. Ferrario and D. Wickramasinghe, Mon. Not. Roy. Astron. Soc. 375, 1009 (2007).
  * Hurley  _et al._ (2010) J. R. Hurley, C. A. Tout, D. T. Wickramasinghe, L. Ferrario, and P. D. Kiel, Mon. Not. Roy. Astron. Soc. 402, 1437 (2010).
  * Ruiter  _et al._ (2019) A. J. Ruiter, L. Ferrario, K. Belczynski, I. R. Seitenzahl, R. M. Crocker, and A. I. Karakas, Mon. Not. Roy. Astron. Soc. 484, 698 (2019), arXiv:1802.02437 [astro-ph.SR] .
  * Jones  _et al._ (2019) S. Jones, F. K. Röpke, C. Fryer, A. J. Ruiter, I. R. Seitenzahl, L. R. Nittler, S. T. Ohlmann, R. Reifarth, M. Pignatari, and K. Belczynski, Astron. and Astrophys. 622, A74 (2019), arXiv:1812.08230 [astro-ph.SR] .
  * Miyaji  _et al._ (1980) S. Miyaji, K. Nomoto, K. Yokoi, and D. Sugimoto, Publications of the Astronomical Society of Japan 32, 303 (1980).
  * Takahashi  _et al._ (2013) K. Takahashi, T. Yoshida, and H. Umeda, Astrophys. J. 771, 28 (2013), arXiv:1302.6402 [astro-ph.SR] .
  * Schwab  _et al._ (2015) J. Schwab, E. Quataert, and L. Bildsten, Mon. Not. Roy. Astron. Soc. 453, 1910 (2015), arXiv:1504.05194 [astro-ph.SR] .
  * Fryer  _et al._ (1999) C. Fryer, W. Benz, M. Herant, and S. A. Colgate, Astrophys. J. 516, 892 (1999).
  * Kitaura  _et al._ (2006) F. S. Kitaura, H.-T. Janka, and W. Hillebrandt, Astronomy & Astrophysics 450, 345–350 (2006).
  * Lyne and Lorimer (1994) A. G. Lyne and D. R. Lorimer, Nature 369, 127 (1994).
  * Wongwathanarat, A.  _et al._ (2013) Wongwathanarat, A., Janka, H.-Th., and Müller, E., Astron. and Astrophys. 552, A126 (2013).
  * Bear and Soker (2018) E. Bear and N. Soker, Astrophys. J. 855, 82 (2018).
  * Paczynski (1976) B. Paczynski, in  _Structure and Evolution of Close Binary Systems_ , IAU Symposium, Vol. 73, edited by P. Eggleton, S. Mitton, and J. Whelan (1976) p. 75.
  * de Kool (1992) M. de Kool, aap 261, 188 (1992).
  * Hurley  _et al._ (2000) J. R. Hurley, O. R. Pols, and C. A. Tout, Mon. Not. Roy. Astron. Soc. 315, 543 (2000), arXiv:astro-ph/0001295 [astro-ph] .
  * Hurley  _et al._ (2002) J. R. Hurley, C. A. Tout, and O. R. Pols, Mon. Not. Roy. Astron. Soc. 329, 897 (2002), arXiv:astro-ph/0201220 [astro-ph] .
  * Hurley  _et al._ (2013) J. R. Hurley, C. A. Tout, and O. R. Pols, BSE: Binary Star Evolution (2013), ascl:1303.014 .
  * Kroupa  _et al._ (1993) P. Kroupa, C. A. Tout, and G. Gilmore, Mon. Not. Roy. Astron. Soc. 262, 545 (1993).
  * Ferrario (2012) L. Ferrario, Mon. Not. Roy. Astron. Soc. 426, 2500 (2012), arXiv:1209.1427 [astro-ph.SR] .
  * Duchêne and Kraus (2013) G. Duchêne and A. Kraus, Annu. Rev. Astron Astrophys. 51, 269 (2013), arXiv:1303.3028 [astro-ph.SR] .
  * Eggleton (1971) P. P. Eggleton, Mon. Not. Roy. Astron. Soc. 151, 351 (1971).
  * Pols  _et al._ (1995) O. R. Pols, C. A. Tout, P. P. Eggleton, and Z. Han, Mon. Not. Roy. Astron. Soc. 274, 964 (1995), arXiv:astro-ph/9504025 [astro-ph] .
  * Pols  _et al._ (1998) O. R. Pols, K.-P. Schröder, J. R. Hurley, C. A. Tout, and P. P. Eggleton, Mon. Not. Roy. Astron. Soc. 298, 525 (1998).
  * Ivanova  _et al._ (2008) N. Ivanova, C. O. Heinke, F. A. Rasio, K. Belczynski, and J. M. Fregeau, Mon. Not. Roy. Astron. Soc. 386, 553 (2008).
  * Lyutikov and Toonen (2019) M. Lyutikov and S. Toonen, Mon. Not. Roy. Astron. Soc. 487, 5618 (2019), arXiv:1812.07569 [astro-ph.HE] .
  * Wickramasinghe  _et al._ (2014) D. T. Wickramasinghe, C. A. Tout, and L. Ferrario, Mon. Not. Roy. Astron. Soc. 437, 675 (2014), arXiv:1310.2696 [astro-ph.SR] .
  * Portail  _et al._ (2015) M. Portail, C. Wegg, O. Gerhard, and I. Martinez-Valpuesta, Mon. Not. Roy. Astron. Soc. 448, 713 (2015), arXiv:1502.00633 [astro-ph.GA] .
  * Maraston (1998) C. Maraston, Mon. Not. Roy. Astron. Soc. 300, 872 (1998), arXiv:astro-ph/9807338 [astro-ph] .
  * Note (1) This error range purely reflects the uncertainty arising from the star formation history.
  * Abdo and al. (2013) (Fermi Collaboration) A. A. Abdo and al. (Fermi Collaboration), "Astrophys. J. Suppl." 208, 17 (2013), arXiv:1305.4385 [astro-ph.HE] .
  * Sudoh  _et al._ (2020) T. Sudoh, T. Linden, and J. F. Beacom, arXiv e-prints , arXiv:2005.08982 (2020), arXiv:2005.08982 [astro-ph.GA] .
  * Song  _et al._ (2021) D. Song, O. Macias, S. Horiuchi, R. M. Crocker, and D. M. Nataf, arXiv e-prints , arXiv:2102.00061 (2021), arXiv:2102.00061 [astro-ph.HE] .
  * Gonthier  _et al._ (2018) P. L. Gonthier, A. K. Harding, E. C. Ferrara, S. E. Frederick, V. E. Mohr, and Y.-M. Koh, Astrophys. J. 863, 199 (2018).
  * Note (2) The overall $\chi^{2}$ obtained from the fit of the predicted prompt spectrum to the GCE data points is 346 for 15 data points (and no tunable parameters), accounting for the errors in both the predicted prompt spectrum and the spectral measurements.
  * Yuan and Ioka (2015) Q. Yuan and K. Ioka, Astrophys. J. 802, 124 (2015), arXiv:1411.4363 [astro-ph.HE] .
  * Petrović  _et al._ (2015) J. Petrović, P. D. Serpico, and G. Zaharijas, J. Cosmo. Astropart. Phys. 2015, 023 (2015), arXiv:1411.2980 [astro-ph.HE] .
  * Horiuchi  _et al._ (2016) S. Horiuchi, M. Kaplinghat, and A. Kwa, J. Cosmo. Astropart. Phys. 2016, 053 (2016), arXiv:1604.01402 [astro-ph.HE] .
  * Linden  _et al._ (2016) T. Linden, N. L. Rodd, B. R. Safdi, and T. R. Slatyer, Phys. Rev. D 94, 103013 (2016), arXiv:1604.01026 [astro-ph.HE] .
  * Song  _et al._ (2019) D. Song, O. Macias, and S. Horiuchi, Phys. Rev. D 99, 123020 (2019), arXiv:1901.07025 [astro-ph.HE] .
  * Macias  _et al._ (2021) O. Macias, H. van Leijen, D. Song, S. Ando, S. Horiuchi, and R. M. Crocker, (2021), arXiv:2102.05648 [astro-ph.HE] .
  * Planck Collaboration (2013) Planck Collaboration, Astron. and Astrophys. 554, A139 (2013), arXiv:1208.5483 [astro-ph.GA] .
  * Finkbeiner (2004) D. P. Finkbeiner, Astrophys. J. 614, 186 (2004), arXiv:astro-ph/0311547 [astro-ph] .
  * Zhang  _et al._ (2009) J. Zhang, X.-J. Bi, J. Liu, S.-M. Liu, P.-F. Yin, Q. Yuan, and S.-H. Zhu, Phys. Rev. D 80, 023007 (2009), arXiv:0812.0522 [astro-ph] .
  * Kaplinghat  _et al._ (2009) M. Kaplinghat, D. J. Phalen, and K. M. Zurek, J. Cosmo. Astropart. Phys. 2009, 010 (2009), arXiv:0905.0487 [astro-ph.HE] .
  * Malyshev  _et al._ (2010) D. Malyshev, I. Cholis, and J. D. Gelfand, Astrophys. J. 722, 1939 (2010), arXiv:1002.0587 [astro-ph.HE] .
  * Harding and Abazajian (2010) J. P. Harding and K. N. Abazajian, Phys. Rev. D 81, 023505 (2010), arXiv:0910.4590 [astro-ph.CO] .
  * McQuinn and Zaldarriaga (2011) M. McQuinn and M. Zaldarriaga, Mon. Not. Roy. Astron. Soc. 414, 3577 (2011), arXiv:1004.1189 [astro-ph.HE] .
  * Note (3) The Fermi-LAT has a 68% containment radius of $\sim 0.1^{\circ}$ at $\sim$ 10 GeV Atwood and al. (2009) (Fermi Collaboration) which corresponds to a lineal size of $\sim$ 14 pc in the bulge.
  * Note (4) Including IC emission into the point source flux – certainly an overestimate – we determine $7.3\pm 3.3$ MSPs are expected to be resolvable. Ploeg  _et al._ (2020) identified three resolved MSPs that have non-negligible probabilities of belonging to the bulge.
  * Cherenkov Telescope Array Consortium (2019) Cherenkov Telescope Array Consortium,  _Science with the Cherenkov Telescope Array_ (2019).
  * Maier (2019) G. Maier, in  _36th International Cosmic Ray Conference (ICRC2019)_ , International Cosmic Ray Conference, Vol. 36 (2019) p. 733, arXiv:1907.08171 [astro-ph.IM] .
  * List  _et al._ (2020) F. List, N. L. Rodd, G. F. Lewis, and I. Bhat, Phys. Rev. Lett. 125, 241102 (2020), arXiv:2006.12504 [astro-ph.HE] .
  * Atwood and al. (2009) (Fermi Collaboration) W. B. Atwood and al. (Fermi Collaboration), Astrophys. J. 697, 1071 (2009), arXiv:0902.1089 [astro-ph.IM] .
  * Ploeg and Gordon (2021) H. Ploeg and C. Gordon, arXiv e-prints , arXiv:2105.13034 (2021), arXiv:2105.13034 [astro-ph.HE] .
  * Kiel and Hurley (2006) P. D. Kiel and J. R. Hurley, Mon. Not. Roy. Astron. Soc. 369, 1152 (2006).
  * Briggs  _et al._ (2015) G. P. Briggs, L. Ferrario, C. A. Tout, D. T. Wickramasinghe, and J. R. Hurley, Mon. Not. Roy. Astron. Soc. 447, 1713 (2015), arXiv:1412.5662 [astro-ph.SR] .
  * Moe and Stefano (2013) M. Moe and R. D. Stefano, Astrophys. J. 778, 95 (2013).
  * Mazeh  _et al._ (2003) T. Mazeh, M. Simon, L. Prato, B. Markus, and S. Zucker, Astrophys. J. 599, 1344 (2003).
  * Willems and Kolb (2002) B. Willems and U. Kolb, Mon. Not. Roy. Astron. Soc. 337, 1004 (2002), https://onlinelibrary.wiley.com/doi/pdf/10.1046/j.1365-8711.2002.05985.x .
  * Moe and Di Stefano (2017) M. Moe and R. Di Stefano, "Astrophys. J. Suppl." 230, 15 (2017), arXiv:1606.05347 [astro-ph.SR] .
  * Tout  _et al._ (1997) C. A. Tout, S. J. Aarseth, O. R. Pols, and P. P. Eggleton, Mon. Not. Roy. Astron. Soc. 291, 732 (1997).
  * Webbink (1984) R. F. Webbink, Astrophys. J. 277, 355 (1984).
  * Pfahl  _et al._ (2003) E. Pfahl, S. Rappaport, and P. Podsiadlowski, Astrophys. J. 597, 1036 (2003).
  * Tutukov and Yungelson (1996) A. Tutukov and L. Yungelson, Mon. Not. Roy. Astron. Soc. 280, 1035 (1996).
  * Watts and Andersson (2002) A. L. Watts and N. Andersson, Mon. Not. Roy. Astron. Soc. 333, 943 (2002), arXiv:astro-ph/0110573 [astro-ph] .
  * Lindblom  _et al._ (1998) L. Lindblom, B. J. Owen, and S. M. Morsink, Phys. Rev. Lett. 80, 4843–4846 (1998).
  * Possenti  _et al._ (2000) A. Possenti, M. Colpi, U. Geppert, L. Burderi, and N. D’Amico, Recycling NSs to ultrashort periods: a statistical analysis of their evolution in the mu-P plane, in  _IAU Colloq. 177: Pulsar Astronomy - 2000 and Beyond_ , Astronomical Society of the Pacific Conference Series, Vol. 202, edited by M. Kramer, N. Wex, and R. Wielebinski (2000) p. 585.
  * Gittins and Andersson (2019) F. Gittins and N. Andersson, Mon. Not. Roy. Astron. Soc. 488, 99–110 (2019).
  * Wang (1996) Y. M. Wang, Astrophys. J. Lett. 465, L111 (1996).
  * Lipunov (1992) V. M. Lipunov,  _Astrophysics of Neutron Stars_ (1992).
  * Boldin and Popov (2010) P. A. Boldin and S. B. Popov, Mon. Not. Roy. Astron. Soc. 407, 1090–1097 (2010).
  * Ghosh and Lamb (1979a) P. Ghosh and F. K. Lamb, Astrophys. J. 232, 259 (1979a).
  * Ghosh and Lamb (1979b) P. Ghosh and F. K. Lamb, Astrophys. J. 234, 296 (1979b).
  * Freudenreich (1998) H. T. Freudenreich, Astrophys. J. 492, 495 (1998), arXiv:astro-ph/9707340 [astro-ph] .
  * Note (5) In this study we ignore a separate stellar structure called the nuclear bulge that covers the central $\sim$+.1667em200 pc. Around 15% of the GCE emission has been associated to the nuclear bulge stars Macias  _et al._ (2019). Note that the spectral $\gamma$-ray data points we plot are those associated to the Boxy Bulge templateMacias  _et al._ (2019).
  * Ackermann and al. (2014) (Fermi Collaboration) M. Ackermann and al. (Fermi Collaboration), Astrophys. J. 793, 64 (2014), arXiv:1407.7905 [astro-ph.HE] .
  * Harding (2021) A. K. Harding, arXiv e-prints , arXiv:2101.05751 (2021), arXiv:2101.05751 [astro-ph.HE] .
  * Carretti  _et al._ (2013) E. Carretti, R. M. Crocker, L. Staveley-Smith, M. Haverkorn, C. Purcell, B. M. Gaensler, G. Bernardi, M. J. Kesteven, and S. Poppi, Nature (London) 493, 66 (2013), arXiv:1301.0512 [astro-ph.GA] .
  * Gabici  _et al._ (2007) S. Gabici, F. A. Aharonian, and P. Blasi, Astrophys. Space Sci. 309, 365 (2007), arXiv:astro-ph/0610032 [astro-ph] .
  * Abazajian  _et al._ (2020) K. N. Abazajian, S. Horiuchi, M. Kaplinghat, R. E. Keeley, and O. Macias, Phys. Rev. D 102, 043012 (2020), arXiv:2003.10416 [hep-ph] .
  * Note (6) Although, as emphasised by Ploeg and Gordon (2021), a Gaussian smoothing is not a very good approximation of the type of smoothing obtained from MSP natal kicks.
  * Crocker  _et al._ (2011) R. M. Crocker, D. I. Jones, F. Aharonian, C. J. Law, F. Melia, and J. Ott, Mon. Not. Roy. Astron. Soc. 411, L11 (2011), arXiv:1009.4340 [astro-ph.GA] .
  * Crocker (2012) R. M. Crocker, Mon. Not. Roy. Astron. Soc. 423, 3512 (2012), arXiv:1112.6247 [astro-ph.GA] .
  * Morris (2006) M. Morris, in  _Journal of Physics Conference Series_, Journal of Physics Conference Series, Vol. 54 (2006) pp. 1–9.
  * Crocker  _et al._ (2010) R. M. Crocker, D. I. Jones, F. Melia, J. Ott, and R. J. Protheroe, Nature (London) 463, 65 (2010), arXiv:1001.1275 [astro-ph.GA] .
  * Note (7) $\mu_{\rm th}(l,b)$ is associated with the 4 FGL-DR2 catalog; see https://fermi.gsfc.nasa.gov/ssc/data/access/lat/10yr_catalog/.
  * Leane and Slatyer (2019) R. K. Leane and T. R. Slatyer, Phys. Rev. Lett. 123, 241101 (2019).
  * Leane and Slatyer (2020b) R. K. Leane and T. R. Slatyer, Phys. Rev. D 102, 063019 (2020b), arXiv:2002.12371 [astro-ph.HE] .
  * Ruiter  _et al._ (2006) A. J. Ruiter, K. Belczynski, and T. E. Harrison, Astrophys. J. Lett. 640, L167 (2006), arXiv:astro-ph/0511813 [astro-ph] .
  * Perez  _et al._ (2019) K. Perez, R. Krivonos, and D. R. Wik, Astrophys. J. 884, 153 (2019).
  * Hooper  _et al._ (2018) D. Hooper, I. Cholis, and T. Linden, Physics of the Dark Universe 21, 40 (2018), arXiv:1705.09293 [astro-ph.HE] .
  * Hooper and Linden (2018) D. Hooper and T. Linden, Phys. Rev. D 98, 043005 (2018), arXiv:1803.08046 [astro-ph.HE] .
  * Hooper and Linden (2021) D. Hooper and T. Linden, arXiv e-prints , arXiv:2104.00014 (2021), arXiv:2104.00014 [astro-ph.HE] .
  * Albert and al. (2020) (HAWC Collaboration) A. Albert and al. (HAWC Collaboration), Astrophys. J. 905, 76 (2020), arXiv:2007.08582 [astro-ph.HE] .
  * Abramowski and al. (2016) (HESS Collaboration) A. Abramowski and al. (HESS Collaboration), Nature (London) 531, 476 (2016), arXiv:1603.07730 [astro-ph.HE] .

Supplemental Material:
Millisecond Pulsars from Accretion Induced Collapse naturally explain the
Galactic Center $\gamma$-ray Excess
Anuj Gautum, Roland Crocker, Lilia Ferrario, Ashley J. Ruiter, Harrison Ploeg,
Chris Gordon, and Oscar Macias

##  S1 Binary-star evolution (BSE) code

The BSE algorithm presented by Hurley  _et al._ (2002), with updates as
described by Kiel and Hurley (2006) is a state-of-the-art rapid evolution
algorithm that melds single star evolution (Hurley  _et al._ , 2000, 2002)
with binary evolution. BSE evolves binary stellar components through detached,
semidetached, or contact phases of binary evolution. During the evolution of a
binary system, eccentricity may evolve, and, owing to tidal interactions
(which are modeled in detail within BSE), synchronization of stellar rotation
with the orbital motion may occur. BSE includes all orbital angular momentum
loss mechanisms including gravitational radiation and mass loss. Wind
accretion (wherein the secondary accretes some of the material lost from the
primary in a wind), and Roche lobe overflow (where mass transfer occurs if
either star fills its Roche lobe), are modelled on nuclear, thermal, or
dynamical time-scales – whichever is smallest. BSE dynamically adjusts orbital
parameters while taking into consideration any mass variation in the two
stars. Stars with deep surface convection zones and degenerate stars are
unstable to dynamical time-scale mass loss unless the mass ratio of the system
is less than some critical value. The outcome is a common-envelope event
(Paczynski, 1976). This results in either the merging of the two stars, likely
accompanied by the generation of a very strong magnetic field (e.g.
Wickramasinghe  _et al._ , 2014; Briggs  _et al._ , 2015), or the formation of
a close binary. Conversely, mass transfer on a nuclear or thermal time-scale
is assumed to be a steady process. Prescriptions to determine the type and
rate of mass transfer, the response of the secondary to accretion, and the
outcome of any merger events are in place in BSE; details can be found in
Hurley  _et al._ (2000, 2002, 2013).

###  S1.1 Choice of initial binary parameters

Secondary masses: It is standard practice in numerical binary population
synthesis calculations and supported by observations of binary systems in the
Galaxy (Moe and Stefano, 2013; Mazeh  _et al._ , 2003), that companion
secondary masses be described by a flat distribution relative to primary
masses: $0<q=M_{2,ZAMS}/M_{1,ZAMS}<1$.

Initial eccentricity: In our BPS we initialise all binaries with circular
orbits, a standard approach in BPS (e.g. Willems and Kolb, 2002; Hurley  _et
al._ , 2010) and, for low and intermediate-mass stars, a tolerably accurate
description of reality Moe and Di Stefano (2017), especially in light of the
fact that mass transfer and tides in interacting systems – the progenitor of
AICs – quickly bring them to zero eccentricity. In any case, this assumption
implies an error significantly smaller than those accruing from mass
determination uncertainty: see the discussion in Ruiter  _et al._ (2019),
section 4.1, and Appendix S9 below.

##  S2 Common envelope treatment

BSE’s CE evolution follows closely the description of Tout  _et al._ (1997)
where the outcome depends on the initial binding energy, $E_{\rm bind}$, of
the envelope and the initial orbital energy, $E_{\rm orb}$, of the two cores.
A structure parameter $\lambda$ (also called the envelope binding energy
factor) is introduced so that the binding energy of the envelope of a star of
mass $M$, degenerate core mass $M_{c}$ and radius $R$ is (Webbink, 1984)

$E_{\rm bind}=\frac{GM(M-M_{c})}{\lambda R},$

where $G$ is the gravitational constant. Because the ejection of the envelope
may not be complete, an efficiency factor, $\alpha_{CE}$, is introduced. Thus

$\Delta E_{\rm orb}=\alpha_{CE}E_{\rm bind}$

gives the fraction of the orbital energy that is transferred to the envelope.

Most CE evolution models use $\alpha_{CE}=1$ (Willems and Kolb, 2002; Pfahl
_et al._ , 2003) and this is the value that we also use in our modeling.
However, there is a lot of uncertainty in this value (e.g. Tutukov and
Yungelson, 1996) due to the CE evolutionary phase being poorly understood,
mostly because it happens on dynamical time scales and therefore there are no
observations that can be used to constrain the physical parameters that
dominate the CE evolution. The parameter $\lambda$ is of the order of unity
and depends on the stellar density distribution. We have used the BSE
prescription for $\lambda$ calculated with the Cambridge STARS code (Eggleton,
1971; Pols  _et al._ , 1995) and assumed that half of the envelope’s internal
energy aids its ejection. This setup is very similar to “Model 2” of Ruiter
_et al._ (2019).

###  S2.1 Stellar types and a typical evolutionary pathway towards AIC

We use the stellar types ($k$) correspond to the evolutionary phase designated
by Hurley  _et al._ (2000), which are:

$\displaystyle 0$ $\displaystyle\rightarrow\text{ Main Sequence (MS)}\ M\leq
0.7\rm M_{\odot}$ $\displaystyle 1$ $\displaystyle\rightarrow\text{ MS star}\
M>0.7\rm M_{\odot}$ $\displaystyle 2$ $\displaystyle\rightarrow\text{
Hertzsprung gap (HG)}$ $\displaystyle 3$ $\displaystyle\rightarrow\text{ Red
Giant Branch (RG)}$ $\displaystyle 4$ $\displaystyle\rightarrow\text{ Core
Helium Burning (CHeB)}$ $\displaystyle 5$ $\displaystyle\rightarrow\text{
Early Asymptotic Giant Branch (EAGB)}$ $\displaystyle 6$
$\displaystyle\rightarrow\text{ Thermally Pulsing AGB (TPAGB)}$ $\displaystyle
7$ $\displaystyle\rightarrow\text{ Helium Main Sequence (HeMS)}$
$\displaystyle 8$ $\displaystyle\rightarrow\text{ Helium Hertzsprung gap
(HeHG)}$ $\displaystyle 9$ $\displaystyle\rightarrow\text{ Helium Giant Branch
(HeGB)}$ $\displaystyle 10$ $\displaystyle\rightarrow\text{ Helium White Dwarf
(He WD)}$ $\displaystyle 11$ $\displaystyle\rightarrow\text{ Carbon/Oxygen
White Dwarf (CO WD)}$ $\displaystyle 12$ $\displaystyle\rightarrow\text{
Oxygen/Neon White Dwarf (O-Ne WD)}$ $\displaystyle 13$
$\displaystyle\rightarrow\text{ Neutron Star (NS)}$ $\displaystyle 14$
$\displaystyle\rightarrow\text{ Black Hole (BH)}$ $\displaystyle 15$
$\displaystyle\rightarrow\text{ Massless remnant (after SN Ia explosion)}$

####  S2.1.1 A typical evolutionary pathway towards AIC

(The reader is referred to Figure S6). The primary’s ZAMS mass is $7.636\ \rm
M_{\odot}$ and its radius at $t=0$ is $3.364\ R_{\odot}$. The secondary has a
ZAMS mass of $5.702\ \rm M_{\odot}$ and a radius of $2.841\ R_{\odot}$. The
initial separation at the time of the system’s birth is $2369.2\ R_{\odot}$.
The primary evolves through HG, First Giant Branch, CHeB, EAGB and HeGB phases
before becoming a O-Ne WD.

The first interaction occurs at ${\sim}44$ Myr after starburst. Initially, the
ZAMS primary undergoes Roche Lobe Overflow (RLOF) as it expands on the main
sequence and the main sequence secondary accretes matter stably. At ${\sim}54$
Myrs from birth, the primary evolves through the Hertzsprung gap and then
along the asymptotic giant branch, all the while donating mass to the MS
secondary. This is followed by a common envelope (CE) phase between the late
AGB primary ($6.859\ \rm M_{\odot}$) and the main sequence companion ($5.434\
\rm M_{\odot}$). The CE phase causes the binary orbit to shrink and the
primary moves from the EAGB to become a Giant Branch Naked Helium star and
soon evolves into a O-Ne WD.

The next interaction does not occur until there is stable RLOF (which begins
at time $t{\sim}83$ Myrs) between the Hertzsprung Gap secondary ($5.969\ \rm
M_{\odot}$) and the O-Ne WD ($1.315\ \rm M_{\odot}$). The primary continues
accreting mass for another ${\sim}20$ Myrs ($t=103$ Myrs) until the secondary
(now a Hertzsprung Gap Naked Helium star) collapses into a CO WD ($0.683\ \rm
M_{\odot}$). The system remains dormant till $t=108.4$ Myrs, while spiraling
closer ($a=0.194R_{\odot}$) via gravitational radiation. The dormant phase is
followed by the O-Ne WD accreting mass from the CO WD for ${\sim}0.3$ Myrs.
Once the O-Ne WD ($1.439\rm M_{\odot}$,$5.3\times 10^{3}$ s) approaches the
Chandrasekhar limit for rotationally supported WDs, it undergoes AIC to form a
NS ($1.357\ \rm M_{\odot}$, $P=4.8$ ms) at $t=109$ Myrs from birth. After it
is formed, magnetic dipole braking is the dominant spin-down mechanism for the
NS for ${\sim}10^{5}$ yrs, and soon after that, the NS begins accreting mass
from the CO WD ($0.109\ \rm M_{\odot}$).

Figure S6: Schematic diagram showing the main evolutionary stages towards, and
beyond, accretion induced collapse of a white dwarf.

Figure S7: Cosmological time vs. donor star type at time of AIC for all AIC
events in our simulated population of $2\times 10^{9}\rm M_{\odot}$ born in
stars. For more details on stellar types, the reader is referred to Section
S2.1.

Figure S7 shows the number of systems that undergo an AIC event in our
simulated population and their corresponding companion stellar types. Each
line represents an AIC event. Regions of solid colour are regions of high
number density of systems with respective donor types that undergo an AIC
event at any point in the age of the Bulge. Systems with CO WD, He WD and He-
MS donors dominate the AIC progenitor population. We also observe that a
population of NSs formed with Core Helium Burning (CHeB) stars and WDs as
companions have the longest delay times which extend even beyond a Hubble
time; RG and MS donors show similar trends.

##  S3 Angular momentum evolution in AIC and gravitational radiation

The period of the new-born AIC NS due to conservation of angular momentum of
the collapsing WD is given by,

$P_{NS}=\bigg{(}\frac{R_{NS}}{R_{WD}}\bigg{)}^{2}\omega_{WD}.$ (S4)

The angular momentum loss due to gravitational radiation is accounted for
following Watts and Andersson (2002). The total stellar angular momentum $J$
is the sum of the star’s angular momentum $J$ and the canonical angular
momentum of the r-mode $J_{c}$ (Watts and Andersson, 2002):

$J_{c}=-\frac{3\Omega\alpha^{2}\tilde{J}MR^{2}}{2}$ (S5)

where, for an $n=1$ polytrope, $\tilde{J}=1.635\times 10^{-2}$ and $\alpha$ is
the mode amplitude (Lindblom  _et al._ , 1998). The total torque on the star
$\dot{J}$ can therefore be written as

$\dot{J}=\dot{I}\Omega+I\dot{\Omega}+\dot{J_{c}}.$ (S6)

Here, $\Omega$ is the star’s angular velocity. $I=\tilde{I}MR^{2}$ is the
moment of inertia of the star. Because the star is modelled as an $n=1$
polytrope, then $\tilde{I}=0.261$. For this work, we only account for
radiation from the dominant $l=m=2$ current multipole of the r-mode. Following
Watts and Andersson (2002), we assume that the canonical angular momentum
evolves as

$\dot{J_{c}}=-\frac{2J_{c}}{\tau}$ (S7)

where

$\tau^{-1}=-|\tau_{g}|^{-1}+\tau_{b}^{-1}+\tau_{s}^{-1}$ (S8)

and $\tau_{g}$ is the timescales associated with growth of the linear
perturbation due to gravitational radiation back reaction and $\tau_{b,s}$ are
the timescales required for the dissipation of energy due to bulk and shear
viscosity, respectively. The timescales due to the $l=m=2$ current multipole
are:

$\tau_{g}\approx-47M_{1.4}^{-1}R_{10}^{-4}P_{-3}^{6}\ \rm{s}$ (S9)
$\tau_{b}\approx 2.7\times 10^{11}M_{1.4}R_{10}^{-1}P_{-3}^{2}T_{9}^{-6}\
\rm{s}$ (S10) $\tau_{s}=6.7\times 10^{7}M_{1.4}^{-5/4}R_{10}^{23/4}T_{9}^{2}\
\rm{s}$ (S11)

where $M_{1.4}=M/1.4\ \rm{M_{\odot}}$, $R_{10}=R/10\ \rm{km}$, $P_{-3}=P/1\
\rm{ms}$ and $T_{9}=T/10^{9}\ \rm{K}$.

The evolution of the mode amplitude, $\dot{\alpha}$, and the angular velocity,
$\dot{\Omega}$, are expressed as

$\dot{\alpha}=-\alpha\bigg{(}\frac{1}{\tau}+\frac{\dot{\Omega}}{2\Omega}+\frac{\dot{M}}{2M}\bigg{)}\qquad\mbox{and}\qquad\dot{\Omega}=\bigg{(}\frac{\dot{J}}{I}-\frac{\dot{M}\Omega}{M}+\frac{3\tilde{J}\alpha^{2}\Omega}{\tilde{I}\tau}\bigg{)}$
(S12)

##  S4 Spin evolution of neutron stars post AIC

In this section we outline our modelling procedure to compute the spin
evolution of post-AIC NSs in the non-accretion and accretion phases.

###  S4.1 Non accreting neutron stars

If there is no accretion onto the NS, this spins down via magnetic dipole
radiation (e.g., Possenti  _et al._ , 2000),

$\dot{P}_{\rm MB}\simeq 3.1\times
10^{-8}\mu_{30}^{2}I_{45}^{-1}\bigg{(}\frac{P}{\rm
s}\bigg{)}^{-1}\,\,\frac{\rm s}{\rm yr}\,,$ (S13)

where $P$ is the NS rotation period, $\mu_{30}$ is its magnetic moment
($\mu=BR_{NS}^{3}$) in units of $10^{30}\ \rm G\ cm^{3}$ and $I_{45}$ is the
NS’s moment of inertia ($I\simeq 2/5M_{NS}R_{NS}^{2}$) in units of $10^{45}$ g
cm2. We take into account the effects of dipole braking only on the spin
evolution of non-accreting pulsars, whereas, in the case of accreting NSs in
binaries, we also include the torques applied by the in-falling matter on the
NS.

###  S4.2 Accreting neutron stars

Accreting NSs are called low mass X-ray binaries (LMXBs). In these systems the
companion star overfills its Roche lobe and matter escapes through the inner
Lagrange point and flows toward the NS’s magnetosphere. Because this matter
has a large specific angular momentum, it cannot impact directly on to the
surface of the NS but instead forms an accretion disk. Over time, matter in
the accretion disk drifts gradually inwards as angular momentum is transferred
outwards by viscous torques. Once the material reaches the magnetospheric
radius $r_{m}$ where the magnetic pressure dominates over the gas ram pressure
(Gittins and Andersson, 2019),

$r_{m}=\xi\bigg{(}\frac{\mu^{4}}{2GM\dot{M}^{2}}\bigg{)}^{1/7},$ (S14)

the gas is magnetically channelled along magnetic fields lines to impact near
the magnetic poles of the NS. Here $0.5\leq\xi\leq 1.4$ depends on the
fraction of the star’s magnetic flux threading the disk. (Gittins and
Andersson, 2019); in our work we adopt $\xi=0.5$ (Wang, 1996).

At the magnetospheric boundary, the coupling between the field lines and the
disk results in a torque acting on the NS that can spin it up or down
depending on whether $r_{m}$ is larger or smaller than the co-rotation radius
$r_{c}$ given by

$r_{c}=\bigg{(}\frac{GM}{\Omega^{2}}\bigg{)}^{1/3}.$ (S15)

At $r_{c}$ matter rotates with an angular velocity equal to the spin angular
velocity $\Omega=2\pi/P$ of the NS.

Thus, there are two possible cases that need to be considered (Lipunov, 1992;
Boldin and Popov, 2010; Gittins and Andersson, 2019; Ghosh and Lamb, 1979a,
b), namely:

  1. 1.

If $r_{m}<r_{c}$, the standard torque (neglecting magnetic field effects, such
as distortion and/or pinching of the magnetic field lines at the threading
region) is (Ghosh and Lamb, 1979a)

$\dot{\Omega}=\frac{\dot{M}\ r_{m}^{2}\
\Omega_{K}(r_{m})}{I}=\frac{\dot{M}\sqrt{GMr_{m}}}{I}\,,$ (S16)

where $\Omega_{K}(r_{m})$ is the Keplerian angular velocity at $r_{m}$ and $I$
is the moment of inertia of the NS. The disk-star system is thus regarded as
an accretor. Here the gas channelled on to the NS has greater specific angular
momentum than the NS and thus the star is spun up.

  2. 2.

When the accretion disk is truncated at a radius larger than the co-rotation
radius ($r_{m}>r_{c}$), the disk-star system is regarded as a “propeller”. In
this phase, the sign of $(1-\omega_{s})$ dictates whether the NS accretes the
gas and spins up (the “slow rotator” regime, $\omega_{s}<1$) or ejects the gas
and spins down (the “fast rotator” regime, $\omega_{s}>1$ (Wang, 1996)). The
quantity $\omega_{s}$ is the fastness parameter, defined as the ratio of the
spin frequency of the NS ($\Omega$) to the Keplerian orbital frequency
($\Omega_{K}(r_{m})$) at the magnetospheric boundary:

$\omega_{s}=\frac{\Omega}{\Omega_{K}(r_{m})}\,.$ (S17)

The change in spin period is given by,

$\displaystyle\dot{P}$ $\displaystyle\approx$ $\displaystyle-8.1\times
10^{-5}\ \xi^{1/2}\ M_{1.4}^{3/7}\ I_{45}^{-1}\
\mu_{30}^{2/7}\bigg{[}\bigg{(}\frac{P}{1s}\bigg{)}\
\dot{M}_{-9}^{3/7}\bigg{]}^{2}\ (1-\omega_{s})\,\,\rm s\ yr^{-1}\,,$

where $\dot{M}_{-9}=\dot{M}/10^{-9}$ M⊙ yr-1

###  S4.3 Spin period evolution of bulge MSP population

In this section we analyse the effects of the different spin evolution
processes described above. We show in Figure S8 the impact that turning off
accretion after NS formation has on the spin period evolution of a binary. We
find that, while accretion effects may cause some initial enhancement in the
$\gamma$-ray emission post AIC, in the long term they actually keep the MSP
somewhat dimmer than if no post-AIC accretion occurs.

We show in Figure S9 the histograms of the total changes in rotational kinetic
energy caused by the various processes outlined in Section S4 up to the end of
our 14 Gyr simulation for all model AIC MSPs. Note that some systems are
considerably spun up by accretion and that gravitational wave energy losses at
the time of AIC are always small in comparison with all other processes (which
act over multi-Gyr timescales).

Figure S8: Period and prompt $\gamma$-ray luminosity evolution of a typical
AIC NS. Two cases are compared: where the NS is non-accreting or isolated, and
where the NS is in a binary and accretes from its companion.

Figure S9: A distributional comparison of rotational kinetic energy change due
to propeller braking, magnetic dipole braking (in “ejector” phases), accretion
spin-up, and gravitational wave emission for the entire simulated AIC pulsar
population (born from a total mass of $2\times 10^{9}M_{\odot}$ evolved with
MW Bulge star formation history).

##  S5 Further characteristics of the model bulge MSP population

In Figure S10 we show the $P-\dot{P}$ diagram of our pulsar population at t =
13.8 Gyrs simulation time; colors depict the NSs’ magnetic field amplitudes.

Figure S10: $P-\dot{P}$ diagram of the model Galatic bulge AIC pulsar
population

The period and prompt $\gamma$-ray luminosities of the model bulge AIC MSP
population are shown in Figure S11 and Figure S12.

Figure S11: Period distribution of the model population of AIC MSPs (born from
a total mass of $2\times 10^{9}M_{\odot}$ evolved with Galactic bulge star
formation history) at different times as indicated.

Figure S12: Prompt $\gamma$-ray luminosity distribution of the model
population of AIC MSPs (born from a total ZAMS mass of $2\times
10^{9}M_{\odot}$ evolved with Galactic bulge star formation history) at
different times.

Figure S11 and Figure S12 reveal the slow spin-down and dimming experienced by
the bulge population of MSPs over many Gyrs.

##  S6 Initial calculation of model MSP population prompt flux

Our aim is to determine the $\gamma$-ray spectrum received at the Earth due to
our model population of bulge MSPs to assess whether this is compatible with
the observed GCE. In our BPS modeling we evolve the spin period of a
population of AIC NSs over cosmological timescales. The (time dependent)
masses of these NSs are also determined within the model. The magnetic field
strength of each NS is chosen by sampling from the magnetic field distribution
recently determined by Ploeg  _et al._ (2020). This, in turn, allows us to
determine the period derivative attributable to magnetic dipole braking,
$\dot{P}_{\rm MB}$. With $P$ and $\dot{P}_{\rm MB}$ determined for each NS, we
can calculate the instantaneous spin-down power liberated by magnetic braking
as:

$\dot{E}_{\rm MB}=\frac{4\pi^{2}\ I}{P^{3}}\dot{P}_{\rm MB}\,;$ (S18)

note here that $\dot{E}_{\rm MB}(t)$ corresponds to the blue band in Figure 2.

Next, to determine the prompt luminosity of each model MSP and the resulting
total prompt signal, we apply equation 1 together with measurements of the
numerical coefficients obtained by Ploeg  _et al._ (2020) in a combined
statistical analysis of data for nearby, resolved $\gamma$-ray MSPs and the
GCE signal. Specifically, we draw $\log_{10}(f_{\rm\gamma,prompt})$ from a
normal distribution with mean 41.4 and standard deviation $0.6$ and we set
$\alpha_{\gamma}=-1.9$ and $\beta_{\gamma}=0.69$. We then assume the following
spectral form for the prompt $\gamma$-ray spectrum of the $i$th MSP (Ploeg
_et al._ , 2020):

$\frac{dN_{\rm\gamma,prompt,i}}{dE_{\gamma}}=K_{i}E_{\gamma}^{-\Gamma_{i}}e^{(-E_{\gamma}/E_{\rm
cut,i})^{2/3}}\,,$ (S19)

where, $K_{i}$ is a proportionality constant determined from:

$\int_{0.1\rm{GeV}}^{100\rm{GeV}}E_{\gamma}\frac{dN_{\rm\gamma,prompt,i}}{dE_{\gamma}}dE_{\gamma}=F_{\rm\gamma,prompt,i}\equiv\frac{L_{\rm\gamma,prompt,i}}{4\pi
d_{i}^{2}}\,.$ (S20)

In equation S19 $E_{\rm cut,i}$ and $\Gamma_{i}$ are assigned to each MSP by
drawing from normal distributions with mean and standard deviation determined
by Ploeg  _et al._ (2020). Specifically, $\log_{10}(E_{\rm cut}/\text{MeV})$
and $\Gamma$ are drawn from a bivariate normal distribution with a spin-down
power-dependent mean and a correlation coefficient. The distribution of
$\log_{10}(E_{\rm cut})$ has a mean at $a_{E_{\rm
cut}}\log_{10}(\dot{E}/(10^{34.5}\textrm{ erg s}^{-1}))+b_{E_{\rm cut}}$ and a
standard deviation $\sigma_{E_{\rm cut}}$. For $\Gamma$, the mean is at
$a_{\Gamma}\log_{10}(\dot{E}/(10^{34.5}\textrm{ erg s}^{-1}))+b_{\Gamma}$ and
the standard deviation is $\sigma_{\Gamma}$. Here we use $a_{E_{\rm
cut}}=0.23$, $b_{E_{\rm cut}}=2.99$, $\sigma_{E_{\rm cut}}=0.22$,
$a_{\Gamma}=0.48$, $b_{\Gamma}=0.99$, $\sigma_{\Gamma}=0.32$ with a
correlation coefficient between $\Gamma$ and $E_{\rm cut}$ given by
$r_{\Gamma,E_{\rm cut}}=0.73$. These parameters were the highest likelihood
parameters in Ploeg  _et al._ (2020) for a randomly selected Markov chain (see
their model A6). The functional form for the prompt luminosity is that of
Gonthier  _et al._ (2018). The error in the prompt flux estimation arising
from the uncertainties in these parameters is discussed in Appendix S9.

With the prompt luminosity and the spectral distribution of the prompt
emission of each MSP determined, we then assign, as described by Ploeg  _et
al._ (2020), a distance to each MSP $d_{i}$ by Monte Carlo sampling from a
detailed spatial model of the “Boxy Bulge” (Freudenreich, 1998; Macias  _et
al._ , 2019). The weighted mean distance that we obtain with this procedure is
$\langle d_{\rm
MSP}\rangle\equiv\left(\Sigma_{i}L_{i}\Sigma_{j}d_{j}^{2}/L_{j}\right)^{1/2}=7.89$
kpc 555In this study we ignore a separate stellar structure called the nuclear
bulge that covers the central $\sim$ 200 pc. Around 15% of the GCE emission
has been associated to the nuclear bulge stars Macias  _et al._ (2019). Note
that the spectral $\gamma$-ray data points we plot are those associated to the
Boxy Bulge templateMacias  _et al._ (2019).. Furthermore, once we have
assigned a 3D position to each MSP, we can also derive their Galactic
longitudes and latitudes which are required because the flux sensitivity of
Fermi-LAT is dependent on pointing direction (see Section S10).

To recapitulate: We compute the prompt spectral luminosity of each MSP via
equation 1 and solve for $F_{\rm\gamma,prompt,i}$, thence $K_{i}$, to obtain
the prompt differential flux (equation S19) for each MSP. The predicted,
current-day synthetic flux density ($E_{\gamma}^{2}$-weighted) is acquired by
adding the contribution from all MSPs and is shown as the solid blue line in
Figure S13 together with the observed spectrum of the GCE Macias  _et al._
(2019).

##  S7 Determining Inverse Compton and Total Spectra

Unfortunately, unlike the case presented by the prompt emission, the amplitude
and spectrum of the delayed IC emission from each MSP is subject to many
uncertainties. These include the efficiency with which magnetosphere-escaping
$e^{\pm}$ are reaccelerated at MSP wind shocks and/or intra-binary shocks
formed between the pulsar wind and a close companion’s stellar wind, and the
relative energy densities of ISM magnetic field and radiation field (according
to which $e^{\pm}$ energy losses are apportioned into either synchrotron or IC
radiation). On the other hand, this signal is significantly constrained (as
explained below) and we may therefore estimate it using a fitting procedure.
Physically, our model of the IC component assumes emission from a power-law
cosmic ray $e^{\pm}$ population, with an exponential cut-off, IC-scattering a
photon population representative Ackermann and al. (2014) (Fermi
Collaboration) of the inner Galaxy. The normalisation, spectral index, and
cut-off energy of the overall $e^{\pm}$ population are variables in the fit.

###  S7.1 Initial spectrum determination

Given the (preliminary) prompt spectrum is now prescribed up to some error
(see Appendix S9), our first approach to determining the IC spectrum is simply
to find the residuals obtained from subtracting the prompt from the measured
spectral data points, then fit these residuals via an IC emission model; cf.
Figure S13.

The normalisation, spectral index, and cut-off energy of the overall $e^{\pm}$
population are variables in the fit; for the latter two we find best-fit
values via this exercise of $\gamma_{e}=-2.9_{-0.4}^{+0.3}$ and
$\log_{10}({E_{\rm e,cut}/{\rm eV}})=12.1_{-0.1}^{+1.3}$. Here the error
ranges are given by the best-fits to the residuals defined by adopting
extremal values of the prompt flux, i.e., the top and bottom of the blue error
band on the prompt spectrum displayed in Figure S13. These fit values are
physically plausible and thus entirely consistent with a putative MSP origin
of the IC-emitting $e^{\pm}$ (e.g.,Harding (2021)). Fractionally, the overall
luminosity ($>0.1$ GeV) is $0.46^{+0.26}_{-0.19}$ due to prompt emission and
the total $\gamma$-ray luminosity is $2.3^{+3.1}_{-1.1}\times 10^{37}$ erg/s
(cf. the spin-down power liberated by magnetic dipole braking at 13.8 Gyr
according to our BPS of $5.8\pm 0.6\times 10^{37}$ erg/s) and, formally, we
have found a statistically acceptable fit with $\chi^{2}/dof=5.1/(15-3)=0.43$.

Figure S13: Bulge prompt $\gamma$-ray (blue) flux at Earth of the model MSP
population at 13.8 Gyr predicted using equation 1 and adopting central values
of bulge stellar mass, etc. with error band encompassing the overall
fractional error of $1\pm 0.52$ on the prompt flux explained in Appendix S9.
The other curves are for IC (yellow) and total $\gamma$-ray emission (green);
these have been obtained by fitting to the residuals obtained by subtracting
the predicted prompt emission from the measured Macias  _et al._ (2019) GCE
data points. The error bands on the IC and prompt + IC curves encompass the
fits obtained for residuals calculated when assuming values of the predicted
prompt emission at the top and bottom of the blue band.

###  S7.2 Constraints on IC and total $\gamma$-ray emission

Despite achieving a good fit, we note that we should consider and exploit,
where warranted, a number of physical constraints on the fitted IC and total
$\gamma$-ray emission that we have, thus far, ignored. The implications of
these, both singularly and in combination, we now assess in turn. They will
help inform the final spectral decomposition between prompt and IC emission
that we present in the main paper in Figure 3.

Firstly, energy conservation requires that the luminosity of the combined
prompt and IC signals be less than the spin-down power liberated by magnetic
braking. Given the numbers quoted above and with reference to Figure 2, it is
clear that there is sufficient magnetic braking power to energise the total
$\gamma$-ray spectrum, although we come interestingly close in the best-fit
parameter region to saturating this (cf. Sudoh  _et al._ (2020)).

Secondly, IC $\gamma$-rays at 1 GeV are up-scattered from the dominant $\sim$
1.5 eV ISRF of the bulge by parent $\sim$ 13 GeV CR $e^{\pm}$’s. It is
unavoidable that such pairs will synchrotron radiate on the bulge magnetic
field. This has two important consequences: i) the synchrotron luminosity of
the model bulge MSP population must also be accounted for in the overall
energetics and ii) we must assess whether the predicted synchrotron emission
obeys empirical constraints derived from microwave frequency observations of
the bulge (see below).

Thirdly, whilst one might seek to evade these synchrotron constraints by
dialing down the bulge magnetic field amplitude (which we now treat as a free
parameter to be determined), this strategy does not work because i) there is
an empirical lower limit on the magnetic field amplitude of $\sim 9\ \mu$G in
this region Carretti  _et al._ (2013) and ii) the IC signal would then become
too smeared out. In more detail about the latter: because the putative IC-
radiating secondary CR $e^{\pm}$ pairs are launched into the ISM, this
introduces the possibility that, unlike the prompt signal which emerges from
the magnetospheres of individual MSPs, the delayed signal is smeared out over
some angular scale corresponding to the distance over which the pairs are
transported during their energy-loss timescale. We show the equivalent one-
dimensional smearing angular scale for CR $e^{\pm}$ IC-radiating at 1 GeV (and
losing energy via their synchrotron and IC emission) as a function of bulge
magnetic field in Figure S14. These curves have been evaluated for a bulge
ISRF energy density of $\sim$ 2 eV cm-3Ackermann and al. (2014) (Fermi
Collaboration). The different colors illustrate the results obtained for
different assumptions about the transport of radiating pairs in the bulge.

Figure S14: Expected angular transport scales vs. magnetic field amplitude for
CR $e^{\pm}$ IC-radiating 1 GeV $\gamma$-rays off the ISRF of the bulge where
$U_{\rm ISRF,bulge}=2$ eV cm-3Ackermann and al. (2014) (Fermi Collaboration).
The blue curves are for advective transport at 500 km/s nuclear wind, the
yellow for diffusive transport with: i) a CR diffusion similar to that in
plane of the Galaxy Gabici  _et al._ (2007)) and ii) suppressed by $0.01$ with
respect to this value; and green lines are for transport at the Alfvén speed
(i.e., for streaming along the field lines of the bulge magnetic field). The
horizontal line at $0.5^{\circ}$ shows the 68% angular containment radius of
the Fermi-LAT at 1 GeV and the one at $5^{\circ}$ indicates the preferred
smearing scale derived from our template analysis (cf. Figure S16). The
vertical dashed line indicates the lower limit on the magnetic field amplitude
in this region of the Galaxy established in Carretti  _et al._ (2013).

These results can be compared with an empirical determination of the smearing
scale we have made via template analysis of the Fermi-LAT bulge data and
present in Figure S15 and S16. We determined the smearing scale of the boxy
bulge template –used to describe the GCE signal for $|b|\gtrsim 0.5^{\circ}$
in Refs. Macias  _et al._ (2019); Abazajian  _et al._ (2020) – by constructing
a likelihood profile with blurred versions of the original boxy bulge stellar
density map. In particular, we applied both a 2D and a 1D (in the latitude
direction) blurring of the projected boxy bulge map using a Gaussian filter
with different angular scales. Although this approach of applying the Gaussian
filter in 2D on the sky plane ignores the variation (over the line of sight
through the bulge) of the lineal distance corresponding to a given angle, it
provides a reasonable estimate of the range of deviations from the boxy bulge
morphology that are allowed by the GCE data.

The template regression applied here is the same as that introduced in Refs.
Macias  _et al._ (2019); Abazajian  _et al._ (2020). In summary, we performed
a bin-by-bin fitting procedure in which the likelihood function was maximized
at each energy bin independently. To obtain the band fluxes for each template,
we used a power law with a fixed slope of $-2$ and simultaneously varied the
normalization of all the templates in each different energy bin. A description
of all the astrophysical templates included in our bin-by-bin analysis is
presented in Tab. I of Ref. Abazajian  _et al._ (2020). The total log-
likelihood for a given model is obtained by summing the log-likelihoods of all
the energy bins. Figure S15 and S16 show the log-likelihood ratio defined as
$\Delta\log(\mathcal{L})=\log(\mathcal{L}_{BB})-\log(\mathcal{L}_{blurred\;BB})$,
where $\mathcal{L}_{BB}$ is the total log-likelihood of the model containing
the original boxy bulge template, and $\mathcal{L}_{blurred\;BB}$ is the total
log-likelihood of the model containing the blurred boxy bulge.

Fascinatingly, Figure S15 and S16 demonstrate that, while there is evidence
against a 2D smearing of the $\gamma$-ray signal – as would be expected, e.g.,
in the case of isotropic diffusion of the CR pairs or in the case that the
putative astrophysical $\gamma$-ray sources receive isotropic kicks at the
time of formation666Although, as emphasised by Ploeg and Gordon (2021), a
Gaussian smoothing is not a very good approximation of the type of smoothing
obtained from MSP natal kicks. – there is mild evidence for a finite (1D)
smearing in the direction perpendicular to the Galactic plane. In fact, we
find a best fit 1D smearing angular scale of $\sim 5^{\circ}$ degrees, around
700 pc at the bulge distance; the improvement in the overall test statistic at
this angle corresponds to a statistical confidence level of 3.3 $\sigma$
(corresponding to a $2\Delta\log(\mathcal{L})=27.4$ for 15 degrees of
freedom). The tentative detection of this smearing constitutes corroborating
evidence in support of a large IC contribution to the detected GCE
$\gamma$-rays.

Figure S15:  The log-likelihood ratio between “background/foreground”$+$“boxy
bulge” model and the “background/foreground”$+$“blurred boxy bulge” model as a
function of the size of the Gaussian 2D blurring applied to the Boxy Bulge
template. Note that a finite blurring immediately starts to degrade the
statistical significance of the detection of the $\gamma$-ray signal attached
to the template.

Figure S16: The log-likelihood ratio between “background/foreground”$+$“boxy
bulge” model and the “background/foreground”$+$“blurred boxy bulge” model as a
function of the size of the Gaussian 1D blurring applied to the Boxy Bulge
template. The blurring is applied in the vertical axis. Note that such a
blurring in the vertical axis improves the statistical confidence with which a
$\gamma$-ray signal correlated with the template is detected, with the maximum
improvement for a blurring scale of $\sim 5^{\circ}$;
2$\Delta\log(\mathcal{L})\simeq 27.4$ corresponds to 3.3 $\sigma$ for 15
degrees of freedom.

A vertical smearing is consistent with advection in a nuclear wind – for which
there is plenty of independent, multi-wavelength evidence (Crocker  _et al._ ,
2011; Crocker, 2012; Carretti  _et al._ , 2013) – and/or streaming of the CR
pairs at the ion Alfvén speed in the rarefied bulge ISM along a predominantly
poloidal magnetic field, for which there is also independent evidence (at
least for the region quite close to the nucleus; (Morris, 2006)). Given the
measured smearing scale, the scenario of advection in a nuclear outflow with a
characteristic speed of a few 100 km/s requires a magnetic field amplitude of
$\sim 20\ \mu$G; the same magnetic field amplitude works for (vertical)
streaming. The correct smearing scale can also apparently be achieved for a
weak field, $\sim 1-2\ \mu$G, either in the case of energy-dependent (and
necessarily highly-anisotropic given the smearing is 1D) diffusion governed by
a coefficient similar to that measured globally in the Milky Way – or, again
for vertical streaming at the ion Alfvén speed. However, we show below in our
multi-parameter, broadband fitting, that such a low magnetic field amplitude
scenario is disfavoured. In any case, combined analyses of broadband, non-
thermal data covering the inner Galaxy establish that the magnetic field on
the angular scale of the Galactic bulge is $>9\ \mu$G Carretti  _et al._
(2013), rising to $>50\ \mu$G in the inner few hundred pc around the Galactic
center Crocker  _et al._ (2010), ruling out this part of the parameter space.

##  S8 Broadband modeling

To gain further insights and to account mathematically for the constraints
described above, we now perform a classical $\chi^{2}$ analysis wherein we
allow the normalization of the prompt spectrum to float (keeping the spectral
shape fixed) in addition to the normalization, spectral index, and cut-off
energy of the IC-emitting CR $e^{\pm}$ spectrum. Where required, we also
introduce a floating bulge magnetic field amplitude.

In our baseline model, we consider the GCE $\gamma$-ray data points but also
introduce a single additional datum that penalizes the $\chi^{2}$ for
renormalization of the prompt signal amplitude away from the central value we
have already calculated. From this we find a minimum $\chi^{2}/dof$ of
$\sim$0.4 and a $1\sigma$ range for the prompt rescaling of $\sim 0.3-1.9$;
cf. Figure S17.

Figure S17: $\chi^{2}$ per degree of freedom ($dof=12$) for a fit of combined
prompt + IC emission to the GCE data Macias  _et al._ (2019) as a function of
the rescaling of the amplitude of the prompt flux away from the central value
determined in Section S6. For the yellow curve we renormalize the min of
$\chi^{2}/dof$ to 1 in order to establish the $\pm 1\sigma$ confidence
interval on the rescaling parameter.

Note that this preliminary fit does not require either that the implied
synchrotron flux (for the minimum empirically acceptable magnetic field of $9\
\mu$G) obey the empirical upper limits that emerge from microwave observations
Planck Collaboration (2013) nor that the implied total luminosity be less than
the spin down power determined from our BPS. We find, indeed, that over much
of the range of the rescaling parameter within the 1$\sigma$ confidence region
determined above one or both of these constraints is violated.

###  S8.1 MSP contribution to the Microwave Haze

Consider the broadband spectral energy distribution predicted for the best-fit
parameters obtained in the preliminary fitting described above: see Figure
S18.

Figure S18: Spectral energy distribution from radio to $\gamma$-ray energies
for the model bulge MSP population adopting best-fit parameters determine in
Section S8 and a magnetic field amplitude 9 $\mu$G (the empirical lower limit
on the bulge field Carretti  _et al._ (2013)). Curves are as follows:- blue:
prompt emission; yellow: IC; blue: total $\gamma$-ray emission; and red:
synchrotron. The solid curves are evaluated for best-fit parameters, dashed
curves evaluated for parameters at the lower edge of the $1\sigma$ confidence
region (prompt rescaling $\sim 0.2$) and dotted curves are evaluated for
parameters at the top edge (prompt rescaling $\sim 1.9$).

What is tantalising about this figure is that while the predicted synchrotron
radiation is clearly somewhat in excess of the microwave data measured with
the Planck satellite Planck Collaboration (2013), it is nevertheless rather
similar in amplitude. These data points describe the spectrum of the Microwave
Haze Finkbeiner (2004). The haze apparently extends up to somewhat larger
angular scales ($\sim 35^{\circ}$) than the bulge and a morphological
correspondence between it and the even more extended ($\sim 55^{\circ}$) Fermi
Bubbles has been claimed (e.g., Planck Collaboration (2013)). On the other
hand, in its lower reaches the Haze is certainly coincident with the bulge
and, as discussed above, transport effects should carry microwave-emitting CR
$e^{\pm}$ over an angular scale of $\sim 5-10^{\circ}$ after they are launched
by bulge MSPs that are themselves distributed out to $\sim 20^{\circ}$ (cf.
Figure S19). It is also, of course, possible that more than one processes or
population contributes CR $e^{\pm}$ that generate the Haze.

Figure S19: Expected angular transport scales vs. magnetic field amplitude for
CR $e^{\pm}$ synchrotron-radiating at 44 GHz off a bulge magnetic field with
the nominated amplitude and accounting for IC losses on a $U_{\rm
ISRF,bulge}=2$ eV cm-3 ISRF. The blue curves are for advective transport at
500 km/s, the yellow for diffusive transport (with CR diffusion coefficient
suppressed by $\chi=0.01$ with respect to Galactic plane values as motivated
by observation of TeV halos around local pulsars), and green lines are for
transport at the ion Alfvén speed (i.e., for streaming along magnetic field
lines) for $n_{H}=0.01$ cm-3. The vertical dashed line indicates the lower
limit on the magnetic field amplitude in this region of the Galaxy established
in Carretti  _et al._ (2013).

In any case, motivated by these considerations, we perform a new $\chi^{2}$
analysis wherein we introduce the magnetic field amplitude as a further
tunable parameter, and we fit the combination of the Haze microwave data, the
GCE $\gamma$-ray data, a datum for the prompt emission rescaling factor, and a
datum which penalizes the $\chi^{2}$ for the difference between the total
luminosity and the total spin-down power from magnetic braking at 13.8 Gyr
according to our BPS ($1.0\pm 0.1\times 10^{38}$ erg/s). In other words, we
assume here that radiative processes saturate the spin-down power (cf. Sudoh
_et al._ (2020)). The results from this are shown in Figure S20 and subsequent
figures.

Figure S20: $\chi^{2}$ per degree of freedom ($dof=16$) for a fit of prompt +
IC emission to the GCE data Macias  _et al._ (2019) combined with a
synchrotron to the Microwave Haze data Planck Collaboration (2013) as a
function of the rescaling of the amplitude of the prompt flux away from the
central value determined in Section S6. For the yellow curve we renormalize
the min of $\chi^{2}/dof$ to 1 in order to establish the $\pm 1\sigma$
confidence interval on the rescaling parameter.

Not unexpectedly, the region of good fit is now somewhat more constrained to
the region $\sim 0.8-1.7$. The broadband spectral energy distribution obtained
for the new best-fit parameters is shown in the main text as Figure 4. Three
other quantities of interest – the spectral index of the model CR $e^{\pm}$
population, the amplitude of the perpendicular component of the bulge magnetic
field, and the fraction of spin-down power liberated by magnetic braking that
is radiated in non-thermal emission – are shown as functions of the rescaling
in Figure S21, Figure S22, and Figure S23.

Figure S21: Spectral index of the model CR $e^{\pm}$ population as a function
prompt flux rescaling. The vertical dashed lines indicate the best fit and
$\pm 1\sigma$ region.

Figure S22: Amplitude of the perpendicular component of the magnetic field as
a function prompt flux rescaling. The vertical dashed lines indicate the best
fit and $\pm 1\sigma$ region.

Figure S23: Ratio between the bolometric radiated luminosity (in synchrotron +
IC + prompt emission) and the power liberated by magnetic braking at 13.8 Gyr
from the bulge MSP population as a function of the prompt flux rescaling. The
vertical dashed lines indicate the best fit and $\pm 1\sigma$ region.

Figure S24: Ratio between the bolometric power radiated into prompt emission
to all radiation (synchrotron + IC + prompt emission) from the bulge MSP
population at cosmological time 13.8 Gyr as a function of the prompt flux
rescaling. The vertical dashed lines indicate the best fit and $\pm 1\sigma$
region.

These figures demonstrate the overall self-consistency of a scenario that
invokes a combination of both IC and prompt emission from MSPs as an
explanation of the GCE while the CR $e^{\pm}$ launched into the ISM by the
MSPs and requisite to explain the observed IC emission also explain the
observed microwave Haze. Note that the magnetic field amplitude close to the
best-fit region is, self-consistently, that implied by the observed smearing
scale. One question hanging over this analysis concerns the very hard spectral
index, $\sim-2.1-2.2$ of the steady state CR $e^{\pm}$ (required to
simultaneously explain haze and GCE). For reasonable transport and ISM
parameters, this population has to largely cool in situ implying, in the case
of a power law at injection, an injection spectral index of $\sim 1.1-1.2$.
Alternatively, a monochromatic or highly peaked injection spectrum would also
work. Neither of these is excluded by current research into MSPs (see the
discussion in Sudoh  _et al._ (2020); we leave further consideration of this
issue to subsequent work). In any case, the best-fit parameters found in this
analysis are those we adopt in generating the GCE spectral decomposition
illustrated in Figure 3 of the main text.

##  S9 Error budget(s)

The two dominant contributions to the total error in rescaling our results
(for the number of AIC systems) from our BPS modeling of $\sim 6.3\times
10^{7}$ binary systems (equivalent host stellar mass of $2.0\times
10^{9}M_{\odot}$) to the Galactic bulge are due to i) the uncertainty in the
current stellar mass of the bulge and ii) the uncertainty in the bulge
binarity fraction. For the former, Portail  _et al._ (2015) state that the
bulge stellar mass is $(1.4-1.7)\times 10^{10}M_{\odot}$ which we treat as
$(1.55\pm 0.15)\times 10^{10}M_{\odot}$ (1$\sigma$). For the latter, if we
allow that the bulge binarity fraction may range over the (likely extremal)
range $50-100\%$, then the possible fractional change in ZAMS mass of the
parent stellar population relative to our fiducial assumption of 70% binarity
fraction is $\sim 1.0\pm 0.27$. Added in quadrature, the overall fractional
error accruing in scaling from the number of binaries we have modeled to that
appropriate for bulge stellar population is $1\pm 0.30$. Other sources of
error arising from systematic uncertainties in accretion physics on the number
of AIC events, or the correct parameter choices to describe this or, e.g.,
assumptions about initial eccentricity distributions are discussed by Ruiter
_et al._ (2019), section 4.1; these are small in comparison with the
uncertainty accruing from the binarity fraction.

In addition to these two uncertainties, a further significant uncertainty is
encountered in assessing the expected prompt $\gamma$-ray luminosity of the
model MSP population. This stems from uncertainties in the physical values of
the numerical coefficients entering into equation 1
–$f_{\rm\gamma,prompt},\alpha_{\gamma},\beta_{\gamma}$ – which have been
estimated on the basis of a combined statistical analysis by Ploeg  _et al._
(2020) of local, resolved MSPs and the GCE signal. In order to assess the
scale of this uncertainty we generate the distribution of prompt spectra by
simulating 1000 bulge populations of MSPs. For each population we sample at
random a parameter set from the Markov chains produced for model A6 in Ploeg
_et al._ (2020) and assign each pulsar a luminosity, spectral index, $E_{\rm
cut}$ depending on the MSP’s period and the magnetic dipole braking component
of the period derivative. We also simulate a position by sampling from the
boxy bulge density model. The simulated spectrum for the prompt emission of
each population is then the sum of the MSP spectra within the region of
interest. We determine the bin-by-bin means and standard deviations of this
collection of MSP population spectra. From this exercise, the flux-weighted
fractional error averaged over energy bins is $\pm$ 0.43. We add this in
quadrature to the errors due to the binarity fraction uncertainty and the
present-day bulge stellar mass uncertainty to obtain a total (fractional)
prompt $\gamma$-ray luminosity error of $\pm$0.52.

##  S10 MSP detection probability

The Fermi-LAT has a point source flux sensitivity of around $F_{>0.1\rm
GeV}\sim 10^{-12}\ \rm erg\ cm^{-2}\ s^{-1}$, but this varies across the sky
and it is, in general, more difficult to positively identify a $\gamma$-ray
pulsar than a generic $\gamma$-ray point source Abdo and al. (2013) (Fermi
Collaboration). We use the Monte Carlo procedure developed by Ploeg  _et al._
(2020) to estimate the flux for each model MSP such that, given its assigned
sky position and spectrum, it could be resolved as a $\gamma-$ray point source
($F_{\rm sens,point}$) and possibly as an individual MSP ($F_{\rm sens,MSP}$)

Following Ploeg  _et al._ (2020) and references therein, we assume that the
MSP detection threshold sensitivity flux $F_{\rm sens}$ is drawn from a log-
normal distribution (see eq. 2.34 of Ploeg  _et al._ (2020)). Concerning the
threshold sensitivity flux, $F_{\rm sens,MSP}$, for detection and positive
identification of model MSPs as MSPs we assume that $\log_{10}(F_{\rm
sens,MSP})$ is distributed according to a normal distribution with mean
$\log_{10}(\mu_{\rm th})+K_{\rm th}$ and standard deviation $\sigma_{\rm th}$.
Here $\mu_{\rm th}$ is the longitude and latitude dependent Fermi-LAT point
source flux detection threshold 777$\mu_{\rm th}(l,b)$ is associated with the
4 FGL-DR2 catalog; see
https://fermi.gsfc.nasa.gov/ssc/data/access/lat/10yr_catalog/, and $K_{\rm
th}=0.458904$ and $\sigma_{\rm th}=0.251911$ are parameters determined in
Ploeg  _et al._ (2020), invariant across the MSP population; $K_{\rm th}$
accounts for the difficulty inherent in identifying a pulsar as distinct from
a $\gamma$-ray point source (so for calculations related to the probability of
detecting models MSPs as resolved $\gamma$-ray point sources of undetermined
type we set $K_{\rm th}=0$).

##  S11 Model MSP population (sub-threshold) flux distribution

There are ongoing efforts, labeled collectively as point source inference
Collin  _et al._ (2021), to statistically probe the nature of sub-threshold
and individually-unresolvable point sources putatively contributing to the GCE
signal. To date, these efforts have employed a number of different techniques
including non-Poissonian template fitting (“NPTF”; Lee  _et al._ , 2016; Chang
_et al._ , 2020; Buschmann  _et al._ , 2020; Calore  _et al._ , 2021), wavelet
analysis Bartels  _et al._ (2016), and machine learning methods List  _et al._
(2020). All of these techniques at least hint at the existence of a
significant sub-threshold source contribution to the GCE. Moreover, some point
source inference techniques proffer, at least in principle, the ability not
only to evince such sources, but actually to recover their flux distribution,
$dN/dF$. In practice, the results emerging from different groups and/or with
different techniques have not yet converged and do not seem to be yet fully
robust to systematics Leane and Slatyer (2019); Chang  _et al._ (2020); Leane
and Slatyer (2020a, b); Buschmann  _et al._ (2020). In particular, the ability
of such techniques to correctly recover the underlying flux distribution of
putative GCE sources is affected by our lack of a perfect model of the
Galactic foregrounds. Moreover, the statistical reconstruction of a sub-
threshold flux distribution must become increasingly unreliable as the zero
photon limit is approached because such an “ultra-faint” source population is
degenerate with Poissonian emission Chang  _et al._ (2020). The boundary to
this region of degeneracy is roughly given by the single photon flux for which
only a single photon from a given source is expected over the exposure of the
instrument Chang  _et al._ (2020). For Fermi-LAT GCE observations, this
corresponds to $\sim 10^{-11}$ counts cm-2 s-1 (e.g., Chang  _et al._ (2020)).
This is more-or-less where our model MSP distribution peaks (cf. Figure S25)
but given, in particular, that we still predict some brighter sources between
the single photon flux and the point-source sensitivity threshold, the flux
distribution we predict will be subject to empirical testing via point source
inference as these techniques mature. Figure S25 and Figure S26 compare the
predictions of our BPS with the results obtained in the pioneering study of
Lee  _et al._ (2016) and the recent study of Calore  _et al._ (2021). There is
apparently some tension here, but note that the results of Lee  _et al._
(2016) and Calore  _et al._ (2021) are themselves discrepant and, what is
more, the level of disagreement between our results and the point source
inference results is reduced with respect to the more recent study.
Furthermore, while we clearly predict more power in the distribution in some
flux region than determined by Calore  _et al._ (2021), this only occurs for
photon fluxes in the problematic regime around and below the single photon
limit. Also note that, because our distribution cuts off fairly sharply just
above a photon flux of $10^{-10}$ cm-2 s-1, while Calore  _et al._ (2021)
assume a distribution that is continuous to much higher fluxes, these results
are not strictly comparable.

Figure S25: Source count distribution of our model population for $\gamma$-ray
energies 1.9-11.9 GeV and a normalising solid angle of 234 square degrees (as
roughly compatible with Lee  _et al._ (2016).) shown as the blue band; the red
band shows the best-fit source-count function for point sources within the
given solid angle using non-Poissonian template fitting and with 3FGL sources
unmasked as obtained from Lee  _et al._ (2016). The vertical, dashed, blue
line shows approximately the flux equivalent to an expected 1 photon per
source; the vertical, dashed, gray lines subtend the (variable) approximate
point source flux sensitivity threshold Chang  _et al._ (2020).

Figure S26: Source count distribution of our model population for $\gamma$-ray
energies 2-5 GeV and a normalising solid angle of 380 square degrees (as
roughly compatible with Calore  _et al._ (2021).) The orange band shows the
$dN/dF$ distribution obtained by Calore  _et al._ (2021) within the given
solid angle for the fiducial case "pPDF-B, sF-noGCE" (in their notation).
Vertical lines as in Figure S25.

##  S12 Predicted Bulge X-ray phenomenology

Accreting NSs generate X-ray emission. If this is higher than about $10^{36}\
\rm erg/s$, these objects could become detectable as X-ray point sources by
INTEGRAL Haggard  _et al._ (2017). The accretion power (the rate at which
kinetic energy is dissipated at the stellar surface) is given by (e.g. Ruiter
_et al._ (2006))

$L_{X}=\eta\frac{GM_{acc}\dot{M}_{acc}}{R_{acc}},$ (S21)

where $M_{acc}$ and $R_{acc}$ are the accreting NS mass and radius,
respectively, and $\eta$ is the dimensionless, bolometric correction factor
which we determine for each system by randomly sampling from a flat
distribution ranging from 0.1 to 0.6.

Using equation S21 we find a population of ${\sim}6.9\times 10^{3}$ neutron
star accretors with about 60 having current luminosities higher than $10^{36}\
\rm{erg\ s^{-1}}$ (cf. S27). This is entirely consistent with INTEGRAL
observations Haggard  _et al._ (2017) which have uncovered 42 LMXBs and a
further 46 LMXB candidates within $10^{\circ}$ radius around the Galactic
Center.

Figure S27: Number of model Bulge INTEGRAL-observable neutron star accretors
with $L_{X}>10^{36}\rm{erg/s}$ vs. cosmological time. The data points show the
number of identified (green) and identified + candidate (red) LMXBs observed
Haggard  _et al._ (2017) in the MW Bulge.

We also need to establish whether the integrated X-ray emission of the
unresolved point sources of our model is consistent with observations. Thus,
we use the observed Perez  _et al._ (2019) $2-10$ keV X-ray luminosity of the
Galactic Bulge ($6.8\pm 1.2\times 10^{37}\ \rm erg\ s^{-1}$), shown as the
data points in Figure S28 at 13.8 Gyrs. To calculate the model X-ray
luminosity shown in the same figure we use moving average smoothing. The
integrated X-ray luminosity of our model population is, for the central value,
close to, but below, that measured, which is totally acceptable considering
that other X-ray sources may contribute to the total X-ray luminosity of the
Bulge.

Figure S28: Integrated LMXB X-ray luminosity from the simulated Bulge
population over cosmological time. Shown as a band around the simulated X-ray
luminosity time-series is the error due to Bulge mass uncertainties. The data
point shown at t = 13.8 Gyrs is the observed X-ray luminosity from the MW
Bulge (Perez  _et al._ , 2019).

##  S13 TeV halos around bulge MSPs?

An additional argument Hooper  _et al._ (2018); Hooper and Linden (2018, 2021)
that has been raised against an MSP origin for the GCE is that, on the basis
of HAWC data Albert and al. (2020) (HAWC Collaboration) covering a number of
local MSPs, one predicts TeV emission in “halos” around the putative MSP
population in the inner $\sim 0.5^{\circ}$ that would be in excess of that
detected by the HESS instrument Abramowski and al. (2016) (HESS
Collaboration). Formally, because we model here the Boxy Bulge region and not
the Nuclear Bulge (the latter being overlapping with the inner $\sim
0.5^{\circ}$ field), this constraint does not pertain to our model MSP
population. In any case, with $E_{\rm cut}$ typically evaluating to $\sim$ TeV
in our analyses – i.e., with an electron spectrum that cuts off around this
energy – we do not typically predict strong TeV IC emission from our model MSP
population. Moreover, as emphasised in the main text, there are systematic
differences between the bulge and local, disk MSP populations. We will return
to this issue in a subsequent publication.