Cherenkov Telescope Array sensitivity to the putative millisecond pulsar population responsible for the Galactic center excess
Oscar Macias, Harm van Leijen, Deheng Song, Shin'ichiro Ando, Shunsaku Horiuchi, Roland M. Crocker
MMNRAS , 1–22 (2021) Preprint 11 February 2021 Compiled using MNRAS L A TEX style file v3.0
Cherenkov Telescope Array sensitivity to the putative millisecondpulsar population responsible for the Galactic center excess
Oscar Macias, , ★ Harm van Leijen, Deheng Song, Shin’ichiro Ando, , ShunsakuHoriuchi, and Roland M. Crocker Kavli Institute for the Physics and Mathematics of the Universe (WPI), University of Tokyo, Kashiwa, Chiba 277-8583, Japan GRAPPA Institute, Institute of Physics, University of Amsterdam, 1098 XH Amsterdam, The Netherlands Center for Neutrino Physics, Department of Physics, Virginia Tech, Blacksburg, VA 24061, USA Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
The leading explanation of the
Fermi
Galactic center 𝛾 -ray excess is the extended emission froma unresolved population of millisecond pulsars (MSPs) in the Galactic bulge. Such a populationwould, along with the prompt 𝛾 rays, also inject large quantities of electrons/positrons ( 𝑒 ± )into the interstellar medium. These 𝑒 ± could potentially inverse-Compton (IC) scatter ambientphotons into 𝛾 rays that fall within the sensitivity range of the upcoming Cherenkov TelescopeArray (CTA). In this article, we examine the detection potential of CTA to this signature bymaking a realistic estimation of the systematic uncertainties on the Galactic diffuse emissionmodel at TeV-scale 𝛾 -ray energies. We forecast that, in the event that 𝑒 ± injection spectra areharder than 𝐸 − , CTA has the potential to robustly discover the IC signature of a putativeGalactic bulge MSP population sufficient to explain the GCE for 𝑒 ± injection efficiencies inthe range ≈ . 𝐸 − . , a reliable CTAdetection would require an unphysically large 𝑒 ± injection efficiency of (cid:38) 𝑒 ± . Key words: gamma-rays:general – pulsars:general – Dark Matter:general – Cherenkov Tele-scope Array:general – Galactic Center:general
Fermi
Large Area Telescope (
Fermi -LAT) observations (Hooper &Goodenough 2011; Abazajian 2011; Abazajian & Kaplinghat 2012;Gordon & Macias 2013; Macias & Gordon 2014; Calore et al.2015; Daylan et al. 2016; Ajello et al. 2016; Ackermann et al. 2017)of the inner ∼ ◦ of the Galactic center (GC) region have re-vealed extended 𝛾 -ray emission in significant excess with respect tothe astrophysical background model. Early studies of this ‘Galacticcenter excess’ (GCE) appeared to show the signal to be sphericallysymmetric around the GC with a radially-declining intensity. This,together with the fact that the GCE is described by a spectral en-ergy distribution peaked at a few GeV, led to the possibility thatit originates from self-annihilation of weakly interacting massive ★ [email protected] particles (WIMPs) spatially distributed according to something ap-proaching (actually somewhat steeper than) a Navarro-Frenk-White(NFW) density profile (Hooper & Goodenough 2011; Abazajian &Kaplinghat 2012; Gordon & Macias 2013). However, recent reanal-yses of the spatial morphology of the GCE (Macias et al. 2018;Bartels et al. 2018; Macias et al. 2019; Abazajian et al. 2020) havedemonstrated that its spatial morphology is better described bythe (non-spherically symmetric) stellar density distribution of theGalactic bulge than by profiles expected for the WIMP annihila-tion. The Galactic bulge hosts a large variety of stellar populationsfrom old to recently-formed (Garzon et al. 1997; Hammersley et al.2000), and should contain many types of 𝛾 -ray sources, the primeexample of which is the expected large population of millisecondpulsars (MSPs) (Abazajian 2011; Abazajian & Kaplinghat 2012;Gordon & Macias 2013; Ploeg et al. 2017; Gonthier et al. 2018;Ploeg et al. 2020)—old pulsars that have been spun-up due to theirinteraction within a binary system. © a r X i v : . [ a s t r o - ph . H E ] F e b Macias et al.
Pulsars have long been predicted to be sources of electron-positron ( 𝑒 ± ) pairs (Erber 1966; Sturrock 1970, 1971; Aharonianet al. 1995; Atoyan et al. 1995). Highly energetic 𝑒 ± in the windregions or inner magnetospheres of pulsars, also when released intothe general interstellar medium, can inverse-Compton (IC) scatterambient photons to very high energies. The recent HAWC (Abey-sekara et al. 2017) observations of extended TeV-scale 𝛾 -ray emis-sion around Geminga and Monogem have provided indirect ev-idence for TeV 𝑒 ± acceleration in normal pulsar nebulae. Sincea multitude of normal pulsars have been observed in the vicin-ity of the solar system, they have been purported (Hooper et al.2009; Delahaye et al. 2010; Abeysekara et al. 2017; Hooper et al.2017; Profumo et al. 2018; Di Mauro et al. 2019; Jóhannessonet al. 2019) as the likely source of the majority of local cosmicray positrons (Aguilar et al. 2013). Several lines of evidence alsopoint to potentially efficient 𝑒 ± acceleration by MSPs. First, a recentstudy (Sudoh et al. 2020) of the correlation between far-infrared andradio luminosities in star-forming galaxies (SFGs) found that radioemission from MSPs may account for a large fraction of the radio lu-minosity observed in systems with high stellar mass but low star for-mation rate. In particular, the unexpectedly high radio luminositiesof such systems might be explained by the cooling of GeV-scale 𝑒 ± from MSPs via synchrotron radiation in each host galaxy’s interstel-lar medium (ISM). Second, Song et al. (2021) have found evidencefor IC emission produced by MSP populations in globular clustersof the Milky Way. Furthermore, the H.E.S.S. telescope has observedextended TeV 𝛾 -ray emission from the direction of the globular clus-ter Terzan 5 (Abramowski et al. 2013). Such 𝛾 -rays could naturallyarise from Comptonization of the background radiation within theglobular cluster by super-TeV 𝑒 ± (Bednarek et al. 2016). Note thatwhile a similar search conducted by H.E.S.S in another 15 globularclusters (Abramowski et al. 2013), and MAGIC in the M15 globularcluster (Acciari et al. 2019), only produced 𝛾 -ray flux upper limits,it seems likely that forthcoming 𝛾 -ray instruments could resolve alarge number of globular clusters at TeV energies (Ndiyavala et al.2018; Ndiyavala-Davids et al. 2020). Third, several analyses of theGCE (e.g., Horiuchi et al. 2016; Linden et al. 2016; Di Mauro2021) have found evidence for a high-energy tail, possibly relatedto IC emission from 𝑒 ± accelerated by the sources responsible forit. Finally, recent simulation work (Guépin et al. 2020) suggests thatMSPs could efficiently accelerate 𝑒 ± pairs to TeV-scale energies andbeyond.In this paper, we propose that the hypothesized populationof approximately 20 to 50 thousand MSPs (Ploeg et al. 2020) inthe Galactic bulge, could inject large numbers of highly energetic 𝑒 ± into the ISM. The IC signal resulting from these could be de-tectable (Song et al. 2019) by planned TeV-scale 𝛾 -ray telescopessuch as the Cherenkov Telescope Array (CTA). In particular, CTAwill observe the center of the Milky Way with unprecedented spec-tral and spatial detail (Acharya et al. 2018). CTA’s anticipated strat-egy for a survey of the inner Galaxy entails a deep exposure obser-vation of the inner few degrees of the GC region plus an extendedsurvey covering a large fraction of the northern side of the Galacticbulge ( | 𝑙 | (cid:46) ◦ and 0 . ◦ (cid:46) 𝑏 (cid:46) ◦ ). The latter will facilitate thestudy of diffuse very high energy (VHE) sources such as Galac-tic outflows (a possible VHE counterpart to the Fermi bubbles),interstellar gas-correlated 𝛾 -ray emission, a potential dark matter(DM) emission signature (Acharyya et al. 2021), unresolved 𝛾 -raysources (Viana et al. 2020), and, as we will thoroughly explore here,the IC emission from the population of 𝑒 ± launched into the bulgeISM by the MSPs putatively responsible for the GCE.Recently, Song et al. (2019) performed detailed simulations (using galprop V54 Porter et al. 2006) of the IC signature from theputative population of MSPs responsible for the GCE. Comparedto that work, here we recompute all our spatial maps using thelatest version of the code (galprop V56)—which contains newthree-dimensional (3D) models for the interstellar radiation fields,Galactic structures, and interstellar gas maps. Furthermore, nowwe run our simulations in the 3D mode of the galprop framework,abandoning the assumption of Galactocentric cylindrical symmetry.Interestingly, Song et al. (2019) argued that it could be possible forCTA to (i) detect the population of MSPs responsible for the GCE,and (ii) distinguish whether the IC signal emanates from either DMself-annihilation or these unresolved MSPs in the Galactic bulge.Here, we perform a realistic assessment of the sensitivity ofCTA to such an unresolved population of MSPs in the GC. Usingsimulated data, we consider scenarios where the Galactic diffuseemission (GDE) is known accurately, as well as cases where theGDE is mismodeled. In each case, we investigate the necessaryconditions for a reliable CTA detection of the MSPs’ IC signal.Similarly, we study whether CTA could separate the nature of thesource producing the IC emission. As shown in Song et al. (2019),the expected IC emission from the Galactic bulge MSPs has mor-phological differences with respect to the one from DM emission.We now present a systematic study on simulated data that takes intoaccount the latest CTA instrument response function and state-of-the-art models for the GDE to fully address these points.This paper is organized as follows. In Sec. 2 we describe theprocedure used to construct the TeV-scale IC flux maps producedby the putative MSP population responsible for the GCE. In Sec. 3,we present the astrophysical backgrounds in the GC region thatare relevant for the IC searches. In Sec. 4 we provide an overviewof the expected CTA performance, the expected background andsignal rates, and the methodology to estimate the CTA sensitivity.We present our results in Sec. 5, and the discussion and conclusionsin Sec. 6. 𝑒 ± pairs in millisecond pulsars Millisecond pulsars are neutron stars with high rotation rates (spinperiods in the range of tens of ms or less) and relatively low surfacemagnetic fields ( ∼ G). Due to their very rapid rotation, MSPscan produce high electric fields and spin-down power. Despite theirlow surface magnetic fields, their very compact magnetospheresallow for magnetic fields at the light cylinder that are comparableto those of normal pulsars (Harding 2021). Their observed pulsedemission—as seen by an observer located inside the cone scannedby the magnetic field—is explained by the misalignment of the pul-sar rotation axis and the magnetic field axis. If the neutron starsurface temperature ( 𝑇 𝑠 ) is greater than the electron surface bindingtemperature (i.e., 𝑇 𝑠 (cid:38) K), then electrons can be torn from theneutron star surface by the strong electric fields (Michel 1991). Suchelectrons would subsequently move parallel to magnetic field linesand emit 𝛾 -ray photons through curvature radiation. An additionallarge number of “secondary” electrons (and positrons) can be cre-ated through pair production processes at different sites throughoutthe MSPs magnetosphere. Dedicated phenomenological studies ofMSPs demonstrate that most of their high-energy emission occursnear the current sheet outside the light cylinder (Abdollahi et al.2020; Harding 2021) Furthermore, in global magnetosphere mod-els, charge currents, electromagnetic fields, and current sheets scale MNRAS , 1–22 (2021)
TA sensitivity to the putative millisecond pulsar population responsible for the GeV excess with the light cylinder independent of the pulsar surface magneticfield strength or rotation period. Thus, new emission models ofMSPs predict high energy radiation that is similar to that of normalpulsars (the reader is referred to Harding 2021 for further details). A very interesting implication of recent GCE analyses (Macias et al.2018; Bartels et al. 2018; Macias et al. 2019; Abazajian et al. 2020),is that prompt 𝛾 -ray emission from MSPs in the Galactic bulge couldaccount for the bulk of the GCE emission (Abazajian & Kaplinghat2012). Since prompt emission occurs within or very close to themagnetospheres of individual MSPs and 𝛾 rays travel followinggeodesics in space-time, the spatial morphology of the GCE can beused to map the spatial distribution of this putative MSP population.Here, we follow the same approach as Song et al. (2019) whoassumed that the MSP population is smoothly distributed followingthe density of stars in the Galactic bulge (see Sec. 6.2). In par-ticular, the Galactic bulge consists of two main components: thenuclear bulge, and the boxy bulge (see Appendix A). The formercorresponds to the stellar structures residing in the inner ∼
200 pcof the Galactic center (the nuclear stellar disk; and the nuclear starcluster Launhardt et al. 2002), and the latter refers to the stars re-siding in the inner ∼ 𝑒 ± sources are distributed with spherical symmetry as this would bespectrally and spatially similar to the IC signal from DM annihila-tion. More specifically, for the spherical distribution of 𝑒 ± sourceswe use a generalized Navarro-Frenk-White (NFW) profile with amild radial slope. Detailed descriptions of the stellar density maps,and the NFW map used in this work are given in Appendix A.A potentially relevant effect, which our simulations do notconsider, corresponds to the possible smearing effect on the sourcedistribution caused by the MSPs kick velocity distribution. However,it is expected that the kicks experienced by MSPs are lower ( ≈ − ) than for normal pulsars (Podsiadlowski et al. 2005).This comes mainly from the observation that the escape velocityin globular clusters ( (cid:46)
50 km s − , according to Pfahl et al. 2002)is much lower than the average kick velocity of normal pulsars( ∼
250 km s − ; Hobbs et al. 2005; Atri et al. 2019), yet globularclusters have been observed to contain a large number of MSPs.For initial kick velocities (cid:46)
20 km s − , we estimate a negligiblesmearing of the Galactic bulge stellar density maps. However, wewill make a more in-depth investigation of the impact of this effectin a future study of this subject. Given the aforementioned spatial distribution of MSPs in the Galac-tic bulge, Song et al. (2019) thoroughly studied the injection andpropagation of highly energetic 𝑒 ± from MSPs into the interstellarmedium. The focus of that study was the construction of high-resolution IC maps for TeV-scale morphological analyses of future 𝛾 -ray observations. There, it was shown that while the prompt 𝛾 rays from MSPs trace the 𝛾 -ray source distribution, their IC coun-terpart exhibits an energy-dependent spatial morphology. This isdue to the IC emission being the result of the convolution of therelevant photon targets – starlight, infrared (IR) light, the cosmicmicrowave background (CMB) – with the steady state cosmic-ray(CR) distribution, which has itself a spatially-dependent spectrum because it is, in turn, the product of the CR source distribution andenergy-dependent loss and diffusion processes.Similar to the approach of Song et al. (2019) and Ishiwataet al. (2020), here, we solve the CRs transport equation for MSP-accelerated 𝑒 ± using a customized version of the numerical propa-gation code galprop (Porter et al. 2006; Strong et al. 2007). Givena certain CR source distribution, CR injection energy spectrum, andinterstellar medium properties, galprop makes self-consistent pre-dictions for the spatial distributions and spatially-dependent spectraof all CR species. The galprop framework includes pure diffusion,convection, diffusive re-acceleration, and energy losses (e.g., inthe lepton case: Coulomb scattering, synchrotron, inverse Comptonscattering, bremmstrahlung, and nuclear collisions).Song et al. (2019) computed the propagation of 𝑒 ± with gal-prop version 54 (v54). This version of the code makes some simpli-fying assumptions for the propagation of Galactic CRs. Namely, itassumes galactocentric cylindrical symmetry of the CR halo, two-dimensional (2D) interstellar radiation fields (ISRF) model, and 2Dinterstellar gas maps. Importantly, it has been pointed out (Porteret al. 2017; Jóhannesson et al. 2018) that these assumptions couldsignificantly impact the predictions of IC, and gas-correlated 𝛾 -rayemission. In view of this fact, for the present article, we update theIC calculations of Song et al. (2019) by using the latest release of thecode [galprop v56 (Porter et al. 2017; Jóhannesson et al. 2018)].Unlike previous versions of this numerical code, galprop v56contains more sophisticated 3D models of the ISRF, neutral hydro-gen, and molecular hydrogen gas. Such recent improvements of thecode have allowed for detailed predictions of anisotropies (Porteret al. 2017; Jóhannesson et al. 2018) in the Galactic diffuse 𝛾 -ray emission and better, physically-motivated solutions for the CRtransport equation. Currently, there are two alternative 3D ISRFmodels that are available to the user: the first one is based on theGalaxy model proposed by Freudenreich (1998) (F98), and the sec-ond one is based on Robitaille et al. (2012) (R12). Despite that thesetwo models assume different dust densities; and stellar luminosities,their predicted local fluxes are both consistent with near/far-infraredmeasurements. Motivated by previous analyses of the GCE (Maciaset al. 2018; Bartels et al. 2018; Macias et al. 2019), we limit ourstudy to the F98 Galactic structure model within galprop v56.In the present work, we assume a steady-state solution forthe CR transport equation, a homogeneous and isotropic diffusioncoefficient, allow for diffusive reacceleration, and neglect advectionof 𝑒 ± in the Galaxy. All the simulations included in our analysis areperformed with the 3D mode of galprop. We also adopt the CRsource density model labeled as SA50 in Jóhannesson et al. (2018).A brief summary of the main parameter values selected for ouranalysis are shown in Table 1. Jóhannesson et al. (2018) obtainedthese propagation parameters following a similar procedure to thatpresented in Porter et al. (2017). As can be seen in Table 1, for the3D spatial grid we use a bin size of 200 × ×
100 pc . Given thatnuclear bulge region has a radius of ≈
200 pc, this spatial resolutionis far from ideal for modeling the diffusion processes in this region.However, increasing it further would be very difficult due to veryhigh computing memory demands. We note that our simulations arerun in a dedicated computer cluster with 256 GBytes of memoryper node. But, since the code only supports OpenMP environments,it is currently not possible to run fully MPI parallelized simulations.Given that 𝑒 ± could efficiently lose energy through synchrotronradiation, it is important to select well motivated parameters for therandom magnetic field of the Galaxy. We use the results of a modelthat matches the 408 MHz synchrotron data (Strong et al. 2000), andthat agrees with total magnetic field estimates (Heiles 1995; Beck MNRAS , 1–22 (2021)
Macias et al.
Parameter value 𝑋 ℎ [kpc] ± . 𝑌 ℎ [kpc] ± . 𝑍 ℎ [kpc] ± . Δ 𝑋 [kpc] 0 . Δ 𝑌 [kpc] 0 . Δ 𝑍 [kpc] 0 . 𝐷 ,𝑥𝑥 [10 cm s − ] 2 . 𝛿 . 𝑉 Alfven [km s − ] 5 . 𝛾 . 𝛾 . 𝑅 [GV] 3 . 𝛾 ,𝐻 . 𝛾 ,𝐻 . 𝛾 ,𝐻 . 𝑅 ,𝐻 [GV] 4 . 𝑅 ,𝐻 [GV] 200 𝛾 ,𝑒 . 𝛾 ,𝑒 . 𝛾 ,𝑒 . 𝑅 ,𝑒 [GV] 5 . 𝑅 ,𝑒 [GV] 76 Table 1. galprop propagation parameter setup. The parameters ( 𝑋 ℎ , 𝑌 ℎ , 𝑍 ℎ ) , and ( Δ 𝑋, Δ 𝑌 , Δ 𝑍 ) , give the size of the CR halo; and thespatial resolution (grid size) in Cartesian coordinates, respectively. The dif-fusion coefficient is a function of rigidity (i.e., 𝐷 ( 𝑅 ) ∝ 𝛽𝑅 𝛿 ). The injectionfor CR protons is given by 𝑄 ( 𝑅 ) ∝ 𝑅 𝛾 ,𝐻 for 𝑅 < 𝑅 ,𝐻 , 𝑄 ( 𝑅 ) ∝ 𝑅 𝛾 ,𝐻 for 𝑅 ,𝐻 < 𝑅 < 𝑅 ,𝐻 , and 𝑄 ( 𝑅 ) ∝ 𝑅 𝛾 ,𝐻 for 𝑅 > 𝑅 ,𝐻 . Similarly,for CR 𝑒 − s it is given by 𝑄 ( 𝑅 ) ∝ 𝑅 𝛾 ,𝑒 for 𝑅 < 𝑅 ,𝑒 , 𝑄 ( 𝑅 ) ∝ 𝑅 𝛾 ,𝑒 for 𝑅 ,𝑒 < 𝑅 < 𝑅 ,𝑒 , and 𝑄 ( 𝑅 ) ∝ 𝑅 𝛾 ,𝑒 for 𝑅 > 𝑅 ,𝑒 . The injectionspectrum of heavier CR nuclei is written as 𝑄 ( 𝑅 ) ∝ 𝑅 𝛾 for 𝑅 < 𝑅 , 𝑄 ( 𝑅 ) ∝ 𝑅 𝛾 for 𝑅 > 𝑅 . Other parameters included here are the Alfvenvelocity ( 𝑉 Alfven ). See Jóhannesson et al. (2018) for more details. 𝐵 ( 𝑟, 𝑧 ) = 𝐵 exp (cid:18) − 𝑟 − 𝑅 (cid:12) 𝑅 (cid:19) exp (cid:18) − 𝑧𝑧 (cid:19) , (1)where 𝐵 = 𝜇 G, 𝑅 =
10 kpc, and 𝑧 = We assume that the Galactic bulge MSPs are injecting 𝑒 ± into theinterstellar medium with a fraction ( 𝑓 𝑒 ± ) of their spin-down power( (cid:164) 𝐸 ) converted into 𝑒 ± pairs. Similarly, the prompt 𝛾 -ray luminosity( 𝐿 𝛾, prompt ) is assumed to be proportional to the MSPs’ spin-downpower, whose efficiency 𝑓 𝛾 ≡ 𝐿 𝛾, prompt / (cid:164) 𝐸 , is estimated to be about10% on average (Abdo et al. 2013). We can hence write 𝐿 𝑒 ± = 𝑓 𝑒 ± (cid:164) 𝐸 = 𝑓 𝑒 ± 𝑓 𝛾 𝐿 𝛾, prompt (cid:39) 𝑓 𝑒 ± 𝐿 𝛾, prompt , (2)where 𝐿 𝑒 ± is the 𝑒 ± injection luminosity. The 𝑒 ± efficiency 𝑓 𝑒 ± hasbeen estimated to be between approximately 7% and 29% (com-puted for 𝑒 ± energies greater than 10 GeV) from TeV observationsof normal pulsars (Hooper et al. 2017). Interestingly, a recent phe-nomenological analysis of MSPs in globular clusters (Song et al. Model Name Γ 𝐸 cut (TeV)Baseline 2.0 50Inj1 1.5 50Inj2 2.5 50Inj3 2.0 10Inj4 2.0 100 Table 2.
Injection spectra of 𝑒 ± pairs accelerated by a Galactic bulge popu-lation of MSPs. See also (Song et al. 2019). 𝑓 𝑒 ± ≈ 𝑓 𝑒 ± that can be even greater than 90%.In our study, the 𝑓 𝑒 ± is estimated using Eq. (2), simulatedGC CTA observations, and the 𝛾 -ray luminosity from Fermi -LATobservations of the GC. We use the fit results in Macias et al. (2019),which for the boxy bulge obtained 𝐿 bulge 𝛾, prompt = (2.2 ± × ergs − , and for the nuclear bulge 𝐿 NB 𝛾, prompt = (3.9 ± × erg s − .A fit that instead of the stellar templates included a DM (NFW )template gave a luminosity estimate of 𝐿 NFW 𝛾, prompt = (2.7 ± × erg s − . These two different morphologies are considered inorder to investigate the capabilities of CTA to distinguish the twohypotheses under realistic conditions. The maximum energy that 𝑒 ± can attain from MSPs is very uncer-tain. Here we consider several values that have been explored in theliterature (e.g., Yuan & Ioka 2015; Song et al. 2019; Guépin et al.2020). In particular, we assume an injection distribution that is apower law with an exponential cutoff 𝑑 𝑁𝑑𝐸𝑑𝑡 ∝ 𝐸 − Γ exp (− 𝐸 / 𝐸 cut ) , (3)where Γ is the spectral slope and 𝐸 cut is the energy cutoff. Weselect several possible parameter values for Γ , and 𝐸 cut , as shownin Table 2. The spectra are assumed to be the mean of a populationof MSPs in the Galactic bulge.Using Eqs. (2) and (3), we can obtain the source term 𝑞 ((cid:174) 𝑟, 𝐸 ) [with units; MeV − cm − s − sr − ] in the CR transport equation tobe included to our customized version of galprop. Specifically, thiscan be written as the product of the injection spectrum 𝑑 𝑁 /( 𝑑𝐸𝑑𝑡 ) ,and the MSPs density distribution 𝜌 ((cid:174) 𝑟 ) as follows 𝑞 ((cid:174) 𝑟, 𝐸 ) = 𝑐 𝜋 𝑁 𝑑 𝑁𝑑𝐸𝑑𝑡 𝜌 ((cid:174) 𝑟 ) , (4)where the 𝑐 / 𝜋 term is a convention in the galprop code, and 𝑁 is a normalization factor. We normalize the source function in sucha way that its integration over volume and energy matches Eq. (2), 𝑁 ∫ 𝐸 𝑑 𝑁𝑑𝐸𝑑𝑡 𝑑𝐸 ∫ 𝜌 ((cid:174) 𝑟 ) 𝑑𝑟 = 𝐿 𝑒 ± . (5)This procedure is applied for each injection spectrum described inTable 2, the propagation parameter setup shown in Table 1, andthe spatial distributions given by the stellar mass in the Galacticbulge or DM. Figure 1 (second, and third panels of the bottom row)shows the predicted IC morphology for the Galactic bulge, and DMdistribution, respectively. MNRAS , 1–22 (2021)
TA sensitivity to the putative millisecond pulsar population responsible for the GeV excess L a t i t u d e [ d e g ] IC ring 0
IC ring 1
IC ring 2
IC ring 3 L a t i t u d e [ d e g ] Gas-correlated ring 0
Gas-correlated ring 1
Gas-correlated ring 2
Gas-correlated ring 3
Longitude [deg] L a t i t u d e [ d e g ] Fermi bubbles
Longitude [deg]
MSP
Longitude [deg]
Dark matter
Longitude [deg]
Mask
Figure 1.
Predicted spatial morphology of various galactic diffusive emission components in the inner 10 ◦ × ◦ of the Galaxy and at a 𝛾 -ray energy of ≈ − cm − s − sr − ], have a spatialresolution of 0 . ◦ × . ◦ , and were computed with galprop V56. The top row shows the inverse Compton emission produced by background astrophysicalsources, which are divided in four Galactocentric rings (0–3.5, 3.5–8.0, 8.0–10.0, and 10.0–50.0 kpc). The middle row shows the gas-correlated 𝛾 -ray emissionmaps ( 𝜋 + bremsstrahlung) also divided in rings of the same size. From left to right, the bottom row shows: the Fermi bubbles template introduced in Maciaset al. (2019) (which is based on the one obtained in Ackermann et al. (2014)), the predicted MSP IC signal at 10 TeV, the expected IC signal from a sphericaldistribution of 𝑒 ± sources also at 10 TeV (see text for details), and the mask map, respectively. In this section we describe the different templates used to model thebackground and foreground emission in the direction of the GC. Forthis we use two approaches. First, we make predictions using gal-prop v56. Second, we use phenomenological maps whose spectraare extrapolated to match
Fermi -LAT results. The misidentified CRbackground is obtained from detailed simulations made publiclyavailable by the CTA consortium. 𝛾 -ray background The brightest source of extended 𝛾 -ray emission in the directionof the Galactic bulge corresponds to CR protons and 𝑒 − s that aremisidentified as 𝛾 rays. The interaction of highly energetic CRprotons (or heavier nuclei) with the Earth’s atmosphere producesshowers of neutral hadrons which subsequently decay to 𝛾 rays,an important fraction of which cannot be distinguished from astro-physical 𝛾 rays due to the finite rejection power of CTA (Acharyyaet al. 2021; Rinchiuso et al. 2021). Additionally, energetic CR 𝑒 − / 𝑒 + induce electromagnetic showers that are very similar to those pro-duced by 𝛾 rays of astrophysical origin.The main feature of the resulting irreducible 𝛾 -ray backgroundis its isotropy. Even though this component is stronger than anyother in our region of interest (RoI), the goal is that by using amorphological template fit analysis, the impact of this componentcan be significantly reduced. We compute the template associatedto this background by using dedicated simulation studies performedby the CTA consortium (see Sec. 4 for details). The interaction of energetic CRs with interstellar gas, ambient pho-tons, and Galactic magnetic fields produces the brightest source of 𝛾 -rays in Fermi -LAT data. This Galactic diffuse emission (GDE)component is very difficult to model due to significant uncertain-ties in the interstellar gas column density maps, CR spatial/spectraldistributions, and the ISRF model.In this study, we use two different models for the GDE. Thefirst model, referred to as “GDE model 1”, corresponds to one of the
MNRAS000
MNRAS000 , 1–22 (2021)
Macias et al. L a t i t u d e [ d e g ] H I ring 0 H I ring 1 H I ring 2 H I ring 3 Longitude [deg] L a t i t u d e [ d e g ] H ring 0 H ring 1 H ring 2 Longitude [deg] H ring 3 Longitude [deg] L a t i t u d e [ d e g ] Dust reddening E(B-V) positive 42024
Longitude [deg]
Figure 2.
The morphological maps of the gas-correlated 𝛾 -ray emission constructed in Macias et al. (2018). The maps are divided into a atomic (top panels),molecular hydrogen (middle panels), and residual dust (bottom panels) maps which trace the dark neutral material. The hydrogen maps are divided into thesame four Galactocentric rings of Fig. 1. representative models constructed in Jóhannesson et al. (2018). Thesecond model, called “GDE model 2”, assumes the hydrodynamicgas maps constructed in Macias et al. (2018), and the IC templatesused in Abazajian et al. (2020).Specifically, “GDE model 1” assumes the same propagationparameter setup shown in Table 1, but we note that this modeldoes not include the population of Galactic bulge MSPs that ac-counts for the GCE. We run galprop v56 in its 3D mode, anddivide all the predicted TeV-scale 𝛾 -ray maps in four Galactocentricrings. Since the bremsstrahlung and hadronic 𝛾 -ray maps share thesame spatial morphology, we combine these two components intoone (gas-correlated 𝛾 -ray emission). Dividing the IC/gas-correlatedcomponents in different rings allows us to account for the system-atic uncertainties associated to the interstellar medium propertiesand to the somewhat uncertain spatial variation of the CRs. To ob-tain the flux normalization of each GDE component, we fitted thegalprop maps to Fermi -LAT observations of the inner 15 ◦ × ◦ of the GC (Macias et al. 2018), separately varying the IC and gas- correlated rings. Figure 1 (first and second rows) shows the pre-dicted IC, and gas-correlated ring templates at 𝐸 = . | 𝑙 | < ◦ an interpolation method must be used to obtain thegas distribution. Here, instead, we use the hydrodynamic method(Pohl et al. 2008; Macias et al. 2018) which makes direct predic- This is similar to the method used in Rinchiuso et al. (2021), except thathere, the GDE components are divided in different rings.MNRAS , 1–22 (2021)
TA sensitivity to the putative millisecond pulsar population responsible for the GeV excess tions for the gas velocities in the region | 𝑙 | < ◦ , thus providingkinematic resolution in our RoI. Figure 2 displays all the hydro-dynamic templates assumed in this work, which follow the sameannular subdivisions of our galprop templates (Fig. 1).Based on detailed statistical tests with Fermi -LAT observationsof the GC, Macias et al. (2018, 2019) and Buschmann et al. (2020)established that the hydrodynamic maps provide a better fit to thedata than the standard gas maps. Although in this article we are in-troducing these maps with the purpose of estimating the systematicuncertainties in the GDE model, we expect that the hydrodynamicgas maps will be extremely useful once CTA (Acharyya et al. 2021)performs the GC survey.In order to have a physically motivated spectrum for “GDEmodel 2”, we assume the same spectra obtained for each ring ofthe gas-correlated components in “GDE model 1”. The impact ofthis assumption in our uncertainties estimates is expected to besmall given that the gas-correlated maps have very different spatialmorphologies, the templates are divided in rings, and the simulatedCTA data is fitted using an analysis procedure that works bin-by-binin energy.
The
Fermi bubbles (FBs) are giant 𝛾 -ray lobes that were discoveredin Fermi -LAT data (Su et al. 2010) using template fitting techniques.The FBs stretch out to high latitudes above and below the Galacticplane ( | 𝑏 | ≈ ◦ ), while their base (Herold & Malyshev 2019) ispositioned slightly off-set in longitude from the GC ( 𝑙 ≈ − ◦ ). Theyhave an approximately uniform intensity, except for the so-called“cocoon” region in the southern FBs lobe (Su & Finkbeiner 2012;Ackermann et al. 2014). At latitudes | 𝑏 | ≥ ◦ , the FBs intensity iswell described by a flat power-law of the form 𝑑𝑁 / 𝑑𝐸 ∝ 𝐸 − , whichsoftens significantly at energies larger than 100 GeV (Ackermannet al. 2014). At latitudes | 𝑏 | < ◦ (Acero et al. 2016; Ackermannet al. 2017; Storm et al. 2017; Herold & Malyshev 2019), the FBsspectrum also follows a flat power-law but, peculiarly to this region,with no evidence for softening nor an energy cutoff up to energiesapproximately 1 TeV.There is still no consensus on the origin of the FBs. Severalpossibilities that have been discussed in the literature include: recentexplosive outbursts (Su et al. 2010) from the supermassive blackhole Sgr A ★ , and sustained nuclear processes (Crocker et al. 2015)such as star formation activity from the GC. These scenarios requireeither IC, or hadronic gamma-ray emission to explain the observedspectra.In this work, we model the FBs using a similar approach tothe one used by Rinchiuso et al. (2021). In particular, we assumethe best-fit low-latitude FBs spectrum in Ackermann et al. (2017)and then extrapolate it to energies above 1 TeV. We consider twodifferent sets of normalization and energy cutoff, named “FB min”and “FB max” (see Table 3). These are chosen such that the FBsmeasurements (Ackermann et al. 2017) fall within the two models,and they do not overshoot the H.E.S.S diffuse observations of theGalactic ridge (Aharonian et al. 2006). In addition, for the spatialmorphology of the FBs we use the map obtained by Macias et al.(2019). That analysis used an inpainting method to correct for ar-tifacts introduced by the point source mask in Ackermann et al.(2017). Macias et al. (2019) validated the inpainted FBs map witha series of statistical tests. This map is shown in the left bottomcorner of Fig. 1. Model 𝑁 [TeV − cm − s − sr − ] Γ 𝐸 cut [TeV]FB max 1 × − . × − Table 3.
FBs models considered in this work. This component is mod-eled with a power-law with exponential cutoff of the form 𝑑𝑁 / 𝑑𝐸 = 𝑁 ( 𝐸 / ) − Γ exp (− 𝐸 / 𝐸 𝑐𝑢𝑡 ) . The model parameters are the same asin Rinchiuso et al. (2021). We model the point sources in the RoI following the same approachintroduced in Rinchiuso et al. (2021). Namely, we select all thehigh energy point sources included in the third
Fermi high energycatalog (3FHL) (Ajello et al. 2017) that lie within our RoI ( | 𝑙 | ≤ ◦ , | 𝑏 | ≤ ◦ ), and for which the spectrum is given by a simple powerlaw. We note that point sources observed to have an energy cutoff,at GeV-scale energies, in the 3FHL catalog are not included in ouranalysis.Our approach is to mask the aforementioned point sources toavoid potential biases due to extrapolations. We use a disk of radius0.25 ◦ centered at the best-fit position of each selected 3FHL pointsource. In addition, we mask the extended TeV source HESS J1745-303 – one of the brightest sources in our RoI – using a disk of radius0.4 ◦ . Since the angular resolution of CTA will be smaller than 0.1 ◦ ,we anticipate a negligible effect of potential photon leakage. In total,we mask five point sources in our region of interest, representing areduction of approximately 2% of our sky region.Following Rinchiuso et al. (2021), we mask the Galactic planeregion limited by | 𝑏 | ≤ . ◦ , which reduces our sky region by anadditional ≈ . 𝑒 ± propagate over much greater distance scales than the size ofour plane mask. Our total mask is shown in the bottom right cornerof Fig. 1. We adopt the latest publicly available instrument response func-tion (IRF) that is adequate for GC observations (
CTA-Performance-prod3bv1-South-20deg-average-50h.root ). This contains informa-tion of the energy-dependent effective area, point spread function,energy resolution, and irreducible 𝛾 -ray background (in the energyrange 10 GeV–100 TeV). The IRF utilized here was constructed bythe CTA team (Hassan et al. 2015) using a suite of dedicated MonteCarlo simulations. These assume an array of detectors composed of4 large-size telescopes (23 m diameter and sensitive to photons in the20–150 GeV range), 24 medium-size telescopes (11.5 m diameterand sensitive to the 150 GeV–5 TeV range) and 70 small-size tele-scopes (4 m diameter and sensitive to the highest energies). At TeVenergies, CTA has an energy resolution as good as approximately5%. For our analysis we assume the most favorable observation Note that the 0 . ◦ masks end up being just one pixel mask due to thelow resolution of the other astrophysical maps. http://cta.irap.omp.eu/ctools/users/user_manual/response.html MNRAS , 1–22 (2021)
Macias et al. Energy [TeV]10 E d N / d E [ T e V c m s s r ] ICGas-correlatedFB min FB max Gas-correlated (Fermi coll. 2017)IC (Fermi coll. 2017) Energy [TeV]10 C o un t s CR backgroundICGas-correlatedFB min FB max Figure 3.
The predicted 𝛾 -ray spectra for various Galactic diffuse emission components and their corresponding CTA count rates for 500 hours of observationsof the inner 10 ◦ × ◦ of the Galactic center region. The left panel displays the bin-by-bin fluxes of the IC (green points) and gas-correlated (red points)components of the Galactic diffuse emission measured by Fermi -LAT (Ackermann et al. 2017). The green solid and red dashed lines are the spectra predictedby galprop V56. These are normalized to match the
Fermi -LAT observations at GeV-scale energies. The blue dotted line corresponds to the best fit low-latitudeFermi bubbles spectra in (Ackermann et al. 2017) extrapolated to higher energies (this is our FB min model in Table 3). The cyan dotted line is the maximumpossible emission that the FBs can take at TeV-scale energy [FB max model in Table 3—this is the same as in Fig.1 of Rinchiuso et al. (2021)]. The right plotshows the result of convolving the Galactic diffuse emission components in the left panel with the CTA instrument response functions. The irreducible cosmicray background (CR background) is also shown, see Sec. 4.1 for details. conditions; a 500 h on-axis observation from the southern site ata mean zenith angle of 20 ◦ . Additionally, we consider events withenergies between 16 GeV and 158 TeV. The RoI is selected to be asquare of size 10 ◦ × ◦ —centered at Galactic coordinates ( 𝑙, 𝑏 ) = ( ◦ , ◦ ) —which is further binned into pixels of size 0 . ◦ × . ◦ .This is the same binning scheme introduced in Rinchiuso et al.(2021). These authors showed that the choice of bin size had noimpact on the results given the high photon statistics obtained ineach spatial bin. We note that all our background and signal modelsconstructed with galprop v56 also have a resolution of 0 . ◦ × . ◦ .Increasing the resolution further is limited by computational costs. In order to get the expected counts for a given sky model, we con-volve the signal/background templates with the CTA IRFs describedabove. The convolution is done with the Gammapy (Deil et al. 2017;Nigro et al. 2019) analysis tools, and the model templates are shownin Fig. 1 and 2.Since the CTA’s point spread function is smaller than ourpixel size (0 . ◦ x 0 . ◦ ), it can be neglected in our calculations.The function that describes the expected counts Φ 𝑚𝑖 𝑗 for a certain This could be obtained with an optimized observation strategy plannedby the CTA consortium (Acharyya et al. 2021). https://gammapy.org/ astrophysical component 𝑚 , at the 𝑖 𝑡ℎ longitudinal, 𝑗 𝑡ℎ latitudinal,and Δ 𝐸 energy bin can then be written as Φ 𝑚𝑖 𝑗 = 𝑇 𝑜𝑏𝑠 ∫ Δ 𝐸 𝑑𝐸 𝛾 ∫ ∞ 𝑑𝐸 (cid:48) 𝛾 𝑑𝜙 𝑚𝑖 𝑗 𝑑𝐸 𝛾 𝐴 𝛾𝑒 𝑓 𝑓 ( 𝐸 (cid:48) 𝛾 ) 𝐷 ( 𝐸 𝛾 , 𝐸 (cid:48) 𝛾 ) , (6)where 𝐴 𝛾𝑒 𝑓 𝑓 ( 𝐸 𝛾 ) is the energy-dependent effective area, 𝐷 ( 𝐸 𝛾 , 𝐸 (cid:48) 𝛾 ) is the energy dispersion function, and 𝑑𝜙 𝑚𝑖 𝑗 / 𝑑𝐸 𝛾 is theincoming flux spectrum from the 𝑚 source. The whole functionis integrated over the energy bin Δ 𝐸 𝛾 , and multiplied by the totalobservation time 𝑇 𝑜𝑏𝑠 . The two different energies stand for the trueincoming energy 𝐸 (cid:48) 𝛾 , and the reconstructed energy 𝐸 𝛾 . Assuminga homogeneous sky exposure in our RoI, we set 𝑇 𝑜𝑏𝑠 =
500 hoursin our analysis.As an example of the IRFs convolution results, we show inFig. 3 (left) the predicted spectra for each of the background modelcomponents, and in Fig. 3 (right) the same background templatesafter convolution with the CTA IRFs. We have carefully checkedthat our pipeline reproduces well previous works in the literature[e.g., (Rinchiuso et al. 2021)].
The conventional analysis method for TeV-scale 𝛾 -ray observationsinvolves selecting two regions of the sky with approximately thesame backgrounds, but different expected signals. The region withlarger signal is called the “ON” region, while the other the “OFF” The irreducible CR background is formally included as an additionalcomponent 𝑚 in Eq. 6. MNRAS , 1–22 (2021) TA sensitivity to the putative millisecond pulsar population responsible for the GeV excess region. In this approach, the hypothesis testing is done using a teststatistic (TS) defined as the difference in photon counts betweenthe ON and OFF regions. However, the study by Silverwood et al.(2015) demonstrated that the template-fitting procedure—which isthe standard analysis method for e.g., Fermi -LAT data—greatlyimproves the CTA sensitivity to extended GC signals because itallows for the full exploitation of the morphological differencesbetween background and signal templates.In this work, we adopt a template-fitting approach to studyCTA sensitivities to a putative IC signal from Galactic bulge MSPs.In particular, we divide the mock data in 11 bins logarithmicallyspaced from 16 GeV to 158 TeV. This makes the bins larger thanthe energy resolution of the CTA data. For each independent energybin Δ 𝐸 , we define the likelihood function as L( 𝝁 | 𝒏 ) = (cid:214) 𝑖 𝑗 𝜇 𝑛 𝑖𝑗 𝑖 𝑗 𝑒 − 𝜇 𝑖𝑗 𝑛 𝑖 𝑗 ! (7)where 𝒏 = { 𝑛 𝑖 𝑗 } is the simulated data, and the model 𝝁 = { 𝜇 𝑖 𝑗 } isa linear combination of model templates 𝜇 𝑖 𝑗 ( 𝜶 ) = ∑︁ 𝑚 𝛼 𝑚 Φ 𝑚𝑖 𝑗 , (8)with Φ 𝑚𝑖 𝑗 as given in Eq. 6, and the flux normalizations 𝛼 𝑚 are thefitting parameters. The optimization is done with the minuit (Nelder& Mead 1965) algorithm contained in the iminuit package. For each independent energy bin, we simultaneously fit the fluxnormalizations of all the background and signal components withinour 10 ◦ × ◦ and | 𝑏 | ≥ . ◦ RoI. The normalizations of the sourcesconsidered are insensitive to the spectral shape assumed at eachenergy bin due to the bins being small. This fitting approach, alsoknown as a bin-by-bin analysis, has been widely used in the analysisof
Fermi -LAT data (e.g., Ackermann et al. 2015). The advantageof this method over the more traditional broad-band analysis (e.g.,Rinchiuso et al. 2021; Acharyya et al. 2021), is that the bin-by-binmethod is much less affected by assumptions about the spectralshape of the model templates.Since CTA is not in operation yet, we create synthetic data foreach energy bin and each component. We do this by drawing froma Poisson distribution with mean 𝜇 𝑖 𝑗 ( 𝛼 𝑚 ) , and then summing theresulting maps to obtain the total mock data set for an energy bin.We evaluate the significance of the MSPs hypothesis at eachenergy bin Δ 𝐸 𝑘 using a test statistic defined asTS k = − (cid:18) L( 𝝁 , ˆ 𝜽 | 𝒏 )L( ˆ 𝝁, ˆ 𝜽 | 𝒏 ) (cid:19) (9)where 𝝁 are the normalizations of the background-only hypothesis,and ˆ 𝝁 and ˆ 𝜽 are the best-fit parameters under the background plusMSPs hypothesis. The total TS of the MSPs IC template can beobtained as TS = (cid:205) 𝑛𝑘 = TS 𝑘 where TS 𝑘 is given in Eq. 9, and thesum runs up to 𝑛 =
11 (number of energy bins). In our case, theMSPs IC template has 11 degrees of freedom (flux norm at eachenergy bin). Hence, in order to get the p-value of this template we https://iminuit.readthedocs.io/en/stable/ are required to use the mixture distribution formula (Macias et al.2018) 𝑝 ( TS ) = − 𝑛 (cid:32) 𝛿 ( TS ) + 𝑛 ∑︁ 𝑘 = (cid:18) 𝑛𝑘 (cid:19) 𝜒 𝑘 ( TS ) (cid:33) (10)where (cid:0) 𝑛𝑘 (cid:1) is the binomial coefficient with 𝑛 = 𝛿 is the Dirac deltafunction, and 𝜒 𝑘 is a 𝜒 distribution with 𝑘 degrees of freedom. Wecan compute the detection significance in 𝜎 units corresponding tothe addition of one new MSP IC norm parameter by using the totalTS value and the p-value shown in the above equation. Specifically,we do this with the following recipe (Macias et al. 2018)Number of 𝜎 ≡ √︂ InverseCDF (cid:16) 𝜒 , CDF (cid:2) p ( TS ) , ˆTS (cid:3)(cid:17) (11)where (InverseCDF) CDF is the (inverse) cumulative distribution.The first argument of each of these functions is the distributionfunction and the second is the value (inverseCDF) at which theCDF is evaluated. The total TS value is denoted by ˆTS. From Eq. 11we obtain that a 5 𝜎 detection corresponds to TS = .
1. It followsthat we may claim a detection of the IC template for total TS valueslarger than this threshold value.
Having introduced the fitting procedure and explained how we cre-ate the CTA simulated data, we now present the results of our CTAsensitivity analysis. Specifically, we start by working out the num-ber of independent gas-correlated rings for which we can get stablefits for the background only hypothesis. Then, we present the mini-mum IC (and 𝑒 ± ) luminosities necessary for a reliable MSPs signaldetection, and the ability of the method to accurately distinguishbetween the MSPs and dark matter origin of the radiating 𝑒 ± . Fi-nally, we show the impact of the Galactic diffuse emission modeluncertainties on our results. To the best of our knowledge, this is the first time that a Galacticdiffuse background model divided in multiple galactocentric rings(see Fig. 1) is used for TeV-scale 𝛾 -ray analyses. Though this methodhas been implemented with great success in studies of Fermi data[e.g., Abdollahi et al. (2020)], it is not a priori obvious that the sametechnique can be applied to CTA, given its smaller field of view.To address this concern we performed a fitting procedure wherewe carefully checked for potential degeneracies between the differ-ent gas-correlated and IC rings shown in Fig. 1. From this test,we obtained that the morphological differences between the fourgalactrocentric IC emission maps were not significant enough forthe pipeline to distinguish them in the simulated data. This is be-cause most of these spatial differences lay in the parts that are furtheraway from our RoI (central 10 ◦ × ◦ and | 𝑏 | ≥ . ◦ of the Galaxy).We therefore decided to combine the four IC rings into one singlemap and only keep the gas-correlated emission maps split into fourgalactocentric rings.It is worth noting that using different IC rings for GC analy-ses could still be viable with a different observational strategy. Inparticular, the GC survey plan proposed in Fig.1 of Acharyya et al.(2021) covers a region that is almost two times larger than the RoI MNRAS000
Having introduced the fitting procedure and explained how we cre-ate the CTA simulated data, we now present the results of our CTAsensitivity analysis. Specifically, we start by working out the num-ber of independent gas-correlated rings for which we can get stablefits for the background only hypothesis. Then, we present the mini-mum IC (and 𝑒 ± ) luminosities necessary for a reliable MSPs signaldetection, and the ability of the method to accurately distinguishbetween the MSPs and dark matter origin of the radiating 𝑒 ± . Fi-nally, we show the impact of the Galactic diffuse emission modeluncertainties on our results. To the best of our knowledge, this is the first time that a Galacticdiffuse background model divided in multiple galactocentric rings(see Fig. 1) is used for TeV-scale 𝛾 -ray analyses. Though this methodhas been implemented with great success in studies of Fermi data[e.g., Abdollahi et al. (2020)], it is not a priori obvious that the sametechnique can be applied to CTA, given its smaller field of view.To address this concern we performed a fitting procedure wherewe carefully checked for potential degeneracies between the differ-ent gas-correlated and IC rings shown in Fig. 1. From this test,we obtained that the morphological differences between the fourgalactrocentric IC emission maps were not significant enough forthe pipeline to distinguish them in the simulated data. This is be-cause most of these spatial differences lay in the parts that are furtheraway from our RoI (central 10 ◦ × ◦ and | 𝑏 | ≥ . ◦ of the Galaxy).We therefore decided to combine the four IC rings into one singlemap and only keep the gas-correlated emission maps split into fourgalactocentric rings.It is worth noting that using different IC rings for GC analy-ses could still be viable with a different observational strategy. Inparticular, the GC survey plan proposed in Fig.1 of Acharyya et al.(2021) covers a region that is almost two times larger than the RoI MNRAS000 , 1–22 (2021) Macias et al. assumed in our work. Another interesting possibility could be to usethe CTA divergent pointing mode (Gérard 2016) in which CTA cansurvey a region as large as 20 ◦ × ◦ [see for example the “Deepexposure scenario” proposed in Coronado-Blázquez et al. (2021)].We leave these interesting alternatives for future studies, and stickonly with our survey strategy, which was suggested in Rinchiusoet al. (2021). Given that the systematic uncertainties in the GDE model is oneof the most difficult problems for GeV-scale 𝛾 -ray analyses of theGC, it is reasonable to assume that this component will also bevery challenging for forthcoming CTA observations of similar skyregions. With this in mind, we started our sensitivity analysis bytesting how well our pipeline recovers the properties of the simu-lated MSP IC signal in different case scenarios. In particular, weconsidered various alternative GDE models, different signal spec-tra and spatial morphologies, and we further employed a method tostudy the impact of mismodeling the Galactic diffuse backgrounds.Details of each of our GDE model components are given in Sec. 3.As for the MSPs IC signal, we considered a range of physicallyreasonable injection spectra presented in Table 2.Figure 4 shows the results of our signal recovery tests in thecase where the diffuse model is given by “GDE model 1” plusFB min , and we assume that the backgrounds are perfectly modeled;this is accomplished by fitting the mock data with the same tem-plates used in the generation of the simulations. This figure is madefrom 5000 realizations of synthetic data that contains irreduciblecosmic ray background photons plus astrophysical 𝛾 rays sampledfrom the aforementioned diffuse model. Each row corresponds tothe five different injection models introduced in Table 2. The lumi-nosity of the signal that is injected (red solid line) is displayed alongwith the 68% containment on the recovered luminosity (blue region)in the right hand side panels. We also include the 5 𝜎 (TS = . 𝛾 -ray luminosity increases very similarlyfor all the injection models that were included (see Table 2).In the left panels of Fig. 4 we present the TS distributions ofthe signal templates. We stress that while fitting the mock data weallow all components (background and signal templates) to varyin the fits. The filled regions denote the TS distributions of theMSPs’ IC signal (light blue) and their mean value (blue solid line)is also displayed. These demonstrate that with 500 h of GC (central10 ◦ × ◦ and | 𝑏 | ≥ . ◦ of the Galaxy) observations with theCTA and an accurate knowledge of the astrophysical backgrounds(corresponding to the most optimistic scenario considered in thiswork) the technique presented here should be able to detect the ICsignal from the putative GC MSP population, thereby allowing usto constrain the source populations generating the high-energy 𝛾 rays in this sky region.For this same case scenario, we present details of the recoveredspectra in Fig. 5. Each panel in this figure shows the characteristicsof the spectra that are recovered with a 5 𝜎 detection significance.The five different panels correspond to each 𝑒 ± injection modelpresented in Table 2. As explained in Sec. 2, we propagate such 𝑒 ± within the galprop V56 framework, then produce IC spectral templates (red solid lines), and lastly inject these signal maps intothe simulated data. We obtained the results in this figure by fittingthe mock data with a bin-by-bin fitting procedure that allows us towork out the fluxes and corresponding 68% confidence intervals(blue filled regions) at each independent energy bin. Fluxes largerthan those shown in these panels should be successfully detectedby CTA, provided the parent 𝑒 ± ’s are injected by an unresolvedpopulation of MSPs (tracing the distribution of stellar mass in theGalactic bulge) and the Galactic diffuse background is perfectlymodeled. We also present a summary of the minimum 𝐿 𝛾, IC (and 𝐿 𝑒 ± ) required for a CTA detection in Table 4. In the previous section (Sec. 5.2) we created mock data by samplingfrom the “GDE model 1”, FB min , and the irreducible CR back-ground. We then injected various different MSP IC signals into thedata and applied a fitting procedure to recover the injected signals.In particular, the fit included all the templates used in the generationof the mock data.In this subsection, we applied the same pipeline, except thatthis time we also added to the fit an IC template generated by DMemission. Namely, we injected a MSP IC signal into the mock data,and subsequently attempted to recover it by including both the MSPIC and DM IC templates in the fit. The main objective of this testis to figure out the conditions under which CTA would be able todisentangle a new extended 𝛾 -ray source in the GC based on themorphological characteristics of the IC radiation emitted by thesource.We present the results of the tests for degeneracy betweenthese two competing hypotheses in Fig. 6. The left panel showsthe TS distribution of the DM IC template as a function of theinjected MSPs IC luminosity. In this panel, we display the meanTS values (green solid line) and their respective 68% containmentband (green filled region). As can be seen, for all the evaluatedluminosities, the DM IC template was found to have TS (cid:46)
10 (ora statistical significance of (cid:46) . 𝜎 for 11 degrees of freedom).The right panel shows the injected MSP IC luminosity versus theluminosities recovered for the MSPs IC (blue filled region) and theDM IC (green filled region) templates, respectively. It is clear thatfor IC luminosities larger than the minimum 𝐿 𝛾 ’s given in Table 4(see the row corresponding to FB min , mismodeling of the GDE),only a small fraction of the injected MSP IC luminosity is absorbedby the DM IC template. Due to uncertainties in the Galactic diffuse background, it will bechallenging for future analyses of actual CTA data to model thiscomponent perfectly (Acharyya et al. 2021). We re-create this real-world situation by constructing simulated data with “GDE model1” and analyzing it with “GDE model 2”, which allows us to testwhether mismodeling of the GDE could originate in a false positivedetection of an MSP IC signal. This methodology is inspired by arecent study on simulated
Fermi data by (Chang et al. 2020), anda similar one performed by the CTA consortium (Acharyya et al.2021) in the context of DM searches in the GC.We thus repeated the same analyses performed in Sec. 5.2and Sec. 5.3, but this time, mimicking the mismodeling of theGalactic diffuse emission as described above. We show the results,
MNRAS , 1–22 (2021)
TA sensitivity to the putative millisecond pulsar population responsible for the GeV excess T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Baseline5 significance T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj15 significance T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj25 significance T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj35 significance Injected Luminosity [erg/s] T e s t s t a t i s t i c
5 significance Injected Luminosity [erg/s] R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj45 significance FB min , Perfect GDE Figure 4.
The CTA sensitivity to the IC signal from an unresolved population of MSPs tracing the distribution of stellar matter in the Galactic bulge. Thesetests assume a perfect knowledge of the Galactic diffuse emission. The latter is modeled with a combination of the FB min model, and “GDE model 1” (seealso Sec. 3.2) .
Left panels:
Detection significance (Test statistic) of the MSPs’ IC signal for a given IC injection luminosity. Each row assumes a different 𝑒 ± spectrum model (Inj1,..., Inj4) shown in Table 2. The 5 𝜎 detection threshold is displayed as a green dotted line (see Eqs. 10 and 11 for details). A summary ofthe minimum 𝐿 𝛾,𝐼𝐶 required for a CTA detection is shown in Table 4. The blue solid line represents the mean of the results, while the light blue region givesits variance. Right panels:
Comparison of the recovered signal luminosity with the injected one. The diagonal red line represents the ideal case in which theextracted signal matches the injected signal perfectly. The blue region shows the recovered IC signals. We used 5000 realizations of synthetic data that containsirreducible cosmic-ray background photons plus astrophysical 𝛾 rays sampled from the aforementioned diffuse model.MNRAS000
Comparison of the recovered signal luminosity with the injected one. The diagonal red line represents the ideal case in which theextracted signal matches the injected signal perfectly. The blue region shows the recovered IC signals. We used 5000 realizations of synthetic data that containsirreducible cosmic-ray background photons plus astrophysical 𝛾 rays sampled from the aforementioned diffuse model.MNRAS000 , 1–22 (2021) Macias et al. E d N / d E [ T e V c m s s r ] injected baseline signalextracted signal injected signal 1extracted signal E d N / d E [ T e V c m s s r ] injected signal 2extracted signal Energy [TeV] injected signal 3extracted signal Energy [TeV] E d N / d E [ T e V c m s s r ] injected signal 4extracted signal FB min , Perfect GDE Figure 5.
Minimum flux to detect the IC signal from a putative population of MSPs responsible for the GCE. The red line shows the injected IC signal whilethe blue regions display the 68% confidence intervals on the recovered signal. The normalization of the spectra corresponds to the 𝛾 -ray luminosity that wouldbe detected with 5 𝜎 significance (see also the green dotted line in Fig. 4) for 500 h of CTA observations of the central 10 ◦ × ◦ and | 𝑏 | ≥ . ◦ of the Galaxy.Higher fluxes than the ones displayed here would be detected by CTA with a statistical significance larger than 5 𝜎 . respectively, for the signal recovery test in Fig. 7, the minimum fluxrequired for a 5 𝜎 significance detection of the MSPs IC signature inFig. 8, and lastly, the degeneracy tests between the MSPs and DMhypotheses in Fig. 9.We found that, in the case where the Galactic diffuse back-ground is mismodeled (“FB min , mismodeling of the GDE” sce-nario), the CTA sensitivity to the MSPs IC signal is reduced byapproximately 38%, 33%, 3%, 4%, and 30%, respectively, for the 𝑒 ± injections scenarios Baseline , Inj1 , Inj2 , Inj3 , and
Inj4 (see Ta-ble 2). A summary of the minimum 𝐿 𝛾, IC required for a CTAdetection is given in Table 4. Fermi bubbles model
We have evaluated the impact of the FBs on our sensitivity byassuming the “FB max” model presented in Table 3 and the per-fect GDE model scenario. The results of this test are presented inFigs. 10, 11, and 12; which can be directly compared Figs. 4, 5, and6, respectively.From these comparisons, it follows that the strongest impacton the CTA sensitivity to a MSPs IC signal in the GC region isdue to assumptions on the FBs model. We found a degradation ofthe sensitivity of approximately one order of magnitude at worst—which corresponds to the MSP 𝑒 ± injection model Inj1 in Table 3—
MNRAS , 1–22 (2021)
TA sensitivity to the putative millisecond pulsar population responsible for the GeV excess Injected Luminosity [erg/s] ) T e s t s t a t i s t i c Detection threshold (5 ) Injected Luminosity [erg/s] R e c o v e r e d L u m i n o s i t y [ e r g / s ] injected signalrecovered MSP signalRecovered DM signal FB min , Perfect GDE Figure 6.
Tests of spectro-morphological degenerecies between the MSPs IC and DM IC templates with simulated CTA observations of the GC region. Theastrophysical background is sampled from "GDE model 1", FB min , and the irreducible CR background (see also caption of Fig. 4). A MSPs IC signal is injectedinto the mock data and subsequently a bin-by-bin fitting procedure is applied using the same background templates used in the generation of the mock data, inaddition to the MSPs IC and a DM IC templates. The left panel shows the mean TS distribution (green solid line) and corresponding 68% confidence region(gree filled area) for the DM IC template. The right panel displays the fraction of MSPs IC (blue filled area) and DM IC (green filled area) luminosities thatare obtained after injection of only the MSPs IC signals of various luminosities.Minimum 𝐿 𝛾, IC for detection [erg s − ] Baseline Inj1 Inj2 Inj3 Inj4 FB min , perfect GDE.1 . × . × . × . × . × FB min , mismodeling of the GDE.1 . × . × . × . × . × FB max , perfect GDE.7 . × . × . × . × . × FB max , mismodeling of the GDE.9 . × . × . × . × . × Table 4.
Minimum IC luminosity required for a 5 𝜎 significance detectionwith CTA. The luminosities are computed in the energy range 16 GeV and158 TeV. The columns correspond to different MSPs 𝑒 ± spectra consideredin Table 2. The corresponding minimum 𝐿 𝑒 ± can be obtained by realizingthat the 𝐿 𝑒 ± / 𝐿 𝛾, IC ratios are: 21.0 for the baseline injection model, 14.3for inj1 , 124.7 for inj2 , 23.4 for inj3 , and 21.1 for inj4 . See also Sec. 6.1 fordetails. and a factor of a few for the other scenarios. As in all previous cases,we summarize the minimum luminosities for detection in Table 4.We note that in this section we have examined the degradationof the sensitivity due to uncertainties in the FBs spectrum only. How-ever, in Appendix B we also consider the case “FB max” togetherwith mismodeling of the GDE model. Those results confirm thatuncertainties in the FBs model could be the single most challengingastrophysical background component for analyses of extended 𝛾 -rayemission with CTA in the GC region. Minimum 𝑓 𝑒 ± for detection [%] Baseline Inj1 Inj2 Inj3 Inj4 FB min , perfect GDE.10 .
5% 2 .
9% 158 .
4% 24 .
3% 8 . min , mismodeling of the GDE.14 .
5% 3 .
8% 163 .
4% 25 .
3% 10 . max , perfect GDE.57 .
5% 41 .
3% 259 .
4% 55 .
0% 58 . max , mismodeling of the GDE.72 .
9% 51 .
8% 326 .
7% 70 .
4% 74 . Table 5.
Minimum MSPs 𝑒 ± injection efficiency ( 𝑓 𝑒 ± ) required for a 5 𝜎 significance detection with CTA. The computation of these efficiencies usethe implied 𝑒 ± luminosities (see Sec. 6.1) based on the minimum luminosi-ties reported in Table 4, the measured Fermi
GeV excess 𝛾 -ray luminosity,and Eq. 2. 𝑒 ± injection efficiency ( 𝑓 𝑒 ± ) for detectionwith CTA In the previous section, we investigated the ability of CTA to charac-terize the TeV-scale IC 𝛾 rays produced by an unresolved populationof MSPs in the Galactic bulge region. In this section, we evaluatewhether the signals that can be detected by CTA are physicallypossible.Using the relations presented in Eqs. 2, 5, and assuming thatthe prompt 𝛾 -ray emission from the Galactic bulge population ofMSPs is fully responsible for the GCE, we can obtain the minimum Note that analyses of the GCE (Lacroix et al. 2016) did not detect the ICsignature from the putative population of MSPs in the Galactic bulge. Thismight be because the prompt 𝛾 rays are much more prominent that the IC 𝛾 rays at GeV scale energies.MNRAS000
GeV excess 𝛾 -ray luminosity,and Eq. 2. 𝑒 ± injection efficiency ( 𝑓 𝑒 ± ) for detectionwith CTA In the previous section, we investigated the ability of CTA to charac-terize the TeV-scale IC 𝛾 rays produced by an unresolved populationof MSPs in the Galactic bulge region. In this section, we evaluatewhether the signals that can be detected by CTA are physicallypossible.Using the relations presented in Eqs. 2, 5, and assuming thatthe prompt 𝛾 -ray emission from the Galactic bulge population ofMSPs is fully responsible for the GCE, we can obtain the minimum Note that analyses of the GCE (Lacroix et al. 2016) did not detect the ICsignature from the putative population of MSPs in the Galactic bulge. Thismight be because the prompt 𝛾 rays are much more prominent that the IC 𝛾 rays at GeV scale energies.MNRAS000 , 1–22 (2021) Macias et al. T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Baseline5 significance T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj15 significance T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj25 significance T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj35 significance Injected Luminosity [erg/s] T e s t s t a t i s t i c
5 significance Injected Luminosity [erg/s] R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj45 significance FB min , mismodeling of the GDE Figure 7.
Same as Fig. 4, except that here we assume mismodeling of the Galactic diffuse emission model and the FB min model. Note that mismodeling of theGDE is mocked up by generating the data with “GDE model 1” and then fitting the data with “GDE model 2”. MNRAS , 1–22 (2021)
TA sensitivity to the putative millisecond pulsar population responsible for the GeV excess E d N / d E [ T e V c m s s r ] injected baseline signalextracted signal injected signal 1extracted signal E d N / d E [ T e V c m s s r ] injected signal 2extracted signal Energy [TeV] injected signal 3extracted signal Energy [TeV] E d N / d E [ T e V c m s s r ] injected signal 4extracted signal FB min , mismodeling of the GDE Figure 8.
Same as Fig. 5, except that here we assume mismodeling of the Galactic diffuse emission model. See also the green dotted line in Fig. 7 for thenecessary IC luminosity for a 5 𝜎 significance detection of the signal. efficiencies 𝑓 𝑒 ± ’s for a CTA detection of the MSPs’ IC signal. Inparticular, we can convert the minimum IC luminosities (shown inTable 4) to the corresponding minimum 𝐿 𝑒 ± ’s using Eq. 5, and thenevaluate this value in Eq. 2, along with the inferred nominal GCEluminosity [ 𝐿 𝛾, prompt = . × erg/s obtained in e.g., Maciaset al. (2019)]. However, the threshold IC luminosities presentedin Table 4 are estimated in the energy range 16 GeV to 158 TeV,while the GCE luminosity in Macias et al. (2019) was computedfor 𝐸 𝛾 (cid:38)
700 MeV. So, in order for us to connect the threshold ICluminosities with the 𝑒 ± injection luminosities that were used in ourgalprop runs, we need to extend the 𝑒 ± luminosity calculation to700 MeV . This is a good approximation since the injected MSPs 𝑒 ± ’s can reach In summary, our luminosity computations assume; 𝐸 𝑒 ± ≥ 𝐸𝛾 ≥
700 MeV, a distance from the Sun to the GC of 8 . ◦ × ◦ around the GC. It is useful to compare the fractional luminosities( 𝐿 𝑒 ± / 𝐿 𝛾, IC ) predicted by galprop—estimated by calculating theluminosity in the galprop IC maps, and the 𝐿 𝑒 ± used as input ingalprop—so as to have a better understanding of the impact ofdiffusion and energy losses. We obtain that 𝐿 𝑒 ± / 𝐿 𝛾, IC is 21.0 forthe baseline injection model, 14.3 for inj1 , 124.7 for inj2 , 23.4 for inj3 , and 21.1 for inj4 (see also Table 2). The very large luminosity roughly the same minimum energies as the 𝛾 -rays. Also, note that by com-paring the 𝑒 ± luminosities included in galprop— before propagation —withthe threshold IC luminosities, we automatically account for the effects ofpropagation and other energy losses (like synchrotron) for the MSPs 𝑒 ± .MNRAS000
700 MeV, a distance from the Sun to the GC of 8 . ◦ × ◦ around the GC. It is useful to compare the fractional luminosities( 𝐿 𝑒 ± / 𝐿 𝛾, IC ) predicted by galprop—estimated by calculating theluminosity in the galprop IC maps, and the 𝐿 𝑒 ± used as input ingalprop—so as to have a better understanding of the impact ofdiffusion and energy losses. We obtain that 𝐿 𝑒 ± / 𝐿 𝛾, IC is 21.0 forthe baseline injection model, 14.3 for inj1 , 124.7 for inj2 , 23.4 for inj3 , and 21.1 for inj4 (see also Table 2). The very large luminosity roughly the same minimum energies as the 𝛾 -rays. Also, note that by com-paring the 𝑒 ± luminosities included in galprop— before propagation —withthe threshold IC luminosities, we automatically account for the effects ofpropagation and other energy losses (like synchrotron) for the MSPs 𝑒 ± .MNRAS000 , 1–22 (2021) Macias et al. Injected Luminosity [erg/s] ) T e s t s t a t i s t i c Detection threshold (5 ) Injected Luminosity [erg/s] R e c o v e r e d L u m i n o s i t y [ e r g / s ] injected signalrecovered MSP signalRecovered DM signal FB min , mismodeling of the GDE Figure 9.
Same as Fig. 6, except that here we assume mismodeling of the Galactic diffuse emission. fraction obtained for inj2 is explained by the fact that this injectionspectrum is very soft.Using the prescription described above, we are now able tocompute the threshold 𝑓 𝑒 ± values for the cases considered in ourstudy. We show the results of this calculation in Table 5. Dependingon assumptions about the astrophysical background componentsand the 𝑒 ± injection model, we obtain threshold 𝑒 ± efficiencies inthe range 𝑓 𝑒 ± ≈ . − . inj2 —which is the softest 𝑒 ± injection spectra considered in our sample. Indeed, we obtainthat the 𝑒 ± luminosity needed for CTA to detect a soft spectrumlike inj2 would exceed the total budget of the MSPs spin-downenergy. This means that, if the most pessimistic GDE mismodelingscenario (“FB max” and mismodeling of the GDE) considered hereis realized in nature, CTA will only be able to reliably detect theGalactic bulge population of MSPs if the overall efficiency of thispopulation satisfies 𝑓 𝑒 ± (cid:38) .
8% (see the last row of Table 5).Notice that CTA might still be suited to detect this signal withpercentage-level 𝑓 𝑒 ± ’s under some specific conditions consideredin Table 5.Interestingly, the recent study of Song et al. (2021) obtained 𝑓 𝑒 ± ≈ 𝑓 𝑒 ± ≈ ( . − . ) %, for MSPs in M15.However, very likely these strong constraints cannot be directlyextrapolated to other systems containing MSPs. This is becausethe overall apparent efficiency in globular clusters can be stronglydecreased by rapid winds from Red Clump giants in globular clus-ters (Bednarek et al. 2016). Note that winds can advect cosmic-ray 𝑒 ± out of the globular clusters systems before they can radiate.Additional clues about the 𝑓 𝑒 ± in systems containing MSPpopulations have been obtained in the recent analysis by Sudohet al. (2020). Based on the break-down of a correlation between far-infrared and radio luminosities in SFGs, Sudoh et al. (2020) positedthat radio emission from MSPs could account for a large fractionof the radio luminosity observed in systems with high stellar massand low star formation rate. This led the authors to conclude that the MSP populations in their sample of SFGs could have 𝑓 𝑒 ± inexcess of 90%. They also noted, however, that several observationaland theoretical uncertainties could lower their inferred efficiency to 𝑓 𝑒 ± ≈ 𝑓 𝑒 ± for MSPsare currently not very well constrained. In case that CTA makes anactual observation of the MSP IC signal in one of the scenariosdisfavoured by our analysis, it could still be possible to reconcilesuch an observation with the MSP emission models. In particular,if the 𝛾 -ray emission from the Galactic bulge MSPs magnetosphereis beamed, only some fraction of their prompt 𝛾 -ray luminosity canbe observed from the Earth. This could decrease our inferred 𝑓 𝑒 ± by factor of a few (Sudoh et al. 2020).In conclusion, we have demonstrated that—even under theassumption of high background uncertainties—if the 𝑒 ± injectionspectra is harder than our inj2 model (slope of 𝑑𝑁 / 𝑑𝐸 ∝ 𝐸 − . ),CTA has the potential to robustly discover the IC signal produced bya new population of MSPs in the GC. However, given observationaland theoretical uncertainties on the 𝑓 𝑒 ± parameter, a detection of asignal described by our inj2 model could still be possible. In this work, we have utilized a spatio-spectral template regressionmethod with simulated CTA data from the GC region. In particular,we have run the signal recovery tests using a bin-by-bin analysis,which has been utilized with great success in analyses of
Fermi -LAT data [e.g., (Ackermann et al. 2015, 2017)]. This methodologyallows us to reduce the impact of potential biases introduced byassumptions about the spectrum of the MSPs and/or DM templates.Using the aforementioned method, we tested whether CTAcould disentangle the IC signal produced by a unresolved popula-tion of MSPs in the Galactic bulge from the one produced by aspherical distribution of DM (or a population of pulsars following aspherical distribution). Given that the 𝑒 ± sources would follow ei- MNRAS , 1–22 (2021)
TA sensitivity to the putative millisecond pulsar population responsible for the GeV excess T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Baseline5 significance T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj15 significance T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj25 significance T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj35 significance Injected Luminosity [erg/s] T e s t s t a t i s t i c
5 significance Injected Luminosity [erg/s] R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj45 significance FB max , Perfect GDE Figure 10.
Same as Fig. 4, except that here we assume the FB max model.MNRAS000
Same as Fig. 4, except that here we assume the FB max model.MNRAS000 , 1–22 (2021) Macias et al. E d N / d E [ T e V c m s s r ] injected baseline signalextracted signal injected signal 1extracted signal E d N / d E [ T e V c m s s r ] injected signal 2extracted signal Energy [TeV] injected signal 3extracted signal Energy [TeV] E d N / d E [ T e V c m s s r ] injected signal 4extracted signal FB max , Perfect GDE Figure 11.
Same as Fig. 5, except that here we assume the FB max model. ther of these two distributions , the IC maps produced by such 𝑒 ± were predicted to have discernible spatial differences in Song et al.(2019). We have injected signals of MSP IC emission with varyingstrengths, and then run the signal recovery pipeline including both,an MSP IC template and a DM IC template in the fit. As a result,we have found that CTA has the capability of robustly disentanglingthese two sources, even in the presence of Galactic diffuse emissionmismodeling. Overall, we have demonstrated that if CTA discoversa diffuse IC signal under the conditions considered in Table 4, thespatial morphology of the IC signal will reveal whether it is relatedto stellar mass or a spherically-symmetric source distribution as, forinstance, would be expected for DM. The reader is referred to Sec. 2.2 for details on the spatial morphologiesconsidered in this work.
We note that TeV-scale 𝑒 ± pairs could potentially lose mostof their energy very close to parent MSPs’ magnetospheres. If thenumber of MSPs responsible for the GCE is relatively small, thenthe predicted IC templates could exhibit a clustering-of-photonseffect (Acharyya et al. 2021), which could facilitate the detec-tion of the MSP population in the Galactic bulge. Very promis-ing methodologies to study these effects have been explored in theliterature and include: the non-Poissonian template fitting proce-dure Lee et al. (2016); Leane & Slatyer (2019); Chang et al. (2020);Leane & Slatyer (2020a,b); Buschmann et al. (2020), wavelet tech-niques Bartels et al. (2016); Balaji et al. (2018); Zhong et al. (2020),deep learning methods (Caron et al. 2018; List et al. 2020), radiodetection (Calore et al. 2016; Macquart & Kanekar 2015; Rajwadeet al. 2017; Hyman et al. 2019), and X-ray detection (Berteaud et al.2020) of point sources responsible for the Galactic bulge emission.Importantly, recent population synthesis models of MSPs (Ploeg MNRAS , 1–22 (2021)
TA sensitivity to the putative millisecond pulsar population responsible for the GeV excess Injected Luminosity [erg/s] ) T e s t s t a t i s t i c Detection threshold (5 ) Injected Luminosity [erg/s] R e c o v e r e d L u m i n o s i t y [ e r g / s ] injected signalrecovered MSP signalRecovered DM signal FB max , Perfect GDE Figure 12.
Same as Fig. 6, except that here we assume the FB max model. et al. 2020) predict anywhere between 20 and 50 thousand MSPs inthe Galactic bulge. We can use the selected spatial resolution of oursimulations, and Eqs. A1; A2; and A3, to obtain an estimate of thenumber of MSPs in each spatial 3D bin of our galprop simulations.We find that every bin of size 200 × ×
100 pc is expected tocontain ∼
300 MSPs in the center of the boxy bulge and ∼
30 MSPsat 3 kpc along the long-axis of the Galactic bar. Furthermore, forthe nuclear bulge region we estimate ∼ We have estimated the impact of various diffuse astrophysical com-ponents on the sensitivity to the signal. Using Monte Carlo simu-lations, we recreated a scenario in which the GDE is mismodeled,obtaining that if this scenario is realized in actual data, it will de-teriorate the CTA sensitivity to the MSPs’ IC signal luminosity atthe (cid:46)
40% level (depending on the characteristics of the injectionspectrum and assumptions about the FBs model). One possible ex-planation as to why this effect is not larger in our analyses, is thatthe intensity of the predicted GDE spectrum rapidly falls off withenergy, becoming comparable to that of the expected MSP IC signalat TeV-scale 𝛾 -ray energies. We note that this is in stark contrastto analyses of Fermi -LAT data from the GC region in which, forexample, the
Fermi
GeV excess signal has an intensity that is just asmall fraction of the GDE emission [e.g., Abazajian et al. (2020)].On the other hand, we observed a drastic deterioration (upto one order of magnitude) of the IC flux sensitivity when weswitched from the “FB min” to the “FB max” model (see Table 4).The latter was proposed in Rinchiuso et al. (2021) (see the left panelof Fig. 1 in that article), as a way of conservatively accounting forthe maximum 𝛾 -ray intensity that the FBs can take at the CTAenergy range. This was accomplished in that work by ensuring thatforthcoming measurements of diffuse 𝛾 rays from this sky regioncannot overshoot current H.E.S.S. measurements of the Galacticridge region (inner ≈
200 pc of the GC). As shown in Fig. 3, the“FB max” model is the dominant astrophysical 𝛾 -ray component for energies greater than ≈
30 GeV. This explains why the FBs modelproduces the strongest impact on the sensitivity to the MSPs’ ICsignal.Another important source of uncertainty corresponds to thespatial morphology of the FBs. The FBs template assumed in ourwork is a inpainted version of the spatial map in Ackermann et al.(2017), which in turn was constructed using a spectral componentanalysis using residual 𝛾 -ray data in the 1 −
10 GeV range. However,the exact details of the FBs map might be susceptible to the energyrange that is included in the spectral component analysis of Acker-mann et al. (2017). In a future study, we will address the impact ofspatial uncertainties in the FBs model, on the CTA sensitivity to anMSP IC signal in the GC. 𝑒 ± injection spectrum on thesensitivity to the IC signal As discussed in the previous section, one of the most importantfactors determining the characteristics of the expected IC signal,corresponds to the 𝑒 ± injection spectrum [see also (Song et al.2019)]. In our analysis we considered five different 𝑒 ± injectionspectra, finding that these produce significantly different results.For example, the injection model Inj1 (see Table 4) has a detectionthreshold luminosity that is a factor of ∼ Inj2 in both “FB min” cases under consideration.The most likely explanation for this is that signals that mirror thespectral shape of the GDE—see the left-middle panel of Fig. 5in comparison to the left panel of Fig. 3—are more difficult todisentangle than the signals that have a distinct spectral shape.Future advances in our understanding of the MSPs’ 𝑒 ± injectionspectrum will help reduce the impact of this component on the sensi-tivity. On the theory side, future global magnetosphere simulationsthat include non-dipolar fields could provide a clearer picture of the 𝑒 ± injection mechanisms and predicted spectrum (Harding 2021).On the phenomenological side, multiwavelength measurements andmodeling of unresolved MSPs in globular clusters (Ndiyavala et al.2018; Song et al. 2021) and SFGs (Sudoh et al. 2020), could revealthe characteristics of the MSPs 𝑒 ± injection spectra. MNRAS000
10 GeV range. However,the exact details of the FBs map might be susceptible to the energyrange that is included in the spectral component analysis of Acker-mann et al. (2017). In a future study, we will address the impact ofspatial uncertainties in the FBs model, on the CTA sensitivity to anMSP IC signal in the GC. 𝑒 ± injection spectrum on thesensitivity to the IC signal As discussed in the previous section, one of the most importantfactors determining the characteristics of the expected IC signal,corresponds to the 𝑒 ± injection spectrum [see also (Song et al.2019)]. In our analysis we considered five different 𝑒 ± injectionspectra, finding that these produce significantly different results.For example, the injection model Inj1 (see Table 4) has a detectionthreshold luminosity that is a factor of ∼ Inj2 in both “FB min” cases under consideration.The most likely explanation for this is that signals that mirror thespectral shape of the GDE—see the left-middle panel of Fig. 5in comparison to the left panel of Fig. 3—are more difficult todisentangle than the signals that have a distinct spectral shape.Future advances in our understanding of the MSPs’ 𝑒 ± injectionspectrum will help reduce the impact of this component on the sensi-tivity. On the theory side, future global magnetosphere simulationsthat include non-dipolar fields could provide a clearer picture of the 𝑒 ± injection mechanisms and predicted spectrum (Harding 2021).On the phenomenological side, multiwavelength measurements andmodeling of unresolved MSPs in globular clusters (Ndiyavala et al.2018; Song et al. 2021) and SFGs (Sudoh et al. 2020), could revealthe characteristics of the MSPs 𝑒 ± injection spectra. MNRAS000 , 1–22 (2021) Macias et al.
ACKNOWLEDGEMENTS
O.M. acknowledges support by JSPS KAKENHI Grant NumberJP20K14463. The work of D.S. is supported by the U.S. Depart-ment of Energy under the award number DE-SC0020262. S.A. issupported by JSPS/MEXT KAKENHI Grant numbers JP17H04836and JP18H04340. The work of S.H. is supported by the U.S. De-partment of Energy under the award number DE-SC0020262 andNSF Grant numbers AST-1908960 and PHY-1914409. This workwas supported by World Premier International Research Center Ini-tiative (WPI Initiative), MEXT, Japan. R.M.C. acknowledges sup-port from the Australian Government through the Australian Re-search Council for grant DP190101258 shared with Prof. MarkKrumholz at the ANU. D.M.N. acknowledges support from NASAunder award Number 80NSSC19K0589. The authors acknowledgeAdvanced Research Computing at Virginia Tech for providing com-putational resources and technical support that have contributed tothe results reported within this paper.
REFERENCES
Abazajian K. N., 2011, JCAP, 1103, 010Abazajian K. N., Kaplinghat M., 2012, Phys. Rev., D86, 083511Abazajian K. N., Horiuchi S., Kaplinghat M., Keeley R. E., Macias O., 2020,Phys. Rev. D, 102, 043012Abdo A. A., et al., 2013, Astrophys. J. Suppl., 208, 17Abdollahi S., et al., 2020, Astrophys. J. Suppl., 247, 33Abeysekara A., et al., 2017, Science, 358, 911Abramowski A., et al., 2013, Astron. Astrophys., 551, A26Acciari V. A., et al., 2019, Mon. Not. Roy. Astron. Soc., 484, 2876Acero F., et al., 2016, Astrophys. J. Suppl., 223, 26Acharya B., et al., 2018, Science with the Cherenkov Telescope Array. WSP( arXiv:1709.07997 ), doi:10.1142/10986Acharyya A., et al., 2021, JCAP, 01, 057Ackermann M., et al., 2012, ApJ, 750, 3Ackermann M., et al., 2014, Astrophys. J., 793, 64Ackermann M., et al., 2015, Phys. Rev. Lett., 115, 231301Ackermann M., et al., 2017, Astrophys. J., 840, 43Aguilar M., et al., 2013, Phys. Rev. Lett., 110, 141102Aharonian F. A., Atoyan A. M., Voelk H. J., 1995, A&A , 294, L41Aharonian F., et al., 2006, Nature, 439, 695Ajello M., et al., 2016, Astrophys. J., 819, 44Ajello M., et al., 2017, Astrophys. J. Suppl., 232, 18Atoyan A. M., Aharonian F. A., Völk H. J., 1995, Phys. Rev. D, 52, 3265Atri P., et al., 2019, Mon. Not. Roy. Astron. Soc., 489, 3116Balaji B., Cholis I., Fox P. J., McDermott S. D., 2018, Phys. Rev., D98,043009Bartels R., Krishnamurthy S., Weniger C., 2016, Phys. Rev. Lett., 116,051102Bartels R., Storm E., Weniger C., Calore F., 2018, Nat. Astron., 2, 819Beck R., 2001, Space Science Reviews, 99, 243Bednarek W., Sitarek J., Sobczak T., 2016, Mon. Not. Roy. Astron. Soc.,458, 1083Berteaud J., Calore F., Clavel M., Serpico P. D., Dubus G., Petrucci P.-O.,2020, arXiv e-prints, p. arXiv:2012.03580Buschmann M., Rodd N. L., Safdi B. R., Chang L. J., Mishra-Sharma S.,Lisanti M., Macias O., 2020, Phys. Rev. D, 102, 023023Calore F., Cholis I., Weniger C., 2015, JCAP, 03, 038Calore F., Di Mauro M., Donato F., Hessels J. W. T., Weniger C., 2016,Astrophys. J., 827, 143Caron S., Gómez-Vargas G. A., Hendriks L., Ruiz de Austri R., 2018, JCAP,05, 058Chang L. J., Mishra-Sharma S., Lisanti M., Buschmann M., Rodd N. L.,Safdi B. R., 2020, Phys. Rev. D, 101, 023014Coleman B., Paterson D., Gordon C., Macias O., Ploeg H., 2020, Mon. Not.Roy. Astron. Soc., 495, 3350 Coronado-Blázquez J., Doro M., Sánchez-Conde M. A., Aguirre-SantaellaA., 2021, arXiv e-prints, p. arXiv:2101.10003Crocker R. M., Jones D., Melia F., Ott J., Protheroe R. J., 2010, Nature, 468,65Crocker R. M., Bicknell G. V., Taylor A. M., Carretti E., 2015, Astrophys.J., 808, 107Daylan T., Finkbeiner D. P., Hooper D., Linden T., Portillo S. K. N., RoddN. L., Slatyer T. R., 2016, Phys. Dark Univ., 12, 1Deil C., et al., 2017, in 35th International Cosmic Ray Conference(ICRC2017). p. 766 ( arXiv:1709.01751 )Delahaye T., Lavalle J., Lineros R., Donato F., Fornengo N., 2010, A&A ,524, A51Di Mauro M., 2021, arXiv e-prints, p. arXiv:2101.04694Di Mauro M., Manconi S., Donato F., 2019, Phys. Rev. D, 100, 123015Erber T., 1966, Reviews of Modern Physics, 38, 626Freudenreich H. T., 1998, ApJ, 492, 495Garzon F., Lopez-Corredoira M., Hammersley P., Mahoney T. J., Calbet X.,Beckman J. E., 1997, Astrophys. J., 491, L31Gérard L., 2016, PoS, ICRC2015, 725Gonthier P. L., Harding A. K., Ferrara E. C., Frederick S. E., Mohr V. E.,Koh Y.-M., 2018, Astrophys. J., 863, 199Gordon C., Macias O., 2013, Phys. Rev., D88, 083521Guépin C., Cerutti B., Kotera K., 2020, Astron. Astrophys., 635, A138Hammersley P. L., Garzon F., Mahoney T., Lopez-Corredoira M., Torres M.A. P., 2000, Mon. Not. Roy. Astron. Soc., 317, L45Harding A. K., 2021, arXiv e-prints, p. arXiv:2101.05751Hassan T., et al., 2015, in 34th International Cosmic Ray Conference(ICRC2015). p. 971 ( arXiv:1508.06075 )Heiles C., 1995, in Ferrara A., McKee C. F., Heiles C., Shapiro P. R., eds,Astronomical Society of the Pacific Conference Series Vol. 80, ThePhysics of the Interstellar Medium and Intergalactic Medium. p. 507Herold L., Malyshev D., 2019, Astron. Astrophys., 625, A110Hobbs G., Lorimer D., Lyne A., Kramer M., 2005, Mon. Not. Roy. Astron.Soc., 360, 974Hooper D., Goodenough L., 2011, Phys. Lett., B697, 412Hooper D., Blasi P., Serpico P. D., 2009, JCAP, 01, 025Hooper D., Cholis I., Linden T., Fang K., 2017, Phys. Rev. D, 96, 103013Horiuchi S., Kaplinghat M., Kwa A., 2016, JCAP, 11, 053Hyman S. D., Frail D. A., Deneva J. S., Kassim N. E., McLaughlin M. A.,Kooi J. E., Ray P. S., Polisensky E. J., 2019, ApJ, 876, 20Ishiwata K., Macias O., Ando S., Arimoto M., 2020, JCAP, 01, 003Jóhannesson G. l., Porter T. A., Moskalenko I. V., 2018, Astrophys. J., 856,45Jóhannesson G., Porter T. A., Moskalenko I. V., 2019, Astrophys. J., 879,91Lacroix T., Macias O., Gordon C., Panci P., Bœhm C., Silk J., 2016, Phys.Rev. D, 93, 103004Launhardt R., Zylka R., Mezger P., 2002, Astron. Astrophys., 384, 112Leane R. K., Slatyer T. R., 2019, Phys. Rev. Lett., 123, 241101Leane R. K., Slatyer T. R., 2020a, Phys. Rev. D, 102, 063019Leane R. K., Slatyer T. R., 2020b, Phys. Rev. Lett., 125, 121105Lee S. K., Lisanti M., Safdi B. R., Slatyer T. R., Xue W., 2016, Phys. Rev.Lett., 116, 051103Linden T., Rodd N. L., Safdi B. R., Slatyer T. R., 2016, Phys. Rev. D, 94,103013List F., Rodd N. L., Lewis G. F., Bhat I., 2020, Phys. Rev. Lett., 125, 241102Macias O., Gordon C., 2014, Phys. Rev. D, 89, 063515Macias O., Gordon C., Crocker R. M., Coleman B., Paterson D., HoriuchiS., Pohl M., 2018, Nat. Astron., 2, 387Macias O., Horiuchi S., Kaplinghat M., Gordon C., Crocker R. M., NatafD. M., 2019, JCAP, 09, 042Macquart J.-P., Kanekar N., 2015, ApJ, 805, 172Michel F. C., 1991, Theory of neutron star magnetospheresNdiyavala-Davids H., Venter C., Kopp A., Backes M., 2020, Mon. Not. Roy.Astron. Soc., 500, 4827Ndiyavala H., Krüger P. P., Venter C., 2018, Mon. Not. Roy. Astron. Soc.,473, 897Nelder J., Mead R., 1965, Comput. J., 7, 308 MNRAS , 1–22 (2021)
TA sensitivity to the putative millisecond pulsar population responsible for the GeV excess Nigro C., et al., 2019, A&A , 625, A10Pfahl E., Rappaport S., Podsiadlowski P., 2002, Astrophys. J., 573, 283Ploeg H., Gordon C., Crocker R., Macias O., 2017, JCAP, 1708, 015Ploeg H., Gordon C., Crocker R., Macias O., 2020, JCAP, 12, 035Podsiadlowski P., Pfahl E., Rappaport S., 2005, in Rasio F. A., Stairs I. H.,eds, Astronomical Society of the Pacific Conference Series Vol. 328,Binary Radio Pulsars. p. 327Pohl M., Englmaier P., Bissantz N., 2008, Astrophys.~J., 677, 283Porter T. A., Moskalenko I. V., Strong A. W., 2006, Astrophys. J. Lett., 648,L29Porter T. A., Jóhannesson G., Moskalenko I. V., 2017, Astrophys. J., 846,67Profumo S., Reynoso-Cordova J., Kaaz N., Silverman M., 2018, Phys. Rev.D, 97, 123008Rajwade K. M., Lorimer D. R., Anderson L. D., 2017, MNRAS, 471, 730Rinchiuso L., Macias O., Moulin E., Rodd N. L., Slatyer T. R., 2021, Phys.Rev. D, 103, 023011Robitaille T. P., Churchwell E., Benjamin R. A., Whitney B. A., Wood K.,Babler B. L., Meade M. R., 2012, A&A, 545, A39Silverwood H., Weniger C., Scott P., Bertone G., 2015, JCAP, 03, 055Song D., Macias O., Horiuchi S., 2019, Phys. Rev., D99, 123020Song D., Macias O., Horiuchi S., Crocker R. M., Nataf D. M., 2021, arXive-prints, p. arXiv:2102.00061Storm E., Weniger C., Calore F., 2017, JCAP, 08, 022Strong A. W., Moskalenko I. V., Reimer O., 2000, Astrophys. J., 537, 763Strong A. W., Moskalenko I. V., Ptuskin V. S., 2007, Ann. Rev. Nucl. Part.Sci., 57, 285Sturrock P. A., 1970, Nature , 227, 465Sturrock P. A., 1971, ApJ, 164, 529Su M., Finkbeiner D. P., 2012, ApJ, 753, 61Su M., Slatyer T. R., Finkbeiner D. P., 2010, ApJ, 724, 1044Sudoh T., Linden T., Beacom J. F., 2020, arXiv e-prints, p. arXiv:2005.08982Viana A., et al., 2020, PoS, ICRC2019, 817Yuan Q., Ioka K., 2015, Astrophys. J., 802, 124Zhong Y.-M., McDermott S. D., Cholis I., Fox P. J., 2020, Phys. Rev. Lett.,124, 231103
APPENDIX A: STELLAR DENSITY PROFILESA1 Nuclear bulge
Motivated by recent reanalysis of the GCE [e.g., Macias et al.(2019)], we assumed here that the MSP population that accounts forthe GCE follows the distribution of stars in the nuclear bulge and theboxy bulge. For the nuclear bulge, we used the three-dimensionaldensity function given in Launhardt et al. (2002) [see also Songet al. (2019)]. As mentioned in Sec. 3, the nuclear bulge consists ofthe nuclear stellar disk (NSD) and the nuclear stellar cluster (NSC).The NSD density function is given by 𝜌 NSD ( 𝑥, 𝑦, 𝑧 ) = 𝜌 (cid:16) 𝑟𝑟 𝑑 (cid:17) − . exp (cid:16) − | 𝑧 | 𝑧 𝑑 (cid:17) , for 𝑟 / pc < 𝜌 (cid:16) 𝑟𝑟 𝑑 (cid:17) − . exp (cid:16) − | 𝑧 | 𝑧 𝑑 (cid:17) , for 120 ≤ 𝑟 / pc < 𝜌 (cid:16) 𝑟𝑟 𝑑 (cid:17) − exp (cid:16) − | 𝑧 | 𝑧 𝑑 (cid:17) , for 𝑟 / pc ≥ 𝑟 = 𝑥 + 𝑦 + 𝑧 , 𝑟 𝑑 = 𝑧 𝑑 =
45 pc and 𝜌 = (cid:12) pc − . The total stellar mass of the NSD in the inner 120 pc is setto 8 × M (cid:12) , which then defines 𝜌 , and 𝜌 through continuityconditions. In addition, the NSC is given by 𝜌 𝑁 𝑆𝐶 ( 𝑥, 𝑦, 𝑧 ) = 𝜌 (cid:20) + (cid:16) 𝑟𝑟 𝑐 (cid:17) (cid:21) − , for 𝑟 / pc < 𝜌 (cid:20) + (cid:16) 𝑟𝑟 𝑐 (cid:17) (cid:21) − , for 6 < 𝑟 / pc ≤ , for 𝑟 / pc >
200 (A2)where 𝑟 𝑐 = .
22 pc, and 𝜌 = . × M (cid:12) . The density 𝜌 isobtained by imposing continuity at 𝑟 = A2 The boxy bulge
The boxy bulge stellar density function was derived by Freudenreich(1998). We used the Model S in that reference, which is given by 𝜌 bar ( 𝑅, 𝜙, 𝑧 ) ∝ sech ( 𝑅 𝑠 ) × (cid:40) , 𝑅 ≤ 𝑅 end 𝑒 −[( 𝑅 − 𝑅 end )/ ℎ end ] , 𝑅 > 𝑅 end (A3)where 𝑅 end = .
128 kpc, ℎ end = .
461 kpc, and 𝑅 𝑠 is defined as 𝑅 𝐶 ⊥ ⊥ = (cid:0)(cid:12)(cid:12) 𝑋 (cid:48) (cid:12)(cid:12) / 𝑎 𝑥 (cid:1) 𝐶 ⊥ + (cid:0)(cid:12)(cid:12) 𝑌 (cid:48) (cid:12)(cid:12) / 𝑎 𝑦 (cid:1) 𝐶 ⊥ 𝑅 𝐶 (cid:107) 𝑠 = 𝑅 𝐶 | ⊥ + (cid:0)(cid:12)(cid:12) 𝑍 (cid:48) (cid:12)(cid:12) / 𝑎 𝑧 (cid:1) 𝐶 (cid:107) , (A4)with 𝑎 𝑥 = .
696 kpc, 𝑎 𝑦 = . 𝑎 𝑧 = . 𝐶 (cid:107) = . 𝐶 ⊥ = . ( 𝑋 (cid:48) , 𝑌 (cid:48) , 𝑍 (cid:48) ) correspond to a Cartesian coordinate system in the boxy bulge frameof reference (see Song et al. 2019 for more details). A3 Spherically symmetric template
We used a generalized NFW density profile of the form 𝜌 NFW ( 𝑟 ) = 𝜌 (cid:16) 𝑟𝑅 (cid:12) (cid:17) 𝛾 (cid:16) + 𝑟 / 𝑅 𝑠 + 𝑅 𝑠 / 𝑅 (cid:12) (cid:17) ( − 𝛾 ) (A5)where 𝛼 = 𝛽 = 𝑅 (cid:12) =
20 kpc, and 𝛾 = .
2. Previousstudies of the GCE (Abazajian & Kaplinghat 2012; Gordon & Ma-cias 2013; Daylan et al. 2016; Calore et al. 2016) showed that thesquare of this density profile ( 𝜌 ) was a good fit to the signalmorphology. Although more recent analyses (Macias et al. 2018;Bartels et al. 2018; Macias et al. 2019; Abazajian et al. 2020) usingimproved background models have shown that stellar mass profilesof the Galactic bulge give a much better fit to the data, it is stillinteresting to consider this morphology because it would test theability of CTA to detect the IC signature from DM self-annihilation. APPENDIX B: IMPACT OF THE FERMI BUBBLESMODEL IN THE PRESENCE OF GALACTIC DIFFUSEEMISSION MISMODELING
Here, we present the results of the systematic uncertainties analysisin the case where the FBs component have the maximum possibleintensity (FB max scenario) and the GDE components are mismod-eled. We show the signal recovering tests in Fig. B1, the minimumflux for detection in Fig. B2, and the MSPs vs. DM degeneracy testsin Fig. B3. We also present the threshold IC luminosities and 𝑓 𝑒 ± inTable 4, and Table 5, respectively.We note that this constitutes the most challenging scenario for MNRAS000
Here, we present the results of the systematic uncertainties analysisin the case where the FBs component have the maximum possibleintensity (FB max scenario) and the GDE components are mismod-eled. We show the signal recovering tests in Fig. B1, the minimumflux for detection in Fig. B2, and the MSPs vs. DM degeneracy testsin Fig. B3. We also present the threshold IC luminosities and 𝑓 𝑒 ± inTable 4, and Table 5, respectively.We note that this constitutes the most challenging scenario for MNRAS000 , 1–22 (2021) Macias et al. a reliable CTA detection of the MSPs IC signal. As was shown in theleft panel of Fig. 3, our FB max model outshines all other componentsat the CTA energy range. We remind the reader that this model wasconstructed by imposing that the FBs model does not overshootthe H.E.S.S. measurements from the Galactic ridge region. Toughthe exact spatial morphology and spectra of the low-latitude part ofthe FBs are very uncertain, our very conservative modeling of thiscomponent allows to show that CTA has the capability to reliablydiscover the IC signature from an unresolved population of MSPsin the Galactic bulge.
This paper has been typeset from a TEX/L A TEX file prepared by the author. MNRAS , 1–22 (2021)
TA sensitivity to the putative millisecond pulsar population responsible for the GeV excess T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Baseline5 significance T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj15 significance T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj25 significance T e s t s t a t i s t i c
5 significance R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj35 significance Injected Luminosity [erg/s] T e s t s t a t i s t i c
5 significance Injected Luminosity [erg/s] R e c o v e r e d L u m i n o s i t y [ e r g / s ] Inj45 significance FB max , mismodeling of the GDE Figure B1.
Same as Fig. 4, except that here we assume mismodeling of the Galactic diffuse emission model and the FB max model.MNRAS000
Same as Fig. 4, except that here we assume mismodeling of the Galactic diffuse emission model and the FB max model.MNRAS000 , 1–22 (2021) Macias et al. E d N / d E [ T e V c m s s r ] injected baseline signalextracted signal injected signal 1extracted signal E d N / d E [ T e V c m s s r ] injected signal 2extracted signal Energy [TeV] injected signal 3extracted signal Energy [TeV] E d N / d E [ T e V c m s s r ] injected signal 4extracted signal FB max , mismodeling of the GDE Figure B2.
Same as Fig. 8, except that here we assume the FB max model. MNRAS , 1–22 (2021)
TA sensitivity to the putative millisecond pulsar population responsible for the GeV excess Injected Luminosity [erg/s] ) T e s t s t a t i s t i c Detection threshold (5 ) Injected Luminosity [erg/s] R e c o v e r e d L u m i n o s i t y [ e r g / s ] injected signalrecovered MSP signalRecovered DM signal FB max , mismodeling of the GDE Figure B3.
Same as Fig. 6, except that here we assume mismodeling of the Galactic diffuse emission model and the FB max model.MNRAS000