Inverse Compton emission from millisecond pulsars in the Galactic bulge

TL;DR

This paper models inverse Compton gamma-ray emissions from millisecond pulsars in the Galactic bulge, comparing different spatial distributions to help identify their contribution to the observed GeV excess.

Contribution

It introduces a detailed GALPROP-based calculation of IC emissions from MSPs using a triaxial bulge model, highlighting spatial morphology differences for future detection.

Findings

Spectra are similar for different MSP spatial models

Spatial morphology features can distinguish MSP origins

High-energy IC emissions can inform MSP distribution models

Abstract

Analyses of Fermi Gamma-Ray Space Telescope data have revealed a source of excess diffuse gamma rays towards the Galactic center that extends up to roughly $\pm$20 degrees in latitude. The leading theory postulates that this GeV excess is the aggregate emission from a large number of faint millisecond pulsars (MSPs). The electrons and positrons ($e^\pm$) injected by this population could produce detectable inverse-Compton (IC) emissions by up-scattering ambient photons to gamma-ray energies. In this work, we calculate such IC emissions using GALPROP. A triaxial three-dimensional model of the bulge stars obtained from a fit to infrared data is used as a tracer of the putative MSP population. This model is compared against one in which the MSPs are spatially distributed as a Navarro-Frenk-White squared profile. We show that the resulting spectra for both models are indistinguishable, but that their spatial morphologies have salient recognizable features. The IC component above $\sim$TeV energies carries information on the spatial morphology of the injected $e^\pm$. Such differences could potentially be used by future high-energy gamma-ray detectors such as the Cherenkov Telescope Array to provide a viable multiwavelength handle for the MSP origin of the GeV excess.

Figures

Tables

Table 1. Table 1: The parameter values for the Galactic bar model S of Ref.  Freudenreich ( 1998 ) .

Parameters	Model S
Distance to the Galactic Plane $Z_0$ (pc)	16.46 $\pm$ 0.18
Bar Tilt Angle ${\theta }_0$ (deg)	13.79 $\pm$ 0.09
Bar $X$ Scale Length $a_x$ (kpc)	1.696 $\pm$ 0.007
Bar $Y$ Scale Length $a_y$ (kpc)	0.6426 $\pm$ 0.0020
Bar $Z$ Scale Length $a_z$ (kpc)	0.4425 $\pm$ 0.0008
Bar Cutoff Radius $R_{\text{end}}$ (kpc)	3.128 $\pm$ 0.014
Bar Cutoff Scale Length $h_{\text{end}}$ (kpc)	0.461 $\pm$ 0.005
Bar Face-On Shape $C_{⟂}$	1.574 $\pm$ 0.014
Bar Edge-On Shape $C_{\parallel }$	3.501 $\pm$ 0.016

Table 2. Table 2: Injection spectral parameters of MSP e ± superscript 𝑒 plus-or-minus e^{\pm} adopted in this work. See Eq. ( 13 ).

Model Name	$\Gamma$	$E_{\text{cut}}$
		(TeV)
Baseline	2.0	50
Inj1	1.5	50
Inj2	2.5	50
Inj3	2.0	10
Inj4	2.0	100

Table 3. Table 3: Propagation parameter setups considered in this study. Our baseline model corresponds to the best-fit propagation parameter in the scan for low-mass isotopes in Ref.  Jóhannesson  et al. ( 2016 ) , while the five additional models reflect the 95% confidence contours in their propagation parameter scan Jóhannesson  et al. ( 2016 ) (see text).

	$D_0$	$z_h$	$v_A$	$\delta$
	(1028 cm2 s-1)	(kpc)	(km s-1)
Baseline	6.330	9.507	8.922	0.466
Model 1	3.159	9.507	8.922	0.466
Model 2	7.006	9.507	8.922	0.573
Model 3	8.072	9.507	8.922	0.351
Model 4	2.748	3.000	8.922	0.466
Model 5	7.742	19.280	8.922	0.466

Equations

$$R_{\perp }^{C_{\perp }}$$

$$R_s^{C_{\parallel }}$$

$${\rho }_{\text{bar}}\propto \{\begin{array}{ll}{\text{sech}}^2(R_s),& R\le R_{\text{end}},\\ {\text{sech}}^2(R_s)e^{-\frac{(R-R_{\text{end}}{)}^2}{h_{\text{end}}^2}},& R>R_{\text{end}}.\end{array}$$

$${\rho }_{\text{NSC}}(R)=\frac{{\rho }_0}{1+{\left(\frac{R}{R_0}\right)}^n},$$

$${\rho }_{\text{NSD}}(r)=\{\begin{array}{ll}{\rho }_0{r}^{-0.1},& r<120\text{pc},\\ {\rho }_1{r}^{-3.5},& 120\text{pc}\le r<220\text{pc},\\ {\rho }_2{r}^{-10},& r\ge 220\text{pc}.\end{array}$$

$$\rho (R{)}_{\text{NFW}}=\frac{{\rho }_0}{{\left(\frac{R}{R_{\odot }}\right)}^{\gamma }{\left(\frac{1+R/R_s}{1+R_s/R_{\odot }}\right)}^{(3-\gamma )}},$$

$$\displaystyle\dfrac{d\psi}{dt}=\$$

$$\dot{p}=\int d(\log \nu )\frac{\nu U_{\nu }}{\hslash \nu }\frac{dp}{dt}(\nu ,\gamma ),$$

$$L_{e^{\pm }}=f_{e^{\pm }}\dot{E}=10f_{e^{\pm }}{L}_{\gamma }.$$

$$q(\vec{r},E)=\frac{c}{4\pi }{N}_0\frac{dN}{dEdt}\rho (\vec{r}).$$

$$N_0\int E\frac{dN}{dEdt}dE\int \rho (\vec{r})d{r}^3=L_{e^{\pm }}=10f_{e^{\pm }}{L}_{\gamma }.$$

$$I(l,b)\propto {\int }_{\text{l.o.s}}\rho (s,l,b)ds,$$

$$\frac{dN}{dEdt}\propto E^{-\Gamma }\exp (-E/E_{\text{cut}}),$$

$$B(r,z)=B_0\exp \left(-\frac{r-R_{\odot }}{R_0}\right)\exp \left(-\frac{z}{z_0}\right),$$

Full text

††preprint: APS/123-QED

I Introduction

In the past decade, the Fermi Large Area Telescope (Fermi-LAT) has provided accurate observations of the gamma-ray sky, of which the Galactic center (GC) remains one of the most intriguing and intricate regions. It is of paramount importance to map the gamma-ray emissions from the GC to better understand the properties of cosmic rays (CRs), the interstellar medium (ISM), and tests of dark matter (DM) in the inner regions of our Galaxy. Using template fitting techniques to regress out the Galactic and extragalactic diffuse emissions, multiple studies Goodenough and Hooper (2009); Vitale and Morselli (2009); Hooper and Goodenough (2011); Abazajian and Kaplinghat (2012); Gordon and Macias (2013); Macias and Gordon (2014); Hooper and Slatyer (2013); Abazajian et al. (2014); Daylan et al. (2016); Calore et al. (2015); Zhou et al. (2015); Ajello et al. (2016); Ackermann et al. (2017) have found an excess of gamma rays towards the GC extending up to $\sim\pm$20 degrees in latitude. Often referred to as the Galactic center excess (GCE), this excess emission has a centrally peaked spatial morphology that is roughly spherically symmetric with a radial power law of slope $\sim$2.4, and a curved energy spectrum peaking at $\sim$3 GeV. Some authors have argued that the GCE is consistent with a DM emission Goodenough and Hooper (2009); Abazajian and Kaplinghat (2012); Gordon and Macias (2013); Macias and Gordon (2014); Calore et al. (2015); Daylan et al. (2016) given its similarities to simplified predictions of weakly interacting massive particle DM models.

However, studies have also shown that the origin of the GCE can be explained by astrophysical sources, such as a population of unresolved gamma-ray-emitting millisecond pulsars (MPSs) Abazajian (2011); Abazajian and Kaplinghat (2012); Gordon and Macias (2013); Macias and Gordon (2014); Calore et al. (2015); Daylan et al. (2016). While there is ongoing debate regarding the consistency of the GCE with the luminosity function of MSPs measured elsewhere in the GalaxyCholis et al. (2015); Hooper and Mohlabeng (2016); Ploeg et al. (2017); Bartels et al. (2018a), evidence supporting the MSP hypothesis is mounting. For example, population synthesis simulations Gonthier et al. (2018) show that \sim$$10^{4} MSPs can inhabit the GC region and explain the GCE. Also, subthreshold photon-count statistics have shown detectable features that can be used to distinguish between the MSPs and DM interpretations of the GCE. In particular, Refs. Lee et al. (2016); Mishra-Sharma et al. (2017) introduced a new statistical technique called the non-Poissonian template fit which can be used for characterizing populations of unresolved point sources at fluxes below the detection threshold. They have shown that an unresolved population of point sources just below the sensitivity of the LAT is responsible for the GCE.

More recently, Refs. Macias et al. (2018); Bartels et al. (2018b); Macias et al. (2019) investigated whether the spatial morphology of the GCE is better described by the distribution of stars or by DM. It has been firmly established that the bulk of the stars in the GC region form a so-called box/peanut-bulge structure Dwek et al. (1995); Freudenreich (1998); Lopez-Corredoira et al. (2000); Nataf et al. (2010); Wegg and Gerhard (2013). The density distribution of these bulge stars can be reasonably described by a triaxial geometric function Freudenreich (1998); Nataf et al. (2010); Wegg and Gerhard (2013) extending in length to a few kpc from the GC. In addition to the box/peanut-bulge, there is a distinct stellar population in the innermost $\sim$200 pc of the Galaxy called the nuclear bulge (NB) Nishiyama et al. (2013). Reference Macias et al. (2018) used two different stellar maps for the box/peanut-bulge Freudenreich (1998); Ness and Lang (2016) and also considered the NB of Ref. Nishiyama et al. (2013), and demonstrated that such nonspherical bulge morphologies provide a significantly better fit to the GCE than a spherical Navarro-Frenk-White squared (NFW2) spatial map describing a DM annihilation signal. Importantly, that article showed that once a bulge model is included in the analysis, there is no longer statistically significant evidence for a NFW2 component. These results have been corroborated by Ref.Bartels et al. (2018b) using a different and more flexible analysis method Storm et al. (2017) and by Ref. Macias et al. (2019) using an improved Galactic diffuse emission model.

The population of MSPs in the Galactic bulge would not only produce prompt gamma-ray emission correlated with their spatial distribution, but also inject $e^{\pm}$ into the interstellar environment. These CR $e^{\pm}$ can produce secondary emissions by interacting with the ISM and magnetic field. It has been pointed out that while the prompt gamma rays from MSPs are expected to follow the morphology of the source distribution, secondary emissions—inverse Compton (IC), bremsstrahlung, and synchrotron radiation—are expected to have different morphologies, since they also depend on their relevant targets, i.e., the interstellar radiation field (ISRF), the gas distribution, and the magnetic field of the Galaxy, respectively. The IC component is a result of an energy-dependent convolution of the spatial morphology of the CR $e^{\pm}$ sources with the ambient photon fields from starlight, infrared (IR) light, and the cosmic microwave background (CMB).

Previous works Yuan and Ioka (2015); Petrović et al. (2015) have studied the secondary IC emission at $\sim$GeV-TeV energies from MSP $e^{\pm}$ interacting with the ISRF assuming a spherically symmetric distribution of MSPs. Under such an assumption, Ref. Lacroix et al. (2016) searched for secondary gamma-ray emission from MSPs by performing template fits to the GCE data, and found it to be difficult to detect or constrain the putative secondary IC component from an unresolved population of MSPs using Fermi-LAT data alone. However, as discussed above, there is now growing evidence for a significant departure from spherical symmetry.

Nonspherical source morphologies have been explored on small ($<$1 kpc) scales. For example, in a region overlapping with the NB called the Galactic ridge, Ref. Macias et al. (2015) showed that several different CR scenarios can explain the multiwavelength data taken from this patch of the sky. In addition, Ref. Carlson et al. (2016) evaluated the impact of star-forming activity in the Galactic ridge on the GCE properties. The study by Ref. Gaggero et al. (2017) solved the diffusion equation for CRs with a position-dependent diffusion coefficient and explained the recent H.E.S.S. Abramowski et al. (2016) measurements from this region by the interaction of the CRs with the gas in the central molecular zone (CMZ). Reference Guépin et al. (2018) explained the H.E.S.S. measurements by MSPs accelerating CR protons.

In this paper, we revisit the IC emission from MSP $e^{\pm}$ focusing on a potential nonspherical source morphology related to the Galactic bulge structure. Reference Abazajian et al. (2015) looked for and found a gamma-ray component following the IR distribution in the GC. The authors interpreted this as IC emission correlating with the distribution of optical light, but a propagation of the underlying $e^{\pm}$ was not performed. Here, we use the publicly available propagation code GALPROPStrong et al. (2000); Gal ; Strong et al. and make detailed predictions for the spectrum and morphology of the IC emission at 100 MeV-100 TeV energies. For the spatial morphology of MSP $e^{\pm}$, we use a three-dimensional (3D) stellar distribution model for the Galactic bulge obtained from a fit to IR data Freudenreich (1998) as well as one for the NB Launhardt et al. (2002) stars. We also show detailed comparisons against the spectrum and morphology obtained when the unresolved populations of MSPs are assumed to be spherically distributed. We find salient recognizable differences on large spatial scales and energies above $\sim$TeV, which can be used to distinguish source spatial morphologies.

This paper is structured as follows. In Sec. II we describe the 3D models used for the spatial distribution of the putative MSP population in the GC. In Sec. III we provide details about the propagation setup and configuration of our GALPROP runs. Details about the assumed MSP injection spectra are also given in this section. Our main results are shown in Sec. IV, where we show the predicted spectrum and morphology for the IC emission at $\sim$GeV-TeV energies. We illustrate how future measurements of diffuse gamma-ray emission with the Cherenkov Telescope Array (CTA) Wagner et al. (2009); Actis et al. (2011) have the potential to constrain the main properties of this purported MSP population at the GC. Finally, we conclude our study in Sec. V.

II Spatial distributions

To propagate the $e^{\pm}$ injected by MSPs, we need to model their spectrum and the spatial distribution of MSPs. We first discuss the spatial distribution. Two scenarios are considered: the stellar mass distribution in the Galactic bulge, and the spherically symmetric model for comparison.

II.1 Stellar models

We describe our first scenario, in which an unresolved population of MSPs is created in situ in the inner Galaxy and follows the distribution of stellar mass. Assuming that the same stellar populations responsible for the IR bulge trace the distribution of MSPs, a reasonable starting point for the MSP distribution is the bulge morphology itself. In principle, a proper morphological calculation would need to take into account the kick velocities of the MSP seeds at birth Eckner et al. (2018). However, the kicks experienced by MSPs should be lower than for isolated pulsars, which is also necessary for them to be confined to globular clusters. For example, while isolated pulsars are consistent with a Maxwellian velocity distribution with dispersion of 190 km/s Hansen and Phinney (1997), MSP estimates fall in the ranges $10-50$ Hooper et al. (2013); Cordes and Chernoff (1997), $85\pm 13$ Hobbs et al. (2004), or at the high end $130\pm 30$ km/s Lyne et al. (1998). Thus for our MSP estimates we do not consider the effects of initial kicks, and we consider the spatial distribution of stars using the 3D Galactic bulge model of Ref. Freudenreich (1998) and the NB model of Ref. Launhardt et al. (2002).

II.1.1 Galactic bulge

The mid- and near-IR signals of the inner Galaxy stars reveal a bar structure. The Galactic bar makes a tilt angle $\theta_{0}$ from the GC-Sun direction while the Sun is located at $\sim$10 pc above the Galactic plane. Here we adopt the bar model derived from data taken with the Diffuse Infrared Background Experiment instrument on board the Cosmic Background Explorer Freudenreich (1998). The shape of the bar is a generalized ellipsoid which can be parametrized as

[TABLE]

where $R_{s}$ is the effective radius; $a_{x}$, $a_{y}$, and $a_{z}$ are the scale lengths; $C_{\perp}$ and $C_{\parallel}$ are the face-on and edge-on shape parameters; and $X^{\prime}$, $Y^{\prime}$, and $Z^{\prime}$ are directions in the bar coordinate. Reference Freudenreich (1998) considered three different models for the radial dependence: model S, $\rho\propto\text{sech}^{2}(R_{s})$; model E, $\rho\propto\exp(R_{s}^{-n})$; model P, $\rho\propto[1+(R_{s}/R_{c})^{n}]$. We use model S here since it was the best-fit model found in Ref. Freudenreich (1998).

The models are truncated by a Gaussian function at the radius $R_{\text{end}}$ with scale length $h_{\text{end}}$. For model S, the density of the bar 111We notice that there was a typo in the argument of the $\exp$ function in Eq.(14) of Freudenreich (1998) which has been corrected in our Eq. 1. We have confirmed this in private communication with H. Freudenreich. is given by,

[TABLE]

The bar parameters used in our work are displayed in Table 1. These correspond to the best-fit values for model S using the so-called primary mask Freudenreich (1998). Recent studies of the bulge suggest larger tilt angles of $\sim 30^{\circ}$ Cao et al. (2013); Portail et al. (2017). However, we keep the best-fit angle from Ref. Freudenreich (1998) since it is consistent with the ISRF implemented in GALPROP v54. Also, as we discuss later, we do not expect the IC emissions to be very sensitive to this angle. Overall, the stellar mass of the Galactic bulge is $(1.4-1.7)\times 10^{10}$ solar masses $(M_{\odot})$ Bland-Hawthorn and Gerhard (2016).

II.1.2 Nuclear bulge

The NB refers to a dense stellar structure contained in the innermost region of the Galaxy. Associated with the CMZ, the NB has younger stars and undergoes active star formation, distinguishing it from the old and evolved stars of the Galactic bulge Launhardt et al. (2002). The NB makes up around 10% of the stellar mass in the bulge and its gamma-ray luminosity is comparable with that of the Galactic bulge Macias et al. (2018); Bartels et al. (2018b). The NB resides within the inner 230 pc of the GC and is made of two components:

Nuclear stellar cluster (NSC):

The NSC is a relatively small and very dense spherically symmetric structure in the innermost part of the NB. The stellar density in this region has been shown Launhardt et al. (2002) to be well described by a simple radial power-law function

[TABLE]

with best-fit power-law indices $n=2.0$ for $R\leq 6$ pc and $n=3.0$ for $R>6$ pc, with core radius fixed to $R_{0}=0.22$ pc. The stellar mass of the entire NSC is $(3\pm$ 1.5) $\times$ 107 $M_{\odot}$.

Nuclear stellar disk (NSD):

Surrounding the NSC is the NSD which makes up most of the stellar mass of the NB. The NSD is a cylindrical object with a radial dependence approximately described by a broken power-law function,

[TABLE]

The scale densities $\rho_{0}$, $\rho_{1}$, and $\rho_{2}$ ensure the continuity of the NSD density function. The density variation along the $z$ direction is given by an exponential cutoff with a scale height 45 $\pm$ 5 pc. The stellar mass of the entire NSD is (1.4 $\pm$ 0.6) $\times$ 109 $M_{\odot}$.

II.2 Spherically symmetric source

Although recent reanalyses of the GCE Macias et al. (2018); Bartels et al. (2018b) have shown that the Fermi-LAT data from the inner Galaxy prefer stellar maps to spherically symmetric ones, here for comparison purposes, we also model the putative MSP population at the GC with the square of an NFW density profile, of the form

[TABLE]

where we use a core radius $R_{s}$ = 23.1 kpc, Sun-GC distance $R_{\odot}$ = 8.25 kpc, and an inner slope of $\gamma$ = 1.20 Abazajian and Kaplinghat (2012); Macias and Gordon (2014). We label this as NFW2.

III Propagation

III.1 GALPROP code

We used the publicly available software package GALPROP v54 Gal ; Strong et al. in order to calculate the secondary gamma-ray emission from CR $e^{\pm}$ injected by MSPs. GALPROP is a numerical tool that solves the particle transport equations for a given source distribution and boundary conditions for all species of CRs. In particular, CRs can get accelerated by a multitude of different sources and then propagate long distances in the Galaxy. During propagation they produce secondary particles via interactions with the ISM and ISRF. Galactic diffuse gamma-ray emission is produced via $\pi_{0}$ decay, bremsstrahlung, and IC scattering, while lower-energy emission is produced via synchrotron radiation.

The CR $e^{\pm}$ transport equation is a partial differential equation of the form

[TABLE]

where $\psi=\psi(\vec{r},p,t)$ is the density per unit of the total particle momentum, $D_{xx}=\beta D_{0}(\frac{\rho}{\rho_{0}})^{\delta}$ is the spatial diffusion coefficient where $\rho$ is the rigidity of the $e^{\pm}$ and $\beta=v/c$, $\vec{V}$ is the convection velocity and $q(\vec{r},r)$ is the source term. Reacceleration is introduced as diffusion in momentum space with coefficient $D_{pp}$, which is related to $D_{xx}$ and the Alfvén speed $v_{A}$. We refer interested readers to the excellent review in Ref. Strong et al. (2007) for more details. We keep convection off, $\vec{V}=0$.

Given the source function of MSPs, the IC energy losses at each spatial bin are calculated by

[TABLE]

where $\nu$ and $U_{\nu}$ are the frequency and energy density of the ISRF, and $\gamma$ is the $e^{\pm}$ Lorentz factor. Very high-energy photons can interact with the photon fields and pair produce additional $e^{\pm}$ ($\gamma\gamma\to e^{+}e^{-}$). This is not included in GALPROP v54; however it only becomes important at $\sim$100 TeV, where the survival probability from the Galactic center to Earth reaches a minimum of 75-80% Moskalenko et al. (2006). At 10 TeV, the effect is reduced to the percent level.

GALPROP v54 contains dedicated routines to compute the propagation of DM annihilation/decay products and predict sky maps of secondary emissions. We modify the gen_DM_source.cc routine, which allows for user-defined source functions of the DM yields (DM profile and particle spectra), to model $e^{\pm}$ injected from MSPs. To compute the IC sky maps, we turn off the propagation of non-MSP CRs. As a first step, we made detailed comparisons of the results obtained with our modified GALPROP package against the literature Cirelli et al. (2014); Cirelli and Taoso (2016); Yuan and Ioka (2015); Petrović et al. (2015). We confirm that we are able to reproduce the gamma-ray spectrum and spatial profiles of Ref. Cirelli et al. (2014) as well as the synchrotron sky maps given in Ref. Cirelli and Taoso (2016) in the context of DM annihilations. Of greater relevance to this study are our detailed checks of the results in Refs. Yuan and Ioka (2015); Petrović et al. (2015). Although we were able to reproduce the gamma-ray spectral and spatial profiles in Ref. Petrović et al. (2015), we were only able to obtain the spectra given by Ref. Yuan and Ioka (2015). There are differences between our predicted spatial maps and those given in Fig. 3 of Ref. Yuan and Ioka (2015) for the same propagation setup. However, we believe these differences could be due to their using older two-dimensional (2D) ISRF maps in GALPROP.

III.2 Source function of MSP $e^{\pm}$

The spin-down energy $\dot{E}$ of MSPs is responsible for generating relativistic $e^{\pm}$ winds. The accelerations of these $e^{\pm}$ are limited when they lose energy via curvature radiation in the pulsar magnetosphere. Gamma rays are generated in this process. The gamma-ray efficiency $L_{\gamma}/\dot{E}$ is estimated to be about 10% on average Abdo et al. (2013). The $e^{\pm}$ that escape the magnetosphere via open field lines carry a fraction $f_{e^{\pm}}$ of the spin-down energy into the interstellar environment. The $e^{\pm}$ injection luminosity of MSPs is therefore related to the gamma-ray luminosity by

[TABLE]

The High-Altitude Water Cherenkov Experiment (HAWC) observations of Geminga and PSR B0656+14 in the TeV energy range suggest $f_{e^{\pm}}$ values of $\sim$7.2-29% Hooper et al. (2017). Constraints by the H.E.S.S. observations also indicate $f_{e^{\pm}}$ of the order 10% if about a thousand MSPs reside in the nuclear stellar cluster around Sgr A∗ Bednarek and Sobczak (2013).

We normalize the two spatial templates we implement in GALPROP (stellar and NFW2) by Eq. (9), using the best-fit gamma-ray luminosities obtained in Ref. Macias et al. (2018). The gamma-ray luminosities for the stellar template (which is a linear combination of the Galactic bar and NB models; see Sec. II.1) are $L^{\text{bar}}_{\gamma}$ = (1.4 $\pm$ 0.2) $\times$ 1037 and $L^{\text{NB}}_{\gamma}$ = (4.0 $\pm$ 1.0) $\times$ 1036 erg s*-1* for the bar and NB, respectively, while the NFW2 template (see II.2) has $L^{\text{NFW}^{2}}_{\gamma}$ = (1.7 $\pm$ 0.2) $\times$ 1037 erg s*-1*. These luminosities were taken from an analysis region of size 15∘ $\times$ 15∘ around the GC Macias et al. (2018).

The source term $q(\vec{r},E)$ (in units of MeV*-1* cm*-2* s*-2* sr*-1*) that is included in GALPROP can be written as the product of the injection spectrum $dN/dEdt$ and the source density distribution $\rho(\vec{r})$,

[TABLE]

The factor $c/4\pi$ is a convention in the GALPROP code. The source function is normalized by $N_{0}$, such that the integration over energy and volume matches Eq. (9),

[TABLE]

We will explore a range of $e^{\pm}$ spectra as detailed in Sec. III.3.

To show the source distribution and their relative injection intensities we produce injection intensity maps obtained by integrating along the line-of-sight direction the corresponding source density models,

[TABLE]

where $l$ and $b$ are the galactic longitude and latitude and $s$ is the distance from the Sun along the line-of-sight direction. Figure 1 (top panel) displays the injection intensities of the $\mbox{Galactic bar}+\mbox{NB}$ template in a 20∘ $\times$ 20∘ region around the GC. The resulting sky map is oblate and asymmetric. It is brighter for positive longitudes since the long axis of the bar makes a tilting angle with the GC-Sun direction. In contrast, Fig. 1 (bottom panel) shows that the NFW2 source is spherically symmetric and has a strong peak in the center of the Galaxy. We also present corresponding longitudinal and latitudinal profiles for the two models in Fig. 2. The aforementioned longitudinal asymmetry in the stellar distribution can clearly be seen in the top panel. Here it is also noticeable that the NFW2 intensity profile is more strongly peaked than the stellar template. The oblateness of the $\mbox{Galactic bar}+\mbox{NB}$ model is made manifest by comparing the tails of the profiles in both panels.

III.3 Injection spectrum

The maximum energy that MSPs can accelerate $e^{\pm}$ to is estimated to be $\sim 100$ TeV. In this work we adopt a similar $e^{\pm}$ spectrum to that used in Ref. Yuan and Ioka (2015) and explore a range of parameter values. Namely, we use a power law with an exponential cutoff,

[TABLE]

where $\Gamma$ is the spectral slope and $E_{\text{cut}}$ is the energy cutoff. We vary $\Gamma$ between 1.5 and 2.5, while $E_{\text{cut}}$ is varied between 10 and 100 TeV. The combinations of spectral parameters assumed in this work are listed in Table 2. Note that our choice of spectral parameters represents the entire MSP population in the Galactic bulge.

III.4 Configurations

In order to make simulations of CR propagation that are as realistic as possible, it is crucial to understand the diffusion parameters for all the CR species. Reference Jóhannesson et al. (2016) performed a scan of the parameter space of the CR injection and propagation. The scan was done separately for the low-mass isotopes ($p$, $\bar{p}$ and He) and the light elements (Be, B, C, N, O). Since each set of species has a different lifetime they probe different regions of the Galaxy. They found that the best-fit parameter setup for the low mass isotopes is different from that obtained for the heavier elements. Here, we adopt the propagation parameters for the low-mass isotopes of Ref. Jóhannesson et al. (2016). We account for uncertainties in the propagation parameters by also including the 95% credible contours provided in their 2D marginalized posterior distributions. The propagation setups used in this work are listed in Table 3.

We perform 3D GALPROP simulations using the standard ISRF data available with version 54 of the software package. The calculations are made for a Cartesian spatial grid with the Galactic plane placed in the $X$-$Y$ plane and the GC located at the origin of the coordinate system. The $X$ axis is defined by the GC-Sun direction.

In our simulations, the propagation volume extends to $\pm$20 kpc in both the $X$ and $Y$ direction. However, the CR halo height $z_{h}$ is set to different values depending on the model considered, as shown in Table 3. To trace the physics of the NB, we performed resolution tests for the spatial grid sizes. We fixed $\Delta Z$ = 0.125 kpc and found that the IC spectra and sky maps showed little difference when the resolution was $\Delta X$ = $\Delta Y$ = 0.125 or 0.25 kpc. We therefore adopted $\Delta X$ = $\Delta Y$ = $\Delta Z$ = 0.125 kpc for all our runs except for the $z_{h}=19.28$ kpc run, which was performed at a lower resolution of $\Delta X$ = $\Delta Y$ = 0.25 and $\Delta Z$ = 0.125 kpc due to computing memory demands. As a result, each run takes about 80 hours to finish using a computer cluster node with 500 GB memory and running at $\sim 6\times 10^{11}$ flop/s.

After choosing the resolution, Eq. (3) can be converted to the number of MSPs per spatial bin and normalized to reproduce the total number of MSPs that inhabit the Galactic bulge ($\sim 10^{4}$ Gonthier et al. (2018)). We find that every (0.125 kpc)3 bin corresponds to $\sim 40$ MSPs at the GC, and $\sim 4$ MSPs at 3 kpc along the long axis of the Galactic bar. We thus approximate the MSP distribution as a smooth function in our simulations.

We adopt the best-fit propagation of Ref. Jóhannesson et al. (2016) with the default $e^{\pm}$ spectrum of Ref. Yuan and Ioka (2015) as our baseline setup. In order to evaluate the impact of different propagation and spectral assumptions, we consider different propagation model setups in Table 3 and the different $e^{\pm}$ injection spectra listed in Table 2. For our baseline spectral setup we explore all variations in propagation setup. For our baseline propagation setup we explore all $e^{\pm}$ injection spectral combinations. In all our simulations the efficiency of the MSP $e^{\pm}$ is fixed to $f_{e}^{\pm}$ = 0.1.

III.5 Magnetic field

We adopt the default magnetic field from GALPROP, which is a double-exponential function,

[TABLE]

where $B_{0}=5\ \mu$G is the local magnetic field at the Solar System radius, and the scale parameters $R_{0}$ = 10 kpc, and $z_{0}$ = 2 kpc. The magnetic field strength of this model matches the 408 MHz synchrotron data Strong et al. (2000) and is in agreement with the total Galactic magnetic field estimates in the literature Heiles (1995); Beck (2001). However, the magnetic field at the center of the Galaxy remains uncertain. In particular, a multiband modeling on scales of 400 pc about the GC has produced a lower limit of 50 $\mu$G on the magnetic field strength Crocker et al. (2010), which the default GALPROP magnetic field does not obey (yielding, instead, a field strength of $\sim 10\ \mu$G). To this end, we test a modified magnetic field where we set $B=50\ \mu$G within a 400 pc region around the GC, but otherwise it matches the GALPROP default field everywhere else. The impacts of such a magnetic field will be discussed.

IV Results

Figure 3 shows the predicted IC emissions from the Galactic ridge region (defined as the $6^{\circ}\times 2^{\circ}$ area around the GC) for our stellar template (Galactic bar + NB). The solid black line is the expected flux from the baseline model. The shaded band represents the uncertainties resulting from changing the propagation setups as listed in Table 3. In the same plot we also show the expected flux from the NFW2 template (red dotted), which is indistinguishable from that from the stellar template given the uncertainties from the propagation parameters. This conclusion also holds for larger regions of the bulge of interest ($20^{\circ}\times 20^{\circ}$ around the GC). For the stellar template, we also show the IC emissions from the NB and the Galactic bar separately. For the Galactic ridge, the NB contributes $\sim 1/3$ of the total IC emission in the GeV energy range, and $\sim 1/2$ in the TeV range.

The GCE data Macias et al. (2015) at GeV energies, scaled to the Galactic ridge region, are shown by green triangles. The IC emission from MSPs is predicted to contribute less than $\sim 10$ % of the GCE. This is consistent with previous works that have not found evidence for secondary emission at the $\sim 1$ GeV energy range Lacroix et al. (2016). The IC fluxes extend to the TeV energy range before decreasing steeply at $\sim 10$ TeV. Meanwhile the predicted fluxes are below or within the range of H.E.S.S. observations (blue circles) Aharonian et al. (2006). Our resolution does not probe the smaller-scales of $0.2^{\circ}$–$0.5^{\circ}$ investigated in Ref. Hooper and Linden (2018). Also shown in the same figure are the differential sensitivities for 100 hours of CTA diffuse GC observations Silverwood et al. (2015) (dot-dashed black lines). We note that our predicted IC fluxes could be detected by forthcoming CTA observations. Reference Silverwood et al. (2015) computed the CTA sensitivities to a DM annihilation signal in the GC. Although the spatial morphology considered in their work is different than ours, we use it as a representative of the CTA sensitivities to a large diffuse signal in this sky region. That work also showed that by performing dedicated morphological analyses, the sensitivity of CTA can be improved by up to an order of magnitude for DM annihilation signals. We expect that similar improvements could be obtained by performing such morphological analyses with the models considered in our current study.

The IC spectra from the Galactic ridge region for different MSP injection spectra (Inj1, Inj2, …, Inj4) are shown in Fig. 4. We find that most of our predictions are at or below the corresponding H.E.S.S. data points (assuming $f_{e^{\pm}}=10\%$). However, for the hard injection spectrum model Inj1 ($\Gamma=1.5$ and $E_{\text{cut}}=50$ TeV), the IC spectra overshoots the H.E.S.S. observations. This means that either the Inj1 model is disfavored, or that $f_{e^{\pm}}$ is lower than $\sim 6\%$ in the Galactic ridge region. Contributing to this overshooting could be that the $e^{\pm}$ injection in the NB is overestimated. This is due to the fact that we normalize the $e^{\pm}$ luminosity by their gamma-ray luminosities. The best-fit NB gamma-ray luminosity obtained by Ref. Macias et al. (2018) may be somewhat overestimated, because the NB spatial morphology is similar to that of the CMZ structure which could cause spatial degeneracies in the fits. In this sense, a fraction of the gamma-ray photons that are of hadronic origin (emitted by the CMZ) could have been absorbed by the NB stellar template. Also, we note that at around 100 TeV, pair production will attenuate the predicted flux by $\sim 3/4$.

Figure 6 shows the morphologies of the IC emission for our baseline models in different energy windows: 3 GeV (left), 500 GeV (center), and 30 TeV (right). The top (bottom) row shows the sky maps for the stellar (NFW2) template. The sky maps are normalized by their fluxes at the GC. We find that there are energy-dependent morphological differences between the two IC predictions. These reflect the different $e^{\pm}$ source distribution models considered. In the GCE energy range ($\sim 3$ GeV, left panels), the IC sky maps are similarly elliptical for both the stellar and NFW2 templates. However, the sky maps become less elliptical at 500 GeV and above. At around $\sim 30$ TeV, the morphologies of the IC component start to show the source distributions displayed in Fig. 1. As it can be seen, in the highest energy window the sky maps for the stellar template (top-right panel) are boxy while that for the NFW2 template (bottom-right panel) is close to spherical. However, the left-right asymmetry due to the tilt of the bar is not seen in the IC emissions.

These features can be seen more clearly in the corresponding latitudinal and longitudinal profiles presented in Fig. 6. Here, we also show the variations in predicted IC morphologies when the propagation setups are varied, as in Table 3. We note that the morphological differences between the stellar and NFW2 templates are robust to changes in the propagation parameters. It is clear that at tens of TeV, the IC sky maps are sensitive to the source injection distributions.

At IC photon energies of $\sim 1$ GeV, the IC emission is in the nonrelativistic (Thomson) regime and the IC flux is dominated by up-scattered optical photons (Fig. 7). The optical photons are mainly emitted by stars whose density peaks along the Galactic plane. As a result, in this energy range the IC emissions are elongated along the Galactic plane and have marked elliptical appearances. This is almost invariant of the spatial morphology of the MSP distributions assumed, explaining the similarity between the left panels of Fig. 6. However, at higher energies, the IC emission starts to enter the relativistic scattering (Klein-Nishina) regime and becomes suppressed, causing the decline of the up-scattered optical photon signal from around 100 GeV (red dotted curve in Fig. 7). Since the Klein-Nishina regime is reached at higher energies for lower-energy target photons Ackermann et al. (2014), the IR and CMB photons continue, overtaking the optical photons. Furthermore, as the IR is less concentrated along the Galactic plane than the optical photons (and the CMB is isotropic), the IC emission retains more morphological information of the injected $e^{\pm}$.

A larger magnetic field means greater energy loss due to synchrotron radiation of the injected $e^{\pm}$ and a corresponding reduction in the IC emission. We tested a constant magnetic field $B=50\ \mu$G in the GC region on scales of 400 pc (see Sec. III.5) to explore its impact on our model predictions. Figure 8 shows the predicted IC spectrum for the baseline model assuming the default magnetic field (black solid curve with gray shaded band) and the corresponding IC spectrum for the same model but assuming the modified magnetic field (red solid curve). As can be seen, the IC spectrum normalization for the enhanced magnetic field setup is $\sim 1/2$ of the default one. This represents a change in normalization that is larger than our estimated modeling uncertainties from different propagation models (gray shaded). This reduction is mainly due to synchrotron energy loss of $e^{\pm}$ in the NB. Figure 8 displays the predicted IC spectrum for the NB component (red dashed curve) and Galactic bar (red dot-dot-dashed curve) separately in the enhanced magnetic field case. We observe that for this case the Galactic bar emission dominates. This is different to our predictions obtained using the default magnetic field setup, where the NB and bar contributions were comparable (see Fig. 3). This can be understood by noticing that the NB resides in the region where the modified magnetic field underwent the normalization increase. Consequently, much of the $e^{\pm}$ emitted from this region endure maximum energy loss via synchrotron radiation.

V Discussion and Conclusion

Recent analyses of the GCE Macias et al. (2018); Bartels et al. (2018b); Macias et al. (2019) have revealed its nonspherical nature and have provided further support for an MSP origin. We have revisited the computation of secondary IC emission from the $e^{\pm}$ injected by such MSPs. Compared to previous studies that assumed a spherically symmetric spatial distribution of MSPs, we adopted 3D models of the stellar distributions in the GC and numerically calculated the IC emissions using the GALPROP code. Furthermore, we systematically explored the impact of diffusion parameter uncertainties with additional GALPROP runs. We found that the predicted IC fluxes beyond 100 GeV are within the forecasted sensitivity limits of future gamma-ray telescopes for our baseline parameters (Fig. 3). The very high-energy IC emission from MSPs is nondegenerate with that caused by DM annihilation, from which the $e^{\pm}$ can only reach a few tens of GeV if DM is responsible for the GCE.

Although the IC spectra from the GC provide insufficient information for identifying the spatial model of the source, we found that the spatial morphology of the IC could serve as a discriminant between the spherically symmetric and 3D stellar distribution injection models (Fig. 6). In the GeV energy range, the IC morphologies are equally elliptical for both the stellar and NFW2 models. However, above $\sim$ TeV energies, they reveal morphological differences that trace the injection distributions. They can therefore be used to discriminate the spherically symmetric and 3D stellar injection models.

Our predicted IC fluxes contribute $\lesssim 10\%$ of the GCE emission and are at or below than the H.E.S.S. observations of the Galactic ridge at around a few TeV. These are consistent with the null detection of secondary emissions in the GeV range Lacroix et al. (2016) and the dominantly hadronic origins of the H.E.S.S. measurements  Gaggero et al. (2017); Abramowski et al. (2016). Thus they constitute important consistency checks of the MSP scenario for the GCE. We compared the IC fluxes with the CTA sensitivity from Ref. Silverwood et al. (2015) and found that the IC emission could be detected and potentially reveal a signature of GC MSPs with a specialized spectral and morphological search. The HAWC telescope DeYoung (2012); Abeysekara et al. (2013) operates in similar energy bands, and while not having a full view of the Galactic bulge region may have sensitivity to hard $e^{\pm}$ injection models. To this end, a wide field-of-view TeV gamma-ray observatory in the southern hemisphere is warranted Mostafa et al. (2018).

The detectability of the IC emission depends on the setup, including MSP $e^{\pm}$ spectrum and propagation parameters. We adopted a canonical $e^{\pm}$ power-law slope of 2.0, but softer spectra would make it a challenge to detect the IC component at TeV energies (Fig. 4). Furthermore, an increased magnetic field at the GC would reduce the IC emission via enhanced synchrotron energy losses. In particular, assuming a constant magnetic field of magnitude 50 $\mu$G in the inner 400 pc of the GC, we found a reduction of $\sim 1/2$ of the IC emission obtained with the default magnetic field setup. On the other hand, our model predictions were not very sensitive to changes in the propagation parameters within the 95% credible contours provided in the 2D marginalized posterior distributions of Ref. Jóhannesson et al. (2016).

There are various assumptions in our calculations that warrant future detailed studies. For example, we used the default ISRF from the GALPROP version 54, which is 2D after averaging the angular dependence. This has recently been updated in Refs. Porter et al. (2017); Jóhannesson et al. (2018), where a 3D ISRF was adopted in the context of CR propagation and high-energy gamma-ray emissions in the Galaxy. Their results show nontrivial impacts from employing the 3D ISRF on the propagation parameters of CRs and the gamma-ray intensity maps. Our results may be affected in many ways, including the fluxes and sky maps. Studies with the PICARD code show that new ISRFs increase gamma-ray intensities from the Galactic center, in particular at energies of $\sim 200$ GeV Niederwanger et al. (2019). Note however that our Galactic bar parametrization remains consistent with the ISRF of GALPROP. Even though more recent analyses suggest a larger tilt angle, we found that the left-right bulge asymmetry caused by the tilt is washed out in the IC sky maps, and is certainly smaller than the uncertainty caused by propagation parameters.

For the propagation part, we have adopted the results of a wide scan of propagation parameters Jóhannesson et al. (2016). However, caution must be exercised since the propagation properties around the Galactic center may be unique. We also have not covered all possibilities, e.g., we did not consider the possibility of cosmic rays advected out of the region by large-scale outflows. Such outflows may be related to the Fermi bubbles Crocker et al. (2015) and depending on the velocity would affect the secondary IC morphology. We have also neglected MSPs in the Galactic disk, which would provide additional $e^{\pm}$ injection and IC emission. However, population syntheses show that the MSP contribution to the Galactic diffuse gamma-ray emission is at the few-percent level or less Gonthier et al. (2018) and we do not expect this to substantially affect our results.

We have modeled the MSP population in the Galactic bulge. However, the presence of younger pulsars and the evolution of the pulsar population were not considered. It has been shown that TeV halos from younger pulsars can contribute to TeV emissions Hooper et al. (2018). This may potentially change the spectral property of IC emissions from the NB where active star formation is ongoing.

The magnetic field at the GC is a crucial parameter affecting the IC emission from a putative MSP population in the nuclear bulge. We have shown that an enhanced magnetic field in this region in turn augments the synchrotron energy losses of the MSP $e^{\pm}$, thus decreasing the IC yields. The estimated magnetic fields at $\sim$ 100 pc around the GC has large uncertainties and vary from 10 $\mu$G LaRosa et al. (2005) to 1000 $\mu$G Morris and Zadeh (1989). Here we only tested the original GALPROP model ($\sim 10\ \mu$G at the GC) and a 50 $\mu$G lower limit obtained in Ref. Crocker et al. (2010). An even larger magnetic field at the GC means that the synchrotron radiation would be dominant, especially for the NB component that resides within the 230 pc region around the GC. The spectrum and morphology of the IC emission from the Galactic ridge would potentially be changed by a strong magnetic field in this region. However, the effects on the larger-scale bar/bulge component are expected to be minor. On the other hand, we have only considered the 2D random magnetic field component. A recent study Orlando (2019) showed that the IC spatial maps can be significantly affected when more realistic 3D magnetic fields with both random and ordered components are included. This will apply also in the context of MSP secondary emission but its investigation is beyond the scope of the current study.

The Galactic center of the Milky Way offers a unique window to study novel astrophysical and dark matter signals. We have shown that the TeV energy range offers a new handle on the morphology of putative MSPs in the Galactic bulge responsible for the GeV excess. Telescopes such as CTA and HAWC South can be helpful for detecting these IC emissions and for constraining the origin of the GCE in the future.

Acknowledgements.

We thank Roland Crocker for careful reading of the manuscript and insightful suggestions. We thank Kev Abazajian and Manoj Kaplinghat for comments on the manuscript. D.S. and S.H. are supported by the U.S. Department of Energy under Awar No. de-sc0018327. This work was partially supported by the World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. O.M. acknowledges support by JSPS KAKENHI Grant Numbers JP17H04836, JP18H04340 and JP18H04578. The authors acknowledge Advanced Research Computing at Virginia Tech for providing computational resources and technical support that have contributed to the results reported within this paper. URL: http://www.arc.vt.edu