Strong Evidence that the Galactic Bulge is Shining in Gamma Rays

Strong Evidence that the Galactic Bulge is Shining in Gamma Rays

Oscar Macias Center for Neutrino Physics, Department of Physics, Virginia Tech, Blacksburg, VA 24061, USA Kavli Institute for the Physics and Mathematics of the Universe (WPI), University of Tokyo, Kashiwa, Chiba 277-8583, Japan GRAPPA Institute, University of Amsterdam, 1098 XH Amsterdam, The Netherlands    Shunsaku Horiuchi Center for Neutrino Physics, Department of Physics, Virginia Tech, Blacksburg, VA 24061, USA    Manoj Kaplinghat Center for Cosmology, Department of Physics and Astronomy, University of California, Irvine, Irvine, California 92697 USA    Chris Gordon School of Physical and Chemical Sciences, University of Canterbury, Christchurch, New Zealand    Roland M. Crocker Research School of Astronomy and Astrophysics, Australian National University, Canberra, Australia    David M. Nataf Center for Astrophysical Sciences and Department of Physics and Astronomy,
The Johns Hopkins University,
Baltimore, MD 21218
July 2, 2019

There is growing evidence that the Galactic Center Excess identified in the Fermi-LAT gamma-ray data arises from a population of faint astrophysical sources. We provide compelling supporting evidence by showing that the morphology of the excess traces the stellar over-density of the Galactic bulge. By adopting a template of the bulge stars obtained from a triaxial 3D fit to the diffuse near-infrared emission, we show that it is detected at high significance. The significance deteriorates when either the position or the orientation of the template is artificially shifted, supporting the correlation of the gamma-ray data with the Galactic bulge. In deriving these results, we have used more sophisticated templates at low-latitudes for the Fermi bubbles compared to previous work and the three-dimensional Inverse Compton (IC) maps recently released by the GALPROP team. Our results provide strong constraints on Millisecond Pulsar (MSP) formation scenarios proposed to explain the excess. We find that an admixture formation scenario, in which some of the relevant binaries are primordial and the rest are formed dynamically, is preferred over a primordial-only formation scenario at confidence level. Our detailed morphological analysis also disfavors models of the disrupted globular clusters scenario that predict a spherically symmetric distribution of MSPs in the Galactic bulge. For the first time, we report evidence of a high energy tail in the nuclear bulge spectrum that could be the result of IC emission from electrons and positrons injected by a population of MSPs and star formation activity from the same site.


I Introduction

The center of the Milky Way Galaxy is a complex environment displaying an extremely rich phenomenology of energetic cosmic-ray production and interactions. Many types of transient and persistent cosmic-ray source candidates populate the region, including a supermassive black hole (Genzel et al., 2010), active star formation (Krumholz and Kruijssen, 2015), supernova remnants (Ponti et al., 2015), pulsars and black holes (Muno et al., 2003), as well as new physics involving particle dark matter annihilations or decays (see, e.g., (Bertone et al., 2005) for a review). The bulk of the gamma-ray signal from the Galactic Center arises from the interactions of energetic cosmic rays with the interstellar gas and radiation fields. However, there is growing evidence that there exists an excess of extended gamma rays observed on top of models of cosmic-ray interactions in the Galactic Center. Dubbed the Galactic Center Excess (GCE), this signal has been detected by multiple analyses of the Fermi-LAT data (Goodenough and Hooper, 2009; Hooper and Goodenough, 2011; Boyarsky et al., 2011; Abazajian and Kaplinghat, 2012; Gordon and Macias, 2013; Macias and Gordon, 2014; Huang et al., 2013; Abazajian et al., 2014, 2015; Calore et al., 2015; Daylan et al., 2016; Ajello et al., 2016).

The origin(s) of the GCE remain elusive. The main properties of the GCE—being centrally peaked with a roughly spherically symmetric spatial morphology, having an energy spectrum peaking at a few GeV, and total luminosity of erg/s—can be accommodated within a scenario where weakly interacting massive particle (WIMP) dark matter of mass distributed in a slightly steep Navarro-Frenk-White (NFW)-like profile self-annihilate. However, the GCE might also arise from a population of gamma-ray emitting pulsars, in the form of either old “recycled” millisecond pulsars (Abazajian, 2011; Abazajian et al., 2014), or young pulsar remnants from Galactic Center supernovae (O’Leary et al., 2015, 2016). The gamma-ray spectra of known pulsars are similar to the GCE, and the capture of pulsars onto the Galactic Center region by Globular Cluster disruptions could explain the GCE’s quasi-spherical spatial distribution (Gnedin et al., 2014; Brandt and Kocsis, 2015). Other explanations for the GCE discussed in the literature include outbursts of cosmic-ray production by the central supermassive black hole (e.g., (Cholis et al., 2015a)).

While pulsars are a natural astrophysical candidate, there is ongoing debate regarding the consistency of the luminosity function required to explain the GCE with that measured for the pulsar population observed elsewhere as point sources (Cholis et al., 2015b; Hooper and Mohlabeng, 2016; Ploeg et al., 2017; Bartels et al., 2018a). Nevertheless, there are strengthening indications that the GCE may have an astrophysical origin. Cosmic-ray interactions in the Galactic Center region are notoriously complex to model. New analyses have cast doubt on some of the main properties claimed previously for the GCE. First, the energy spectrum of the GCE is subject to large systematic uncertainties arising from the incomplete understanding of cosmic-ray interactions in the Galactic Center region (e.g., (Abazajian et al., 2014, 2015; Calore et al., 2015; Ajello et al., 2016)). A variety of scenarios might, therefore, describe the GCE based on the energy spectrum alone. Second, it has been argued that the photon count distribution of the GCE is more consistent with arising from a population of faint point sources rather than being truly diffuse (Lee et al., 2016; Bartels et al., 2016). This observation supports an astrophysical origin in the form of faint point source over a microphysics origin in the form of dark matter particle annihilation. Third, there are growing reports of a spatial variation of the GCE energy spectrum with latitude (Horiuchi et al., 2016) showing a high-energy tail away from the plane, similar to that recently reported by Ref. Linden et al. (2016). Again, this disfavors a dark matter origin because particle properties would be expected to be spatially invariant (although secondary gamma-ray mechanisms might remain spatially-dependent).

Finally, it has recently been argued that the spatial morphology of the GCE is better matched to the asymmetric (stellar) bulge of the Milky Way than the spherically symmetric distribution expected for dark matter (Macias et al., 2018; Bartels et al., 2018b). The Galactic bulge is a triaxial bar-like structure in the central region of our Galaxy extending in size to a few kpc (Dwek et al., 1995; Freudenreich, 1998; Lopez-Corredoira et al., 2000). While the exact morphological properties of the bulge still remain under investigation, studies agree that it is asymmetric on our sky plane given that its major axis approaches us at positive Galactic longitudes (see, e.g., (Babusiaux and Gilmore, 2005)). Reference Macias et al. (2018) adopted two bulge maps to analyze the GCE: (i) The X-shaped Galactic bulge template of Ness & Lang (Ness and Lang, 2016), and (ii) the boxy shaped bulge of Freudenreich Freudenreich (1998). They found that the data overwhelmingly preferred either of the bulge maps over the spherically symmetric distribution expected by a dark matter origin. The bulge houses a broad mix of stellar populations from old to star forming (Garzon et al., 1997; Hammersley et al., 2000), and should contain ample candidates of astrophysical gamma-ray emitters.

In this paper, we put the reported bulge-correlation of the Galactic Center Excess to further tests. First, we include the recent GALPROP three-dimensional (3D) Inverse Compton (IC) maps constructed in Ref. Porter et al. (2017) in our pipeline. The new IC maps do not assume the Galactocentric cylindrical symmetry adopted in all two-dimensional (2D) IC maps constructed with previous GALPROP versions (v54 or older) and contain sophisticated model templates for the bulge/bar, spiral arms and stellar disk. Second, compared to previous work in Ref. Macias et al. (2018), our region of interest (RoI) is about a factor of larger. Third, we translate and rotate the bulge template to investigate how robustly they are detected in the Fermi-LAT data. We find that the Galactic bulge is detected strongly in gamma rays only when consistently positioned and aligned within one or two degrees of the known stellar bulge position Freudenreich (1998).

Fourth, we perform a morphological analysis of the stellar bulge distribution in order to test the different MSP formation scenarios that have been discussed in the literature. We find that an admixture formation scenario—in which some of relevant binaries are primordial and the remaining ones are the result of stellar interactions—is preferred to primordial-only formation with a confidence of . Finally, we present the gamma-ray spectra of the boxy bulge and the nuclear bulge. For the first time, we report evidence of a high-energy tail in the nuclear bulge spectrum that could be the result of IC emission from electrons and positrons injected by a population of MSPs and/or star formation activity from the same site. The dark matter implications of our new results will be presented in a separate article.

The structure of this article is as follows. In Sec. II we provide details about the data, model templates, and pipeline used in this study. In Sec. III we show our results for the empirical tests of rotation and translation of the stellar templates as well as our morphological analysis. In Sec. IV we present the gamma-ray properties of the Galactic bulge. Finally, a discussion and our conclusions are given in Sec. V.


width=1.0center Component Brief description Reference Interstellar gas correlated gamma-ray emission Hydrodynamical gas and dust maps divided in various rings: HI and H gas column density templates (4 rings each) plus two dust correction templates  Macias et al. (2018) Inverse Compton emission Seven IC maps are considered: a standard 2D IC map called “Std-SA0”, in addition to six different 3D IC models: F98-, SA0, F98-SA50, F98-SA100, R12-SA0, R12-SA50 and R12-SA100 Porter et al. (2017) Fermi bubbles Three templates are tested: Catenary, structured Fermi bubbles (SFB) and an inpainted version of the latter (SFB (Inp.))  Acero et al. (2016); Ackermann et al. (2017) Loop I Analytical model Wolleben (2007) Point sources Second Fermi Inner Galaxy Catalog (2FIG) Ajello et al. (2017) Sun and Moon templates Obtained with Fermi science tools Isotropic emission isoP8R2ULTRACLEANVETOV6v06.txt Nuclear bulge (NB) Template constructed from star count data  Nishiyama et al. (2013) Boxy bulge Model derived from a fit to near-infrared data  Freudenreich (1998) X-bulge Residual image from WISE data  Ness and Lang (2016) Spherically symmetric template Square of a Navarro-Frenk-White with a mild slope (NFW)

  • We adopt two different Galaxy-wide dust and stellar distribution models as considered in Ref. Porter et al. (2017): Freudenreich (1998) Freudenreich (1998) (F98) and Robitaille et al. (2012) Robitaille et al. (2012) (R12). For each, three different CR propagation setups are also considered: SA0, SA50 and SA100. See Table 3 of Ref. Porter et al. (2017).

Table 1: Summary of the model maps considered in this study.

Ii Data and Methods

We first describe our modeling procedure, starting with a description of data selection, the templates used in our analysis, and finally the analysis method.

ii.1 Data selection

We used years (August 4, 2008September 4, 2015) of Pass 8 ULTRACLEANVETO class photons with reconstructed energy in the 667 MeV158 GeV range. Photons detected at zenith angles larger than 90 were excised to limit the contamination from -rays generated by cosmic-ray interactions in the upper layers of the atmosphere. Moreover, we filtered the data using the recommended specifications (DATAQUAL0)&&(LATCONFIG==1). Fermi Science Tools v10r0p5 and instrument response functions (IRFs) P8R2ULTRACLEANVETOV6 were used for this analysis. In addition, the analysis was restricted to a square region of centered at Galactic coordinates with a spatial binning of .

ii.2 Baseline model templates

We started with a similar baseline model to that developed in Ref. Macias et al. (2018). In summary, this consists of the point-like sources listed in the 2FIG catalog Ajello et al. (2017), specialized templates for the Sun and the Moon, an isotropic component (isoP8R2ULTRACLEANVETOV6v06.txt), a standard 2D IC template generated with GALPROP v56 Porter et al. (2017), and an emission map for Loop I. In addition, the interstellar gas-correlated photons were modeled with a linear combination of atomic and molecular hydrogen gas templates divided into concentric rings (see Section II.7). Furthermore, these gas maps were obtained from a suite of hydrodynamical simulations of interstellar gas flow. This contrasts with gas maps in GALPROP and the official Fermi diffuse emission model that are constructed with an interpolation technique. Reference Macias et al. (2018) showed in detail that there are important morphological differences between the interpolated and hydrodynamic gas maps and that the latter provides a significantly better fit to the gamma-ray data in the Galactic Center region. See Table 1 for a summary of the model templates considered in this work.

ii.3 Bulge templates

Analyses Weiland et al. (1994); Binney et al. (1997); Freudenreich (1998) of the near-infrared emission measured by the DIRBE instrument on board of the COBE satellite were the first to uncover the non-spherical nature of the Galactic bulge. More recent studies of the Galactic bulge stars Cao et al. (2013); Wegg and Gerhard (2013) have firmly established that the bulk of bulge stellar mass forms a Boxy/Peanut (B/P) spatial morphology (this geometry is also sometimes called an “X-shape” in the literature). Using -body simulations of disk galaxy formation Bland-Hawthorn and Gerhard (2016), it has been found that the bulge hosts a rapid initial starburst component, as well as a second component which forms after disk build-up through dynamical instabilities Samland and Gerhard (2003). The former is thought to be associated with a classical bulge structure111Classical bulges are known to have roughly spherically symmetric morphologies. while the latter with the B/P structure. Studies Shen et al. (2010) of the internal kinematics of the stars of the Galactic bulge have determined that these are consistent with at least 90% of the bulge mass being contained in the B/P structure. Interestingly, the detection of a classical bulge component in the Milky Way has been claimed in Kunder et al. (2016). However, as noted above, this is expected to harbor a sub-dominant fraction of the Galactic bulge stellar mass.

Further evidence for an X-shaped or B/P bulge morphology was deduced from the bifurcation in red clump giant counts in the 2MASS and OGLE-III surveys Skrutskie et al. (2006); McWilliam and Zoccali (2010); Nataf et al. (2010). The analysis in Ref. Ness and Lang (2016) found, by constructing an image of the Milky Way bulge from an independent co-adding of publicly available WISE data, that the Milky Way bulge shows a distinct X-shaped structure in the projected stellar distribution that is different from previous works in that it has a more pronounced pinch-in effect at longitude ().

Assuming that the same stellar population responsible for the COBE/DIRBE infrared emission traces the spatial distribution of gamma-ray sources in the Galactic Center, then a reasonable guess for the gamma-ray morphology should be the boxy bulge structure. In particular, these Galactic bulge stellar structures are expected to be good tracers of MSPs and cosmic-ray (CR) acceleration sites in a more general sense.

ii.3.1 Boxy/Peanut bulge

Reference Freudenreich (1998) derived a parametric model for the spatial morphology of the Galactic bulge by fitting to COBE/DIRBE near-infrared ( m) data. The best-fitting model obtained in Ref. Freudenreich (1998) was Model S, which is a -squared function on the bar radial spatial profile. The density distribution of the boxy bulge for Model S is written as222Note that there was a typo in the argument of the function in Eq. (14) of Ref. Freudenreich (1998) which has been corrected in our Eq. (1). We have confirmed this in private communication with H. Freudenreich.


where is the effective radius:


and (, , ) is the Boxy bulge frame of reference. The bulge axis scale-lenghts are given by kpc, kpc and kpc, the bulge cut-off radius is kpc and the scale height kpc. In addition, the bulge face-on and edge-on shape parameters are and respectively.

In order to transform from the Galactocentric frame333The Galactocentric axis is assumed to be positive in the Sun-Galactic Center direction and the -plane is parallel to the plane of the Galaxy. , to the bulge frame (, , ), we need to perform two successive passive rotations of coordinates: the first is a clockwise rotation around the Galactocentric axis, and the second is a clockwise rotation of around the new axis . The Sun is assumed to be at pc and kpc in the Galactocentric frame. Lastly, we performed a line-of-sight integration of the density function in Eq. (1) and normalized the map appropriately for analyses within the Fermi Science Tools.

ii.3.2 X-shaped bulge

We reproduced the method described in Ref. Ness and Lang (2016) and applied it to the raw infrared WISE image of the Milky Way made available in Ref. Ness and Lang (2016). We employed the 3.4 (W1) and 4.6 (W2) micron WISE data and masked out pixels with negative flux, as well as the top and bottom 5% of pixels based on W1W2 color. We fitted an exponential disk model to the W1 and W2 bands and then subtracted the best fit models from the W1 and W2 data, respectively. Namely, we assumed


where   kpc, , , , and . We then applied a median filter of radius to the masked versions of both residual templates. In order to include this template in our gamma-ray pipeline, all pixels values with negative fluxes were set to zero in each median filtered exponential subtracted map. Finally, we computed the average of the two resultant maps and normalized it to flux unity for analyses with the Fermi Science Tools.

ii.3.3 Nuclear Bulge

Another stellar structure that is evident from infrared observations corresponds to the so called Nuclear Bulge (NB); this is thought to consist of a spherical Nuclear Stellar Cluster centered at the position of the supermassive black hole and a Nuclear Stellar Disk with radius pc and scale height pc. The gamma-ray emission from this stellar component was also studied in Ref. Macias et al. (2018). This template was constructed from a near-infrared stellar density measurement Nishiyama et al. (2013) of the central region of our Galaxy ( and ). More details about this template can be found in Ref. Macias et al. (2018).

Figure 1: Fermi bubbles templates considered in this analysis. Left panel: Flat Fermi bubbles map obtained in Ref. Acero et al. (2016). This map consists of two catenary curves describing the edges of the bubbles. This was the Fermi bubbles model considered in our previous GCE analysis Macias et al. (2018). Center panel: Structured Fermi bubbles template recovered with a spectral component analysis in Ref. Ackermann et al. (2017). The green contours represent the point source mask applied in that work. Right panel: Structured Fermi bubbles template after application of a Laplace inpainting method Lap () to correct for artifacts of the point source mask. The maps are normalized to unitary flux in the RoI and the colorbar shows the full range of variations in arbitrary units for the three panels.

width=1.0center Base Source Number of Reference source parameters for FB template baseline Catenary -2486188.1 -2486753.1 1130 15 Acero et al. (2016) baseline Structured FB -2486188.1 -2487322.3 2268 15 Ackermann et al. (2017) baseline Structured FB (Inpainted) -2486188.1 -2487802.1 3228 15 adapted from Ackermann et al. (2017)

  • The baseline model includes the 2FIG Ajello et al. (2017) point sources, Loop I, a standard 2D IC map predicted by GALPROP, hydrodynamic interstellar gas and dust maps divided in several rings, an isotropic component, and a template for the Sun and the Moon that matches the filters and cuts applied to the data (see also Table 1 for details). Three different models for the Fermi bubbles were considered. See Fig. 1 and Sec. II.3 for details. The statistical significance of a new source is given by TS, where maximum likelihood () values are computed independently for the Base and BaseSource models. Note that for both cases all parameters are maximized and so the will have additional parameters whose number is given in the second to last column of the table. The conversion between TS and is discussed in Ref. Macias et al. (2018).

Table 2: Summary of the likelihood analyses for Fermi bubbles maps.

ii.4 Spherically symmetric template

In addition to the stellar bulge templates, we also consider a map describing the potential gamma-ray emission from dark matter annihilation or a spherically symmetric distribution of MSPs. We model this by the square of a NFW density profile with an inner slope of 1.2, which has been shown to describe the GCE well in previous works (see Ref. Gordon and Macias (2013) for details on the precise parameters choices used in this work).

ii.5 Fermi Bubbles

Reference Acero et al. (2016) carried out a reanalysis of the low latitude counterpart of the Fermi bubbles (FB). The residual map found in that study has a spatial morphology whose boundaries are well described by two catenary curves of the form and for the Northern and the Southern bubbles, respectively. This flat FB template (hereafter simply called “Catenary”) is displayed in the left panel of Fig. 1.

Recently, Ref. Ackermann et al. (2017) derived an all-sky template for the Fermi bubbles using a spectral component analysis technique Malyshev (2012). This method assumes that the FB spectrum at high latitudes is the same as the spectrum at low latitudes in the energy range GeV. We show their “Structured FB” 444The Strucutured FB map used in this work corresponds to Fig.12 of Ref. Ackermann et al. (2017). The authors of that reference constructed the map with their 3 component spectral decomposition technique. template for our RoI in the central panel of Fig. 1. it is evident that this FB template contains image processing artifacts (see green contours in Fig. 1) due to point source mask template used in its derivation.

In order to ameliorate the artifacts present in the Structured FB map, we have applied an inpainting method to the FB image. The inpainting algorithm utilizes Laplace Interpolation Lap () which is a specialized interpolation technique for restoring missing data on a grid. The inpainted FB template is shown in the right panel of Fig. 1.

ii.6 Inverse Compton Models

Figure 2: Fractional residual difference map for the 3D IC Model “F98-SA0” with respect to the Standard 2D IC Model “Std-SA0” in the inner of the Galactic Center . The map is displayed at an energy of GeV with a pixel size of . The IC maps were computed with the most recent release of GALPROP (v56) Porter et al. (2017) and the new 3D ISRF data. No smoothing has been applied to the image. These maps reproduce the ones presented in Ref. Porter et al. (2017). See Sec. II.6 for a description of the 3D IC models.

Reference Macias et al. (2018) investigated the GCE systematics associated to the IC component with the use of the 2D Galactocentric cylindrically symmetric interstellar radiation field (ISRF) model attainable in GALPROP v54 Gal (). The authors of Ref. Ackermann et al. (2012) compared 128 different GALPROP models with all sky Fermi-LAT data finding a range of possible propagation parameter setups that are consistent with local CRs and gamma-ray measurements. Reference Macias et al. (2018) reconstructed a subset of the IC maps in Ref. Ackermann et al. (2012) that were found to have the largest morphological and spectral variations in order to estimate the impact of the IC component in the best-fit GCE emission.

In this article we improve upon the pipeline of Ref. Macias et al. (2018) by using the two recently developed 3D ISRF models made publicly available in the most recent release of GALPROP v56 Porter et al. (2017). These are based in two different Galaxy wide dust and stellar distribution models labeled F98 and R12: the F98 model is based on the work by Freudenreich (1998) Freudenreich (1998) and its bulge component is described in more detail in Sec. II.3, while the R12 is based on the study by Robitaille et al. (2012) Robitaille et al. (2012). The corresponding spectral intensities for the ISRF were calculated by Ref. Porter et al. (2017) with the Numerical Kalculation for Interstellar Emission (FRaNKIE) code. Although both models assume different stellar luminosities and dust densities, their predicted local intensities are consistent with near-infrared to far-infrared observations.

The IC maps investigated in Porter et al. (2017) utilized three different CR propagation setups (see Table 3 of that publication). These were labeled SA0, SA50 and SA100 according to the proportion of CR luminosity injected by the spiral arms. We have reproduced the main IC skymaps presented in that study and included them in our step-wise statistical procedure described in the following section. In particular, for each Galaxy model (F98 and R12) we have utilized the three proposed CRs propagation setups. These are named: F98-SA0, F98-SA50, F98-SA100, R12-SA0, R12-SA50 and R12-SA100. Figure 2 displays the fractional residual differences between a standard 2D IC map extracted from Ref. Ackermann et al. (2012) (Std-SA0) and the 3D model F98-SA0.

ii.7 Analysis pipeline

In Ref. Macias et al. (2018) some of us showed that the GCE spatial morphology was better explained by the stellar nuclear bulge and Boxy/Peanut (or X-shaped) bulge templates than by a spherically symmetric excess map given by an NFW profile. This was done by performing a set of maximum-likelihood fits summarized in Table I of that work Macias et al. (2018). In the present analysis we use a similar step-wise statistical procedure to establish whether a certain extended template for the GCE is statistically preferred over another.

The analysis pipeline used here to study the Fermi-LAT data utilizes the inner of the Galactic Center. Similar to Ref. Macias et al. (2018), we employed a bin-by-bin analysis method in which the full Fermi-LAT data were divided into 15 logarithmically spaced energy bins where, for each energy bin, we performed a separate maximum-likelihood fit with pyLikelihood within the Fermi Science Tools. In each of these maximum-likelihood runs we simultaneously fitted the normalization of all included diffuse emission components while keeping their spectral slope as a fixed parameter. Due to the small size of the bins, our results were not sensitive to the assumed spectral shapes of the extended and point-like sources. For simplicity we therefore used power-law spectra with fixed index of negative two.

The point source modeling was done with the 2FIG catalog Ajello et al. (2017). We considered only the point sources that were confidently found with the two different interstellar emission models considered in that work. There is a total of 288 such point sources inside our RoI. Due to limitations in the maximum number of parameters that can be reliably fitted in a given run with Fermi Science Tools, we have opted to follow a hybrid fitting approach analogous to that employed in Ref. Ackermann et al. (2017). Namely, we independently varied the normalizations of the 100 brightest 2FIG point sources in the RoI, while for the remaining point sources we constructed a single point source template in each energy bin. This is a reasonable approximation given that the amount of data utilized in the present study is very similar to the one used by the authors of the 2FIG catalog. To construct the combined point source template, we used the best-fit spectra in the 2FIG catalog and convolved it with the Fermi point spread function at each energy bin with the gtmodel tool. The resulting point source template was appropriately normalized and included in the fit along with the other extended and point-like sources. Other extended sources (W28, W30, RXJ1713, HESSJ1825-137, HESSJ1632-478 and HESSJ1837-069) that exist in our RoI were taken from the 3FGL catalog Acero et al. (2015) and varied independently in the fits.

The significance of each source was evaluated using the test statistic , where is the likelihoods of the background (null hypothesis) and is the likelihood of the source being tested plus background. The total TS of a source is given by where is the TS for bin and there are bins. Since a source has now 15 degrees of freedom (non-negative normalization in each energy bin) instead of two in a standard broad-band analysis (non-negative normalization and index), the detection threshold condition has to be adjusted accordingly. As discussed in Ref. Macias et al. (2018), this can be done by utilizing the formula given in case 9 of Self and Liang (1987):


where is the Dirac delta function, is a binomial coefficient with , and is a distribution with degrees of freedom. It follows (see Ref. Macias et al. (2018)) that a threshold of corresponds to a 4 detection of a new extended source.

Figure 3: Comparison of the log-likelihood values obtained for our six 3D-IC-Models vs a Standard 2D-IC Model. Blue (red) bars indicate which of the 3D-IC-Models considered improve (worsen) the fit to the inner region of the Galactic Center compared to a Standard 2D-IC Model (Std-SA0). The fits to the data include: the gas maps, 2FIG point source catalog, Isotropic, Loop I, Sun & Moon and inpainted Fermi bubbles template plus an alternative IC map at each different run. Notice that the different IC maps were not nested in the fits. See text for a description of the IC models evaluated in this work.
Figure 4: TS-values of the two main Galactic bulge models considered in this work: Navarro-Frenk-White squared map (NFW) (left) and nuclear bulge + boxy bulge (right). The reference model used in the computation of the TS-values is similar to the baseline described in Sec II.2, except that use of the inpainted structured map (far left black circle). Each bar shows the replacement of the Standard 2D map in the baseline model by a 3D IC map. The computation of the TS-values is done with a maximum-likelihood routine in which both the reference model (Base + FB(Inp)) and the bulge model (NFW or nuclear bulge + boxy bulge) are nested in the fit. More details can be found in Table 3.
Figure 5: Significance of the NFW template, in units, after the Galactic bulge template has been included to the baseline model. In other words, the fiducial model used in the computation of the NFW significances consists of the baseline described in Sec II.2, the inpainted structured FB map, and nuclear bulge + boxy bulge stellar map. Each bar displays a different 3D IC map as described in Sec. II.6. The computation of the TS-values is done with a maximum-likelihood routine in which both the fiducial model and the NFW are varied in the fit. We note that in none of the cases considered has the NFW template been significantly detected. More details can be seen in Table 3.

Iii Results


width=1.0center Base Source Number of source parameters baseline SFB (Inp.) -2486188.1 -2487802.1 3228 15 baselineSFB (Inp.) NFW -2487802.1 -2488619.8 1635 15 baselineSFB (Inp.) X-bulgeNB -2487802.1 -2488839.6 2075 baselineSFB (Inp.) Boxy bulgeNB -2487802.1 -2489230.1 2856 baselineSFB (Inp.)Boxy bulgeNB NFW -2489230.1 -2489233.9 8 0.9 15 baseline SFB (Inp.) -2487135.2 -2488729.2 3188 15 baselineSFB (Inp.) NFW -2488729.2 -2489119.4 780 15 baselineSFB (Inp.) X-bulgeNB -2488729.2 -2489242.2 1026 baselineSFB (Inp.) Boxy bulgeNB -2488729.2 -2489341.3 1224 baselineSFB (Inp.)Boxy bulgeNB NFW -2489341.3 -2489354.7 27 3.3 15 baseline SFB (Inp.) -2486191.5 -2488329.8 4276 15 baselineSFB (Inp.) NFW -2488329.8 -2489004.2 1349 15 baselineSFB (Inp.) X-bulgeNB -2488329.8 -2489120.2 1581 baselineSFB (Inp.) Boxy bulgeNB -2488329.8 -2489479.9 2300 baselineSFB (Inp.)Boxy bulgeNB NFW -2489479.9 -2489480.6 1 0.0 15 baseline SFB (Inp.) -2484807.7 -2487555.6 5496 15 baselineSFB (Inp.) NFW -2487555.6 -2488596.9 2083 15 baselineSFB (Inp.) X-bulgeNB -2487555.6 -2488700.0 2289 baselineSFB (Inp.) Boxy bulgeNB -2487555.6 -2489489.4 3868 baselineSFB (Inp.)Boxy bulgeNB NFW -2489489.4 -2489495.1 11 1.3 15 baseline SFB (Inp.) -2486304.7 -2488491.0 4373 15 baselineSFB (Inp.) NFW -2488491.0 -2489006.0 1030 15 baselineSFB (Inp.) X-bulgeNB -2488491.0 -2489093.6 1205 baselineSFB (Inp.) Boxy bulgeNB -2488491.0 -2489285.1 1588 baselineSFB (Inp.)Boxy bulgeNB NFW -2489285.1 -2489298.7 27 3.3 15 baseline SFB (Inp.) -2485197.9 -2488049.6 5703 15 baselineSFB (Inp.) NFW -2488049.6 -2488901.7 1704 15 baselineSFB (Inp.) X-bulgeNB -2488049.6 -2488927.8 1756 baselineSFB (Inp.) Boxy bulgeNB -2488049.6 -2489474.6 2850 baselineSFB (Inp.)Boxy bulgeNB NFW -2489474.6 -2489481.2 13 1.6 15 baseline SFB (Inp.) -2483983.0 -2487246.5 6527 15 baselineSFB (Inp.) NFW -2487246.5 -2488472.5 2452 15 baselineSFB (Inp.) X-bulgeNB -2487246.5 -2488479.3 2466 baselineSFB (Inp.) Boxy bulgeNB -2487246.5 -2489479.6 4466 baselineSFB (Inp.)Boxy bulgeNB NFW -2489479.6 -2489496.1 33 3.9 15

  • The baseline model in the first row set is described in Sec. II.2 and assumes the Standard 2D IC map. In the next row sets we change the 2D IC component by alternative 3D IC maps. This is denoted by baseline where with running through models {F98-SA0, F98-SA50, F98-SA100, R12-SA0, R12-SA50, R12-SA100} respectively. The Fermi bubbles model assumed corresponds to the Structured FB inpainted [SFB (Inp.)], see right panel of Fig. 1 and Sec. II.3 for details. Other bulge templates are: the Nuclear Bulge (NB) Nishiyama et al. (2013), boxy bulge by Freudenreich (1998) Freudenreich (1998), Ness&Lang X-bulge Ness and Lang (2016) and a Navarro-Frenk-White profile with (NFW). See Sec. II.3 for descriptions. The maximized likelihoods () are given for the Base and BaseSource models and the significance of the new source is given by TS. See also Table 2 for definitions.

Table 3: Summary of the main maximum-Likelihood analyses considered in this study.
Figure 6: Results of the rotation analysis. between the non-rotated boxy bulge and the boxy bulge template rotated with respect to the line in steps of () with significance ( for one extra parameter). This indicates that the unrotated infrared boxy bulge morphology is preferred by the gamma-ray data. The gamma-ray background and foreground model used to evaluate the maximized log-likelihood values of rotated boxy bulge templates consists of the baseline (for which the IC map corresponds to F98-SA0) described in Table 3, the inpainted structured FB and nuclear bulge map.

Using the initial baseline model described in Sec. II.2, the first step of our procedure is to find the best-fit normalization parameters for all free extended and point-like sources555The fitting of the 2FIG point sources in carried out with a hybrid approach described in Sec. II.7. in the RoI. The results are shown in Table 2. In this first iteration we seek to find out which of the three alternative Fermi bubbles maps (see Fig. 1) is preferred by the data. Three fits were done, each time adding a different FB template: Catenary Acero et al. (2016), Structured FB Ackermann et al. (2017) and Structured FB corrected with the inpainting method. As summarized in this table, the Structured FB (Inpainted) map was the one found to improve the fit the most. Therefore, for the next iterations we add this FB template and re-optimize the baseline model.

In the next stage of our procedure we perform six fits in which we replace the standard 2D IC map by one of the new 3D IC templates: F98-SA0, F98-SA50, F98-SA100, R12-SA0, R12-SA50 or R12-SA100 (which were recently proposed in Ref. Porter et al. (2017); see Sec. II.6). A comparison of the log-likelihood values obtained for each 3D IC template versus the standard 2D IC is displayed in Fig. 3 (larger is better). Interestingly, gamma-ray data from the inner prefer some of the 3D IC models to the Standard 2D one, with F98-SA0 being particularly favored. However, we note that although the new 3D IC maps are more physically motivated and much more realistic than the Standard 2D IC maps, the assumed CR propagation models and CR luminosities have not been tuned to match the gamma-ray data. Rather, they were obtained from a maximum-likelihood fit to local CR data Porter et al. (2017). So, in the remainder of this study we use the 3D IC maps in two different ways: first, we evaluate the impact of the new 3D IC maps on the results obtained in Ref. Macias et al. (2018). And second, we replace the Standard 2D IC model by F98-SA0 in our baseline when performing empirical tests of the bulge-correlated gamma rays (described following sections).

Table 3 shows the log-likelihood values obtained by adding to the fits the new extended templates (bulge models and preferred Fermi bubbles map). There are seven row sets in which for each different set an alternative 3D IC map is included in the baseline instead of the Standard 2D IC map: the first row set considers the Standard 2D IC model, the second F98-SA0, the third F98-SA50, and so forth, until the last row set that uses R12-SA100. Regardless of the IC map assumed, we find that a bulge model in the form of either the boxy bulge + NB, X-bulge + NB or an NFW is still needed by the data. As can be seen, this holds even after inclusion of the inpainted structured FB map. Figure 4 shows a comparison of the TS values obtained for the boxy bulge + NB and NFW models respectively.

In each row set shown in Table 3, the methodology is to keep the component with highest log-likelihood value and move to the next iteration. In all the cases considered in this study we find that the stellar models have larger TS values than the NFW map. Remarkably, when the boxy bulge + NB is included in the base model, there is no statistically significant need for an NFW template (the significance of the NFW map is always ). The stellar templates evidently absorb much of the residual photons and the data no longer need a map describing a dark matter distribution. Finally, when the NFW template is added to the base model, we find very strong evidence for the nuclear bulge + boxy bulge template. We therefore conclude that the data most prefer the nuclear bulge + boxy bulge. The low statistical significances obtained for the NFW model when different IC maps were assumed can be seen in Fig. 5. We note that this is consistent with our previous work Macias et al. (2018) where a smaller RoI and only the standard 2D IC maps were used.

iii.1 Rotation

Given the statistical preference for the stellar bulge templates under consideration, we continue to investigate whether the likelihoods can be improved by allowing more flexibility to the templates. In particular, we consider rotations and translations of the boxy bulge map. We even allow changes that are inconsistent with stellar data, which serve to test the interpretation of the gamma-ray correlations. Although the nuclear bulge is also robustly detected, we only perform these tests on the boxy bulge as this template is crucial to explain the GCE at higher latitudes. In order to test the rotational symmetry of the boxy bulge morphology we rotate the template about a line defined by (see inset image of Fig. 6) in increments of and compute the maximum log-likelihood at each turn. The results from this test are shown in Fig. 6 and indicate that the template in its original orientation is preferred. Note that the boxy bulge is East-West asymmetric (this is elegantly illustrated in Fig. 9 of Ref. Bland-Hawthorn and Gerhard (2016)), and thus the profile looks different for negative and positive rotations. This indicates that the data has constraining power to statistically distinguish any asymmetries present in the boxy bulge template. The maximum likelihood rotation was found to be with significance ( for one extra parameter). Therefore the data is consistent with the original orientation of the Boxy bulge image.

Figure 7: Results of the translation analysis along the longitudinal direction. values between the non-translated Boxy bulge template Freudenreich (1998) and the template artificially translated along the Galactic plane in steps of . There is evidence for a shifted center at (). The statistical significance for the shift is ( for one extra parameter), but systematic uncertainties (not included here) dwarf the statistical errors. See the caption of Fig. 6 for a description of the gamma-ray background model used.
Figure 8: Results of the translation analysis along the latitudinal direction. Same as Fig. 7, except that this time the translation analysis is applied in the latitudinal direction. There is evidence for a small offset in position of () with a significance ( for one extra parameter). See the caption of Fig. 6 for a description of the background model used.

It is instructive to compare the maximum likelihoods of the original Boxy bulge image for rotations. In the case of a clockwise rotation of , we find that the original orientation is a better model than the rotated one with a confidence of ( for one additional parameter). For the counterclockwise rotation of the same angle ( ), we find that the original orientation is preferred with a statistical significance of ( for one extra parameter). Additionally, the profile displayed in Fig. 6 shows that a rotation around is disfavored at ( for one extra parameter). This demonstrates that the data can also statistically differentiate the North-South bulge asymmetries.

iii.2 Translation

Next, we translate the boxy bulge template along the longitudinal and latitudinal directions in increments of and computed the log-likelihood value at each translated position respectively. The results from these tests are shown in Figs. 7 and 8.

In the longitudinal direction case (Fig. 7), we find evidence for a translation of , with a broad curve encompassing the range , while in the latitudinal case (Fig. 8) we find a translation at with a curve that is comparatively smaller (). The boxy bulge images translated by and are preferred with () and ( for one additional parameter) respectively. The preference for a shifted central position arises due to the presence of residuals on the west and north sides of the Galactic Center that are able to be absorbed by the boxy bulge template when it is translated along the latitude and longitude directions. Therefore, the translations appear real in the sense it is warranted by the gamma-ray data. However, in the construction of the 2FIG catalog Ajello et al. (2017), a group of several new point sources separated by less than were found around the region and . Due to those new point sources not being easily individually resolved, the procedure of the authors of Ref.  Ajello et al. (2017) excluded that group of new point source candidates from the 2FIG catalog. It is likely that this exclusion drives the bulge translation shift and that with their inclusion in the form of, for example, an empirical patch of positive residual emission at and , the boxy bulge latitudinal shift could be ameliorated. This is an interesting possibility that should be pursued in future work since it requires a dedicated focus along the Galactic plane, rather than the larger-scale Bulge emission of this study.

With regards to the longitudinal shift, recent gamma-ray analyses of the inner Galaxy Acero et al. (2016); Ackermann et al. (2017); Balaji et al. (2018) have also found excesses of gamma radiation off the Galactic Center at negative longitudes. Specifically, Refs. Acero et al. (2016); Ackermann et al. (2017) have both found that the base of the FB is a few degrees off the Galactic Center toward negative longitudes (see also central panel of Fig. 1). In addition, Ref. Balaji et al. (2018) has pointed out that the GCE is potentially offset by . However, when the authors run their pipeline with the Galactic plane masked, the offset becomes . The fact that this offset is sensitive to the applied masks indicates that the excess emission is mainly due to excess gamma-ray emission at very low latitudes () but that there is also some evidence for a shift of the gamma-ray bulge from high latitude data.

In this context, we note that the distribution of molecular clouds and young stars in the so-called “central molecular zone” (CMZ) is highly asymmetric Yusef-Zadeh et al. (2009); Immer et al. (2012); Koepferl et al. (2015); Longmore and Kruijssen (2018). Close to 2/3 of the molecular gas is at positive longitude—Galactic East near Sgr A—while about 2/3 of the young stars are at negative longitudes—Galactic West near Sgr A 666See for example the bright 24 m emission shown in Fig.1 of Ref. Longmore and Kruijssen (2018) Represented by blue color in that figure, this channel shows many bright point sources at which correspond to a population of massive young stellar objects and evolved high mass stars.. The asymmetric nature of the disposition of the molecular gas and young stars in the CMZ implies that the ISRF density distribution and the locations of CRs accelerators like supernovae and young massive stars, should be modeled with templates capable of accounting for these asymmetries. However, the triaxial function used in the derivation of the boxy bulge template adopted in this work Freudenreich (1998) is not able to account for deviations of the triaxial symmetry in the stellar bulge data and, by construction, traces mostly the oldest bulge stars Freudenreich (1998). Nonetheless, Ref. Cao et al. (2013) performed a translation analysis in their fits to red giant clump stars and found that the data always preferred an offset toward negative longitudes (see Table 1 in Ref. Cao et al. (2013)). While there are still quantitative differences, with the best-fit offsets of Ref. Cao et al. (2013) being smaller than those we find in the gamma rays, it is very likely that these are mostly due to gamma-ray correlations with star formation sites along the disk or in the far Galactic East region of the nuclear bulge. As a consequence, high level emission from these sites might be biasing our fits to larger shift values than those found by Ref. Cao et al. (2013). The asymmetry of the gamma-ray bulge is a very interesting topic that should be more carefully investigated in future work, using bulge templates that go beyond the functional forms adopted in this study 777 We observe that the 3D IC templates Porter et al. (2017) used in this work do not yet account for the stellar populations of the nuclear bulge and assume triaxial geometry (F98 or R12) for the bulge/bar structure. .

Iv Gamma-ray properties of the bulge

Figure 9: Results of the morphological analysis: Panels (a) and (b) show the latitudinal and longitudinal profiles of the boxy bulge + nuclear bulge best-fitting models summed over three energy bins between 1185.69 MeV and 2811.707 MeV. Note that in all cases the boxy bulge stellar density in Eq. (1) has been modified to include an additional slope parameter of the form . The red solid line () represents the MSPs distribution in the primordial formation scenario while the green dotted line () represents the dynamical formation scenario. Panel (c) displays the profile for different stellar density slopes () with respect to the primordial formation model (). The blue dashed line shows the best-fit slope which was found with a significance of ( for one extra parameter). The grey band corresponds to the statistical only errors.

iv.1 Gamma-ray bulge morphology

Globular clusters are old and extremely dense stellar agglomerations that might provide important clues about the formation mechanisms of MSPs in the Galactic bulge Abazajian and Kaplinghat (2012). There is strong observational evidence that some of the most massive globular clusters in the Milky Way contain large numbers of MSPs Abdo et al. (2013). However, given the relatively low escape velocities in globular clusters ( km/s) Gnedin et al. (2002), it is not straightforward to reconcile the presence of MSPs in these objects with empirical estimates of typical neutron star (NS) kick velocities ( km/s) Hobbs et al. (2005).

In order to solve this “NS retention problem” Pfahl et al. (2002) in globular clusters, two scenarios have been proposed: First is the primordial formation scenario Brandt and Podsiadlowski (1994), which posits that NSs form in massive binaries making them more easy to retain. This is because the kick-momentum transferred to the NS is shared with a massive companion, which leads to a reduction in the velocity of the post-supernova compact object. The number of MSPs is expected to scale linearly with the total stellar mass in this scenario. Second is the dynamical formation scenario, which rests on the observed positive correlations between the stellar encounter rate and the number of low-mass X-ray binaries (thought to be the progenitors of MSPs) in globular clusters Hui et al. (2010). Such observations indicate that after NSs are born, they can be captured by massive stars to form binaries which in turn can evolve to become MSPs. Since , where is the stellar density, it follows that in this picture the number of MSPs scales with the square of the stellar density.

In spite of the very different stellar densities and dynamical histories of globular clusters and the Galactic bulge, similar MSPs formation mechanisms might be operating in the Galactic Center Abazajian and Kaplinghat (2012); Abazajian et al. (2014); Gonthier et al. (2018). The spatial distribution of low-mass X-ray binaries in M31 provides evidence for the view that both dynamical and primordial components exist Voss and Gilfanov (2006), with the dynamical component dominating in the center. The close match between M31 X-ray binaries and the GCE morphology was previously used by Ref. Abazajian et al. (2014) to argue for unresolved MSPs being the dominant source of the GCE.

Here, we perform a morphological analysis of the GCE data to find out whether the spatial distribution of photons is better explained by the primordial or dynamical formation scenario. To this end, we modify the boxy bulge density function in Eq. (1) to include an additional slope parameter in the form . By varying in the range in steps of 0.1 and integrating the resulting density function along the line-of-sight, we create additional stellar maps which we individually pass through our pipeline routine. We note that in doing so the nuclear bulge map is left unmodified as we are focusing on the larger scale bulge component.

The results of the morphological analysis are presented in Fig. 9. Panels (a) and (b) show the longitudinal and latitudinal dependencies of the boxy bulge + NB models for various slope values: the primordial formation (), dynamical formation () and the best-fit slope (). There are appreciable differences between the modified boxy bulge + NB profiles and the NFW, which helps to explain why the maximum-likelihood analysis assigns largely different log-likelihood values to bulge models compared to spherically symmetric ones (e.g., Table 3). Interestingly, we find that the density slope that best fits the data is , which is preferred over the primordial scenario with a statistical significance of ( for one additional parameter), see panel (c) of the same figure. Therefore, we find significant evidence for an admixture formation scenario in which a fraction of the MSPs are formed through the primordial channel and the remainder are formed dynamically. Figure 9-(c) also shows that the fit worsens when stellar density models are assumed with slopes or . For instance, a slope is statistically disfavored with ( for one additional parameter) significance. Likewise, slopes of and are disfavored with () and () significance, respectively. This is an important consistency check that gives further support to the MSPs hypothesis for the GCE.

It has been argued that the putative population of Galactic bulge MSPs could have been the result of depositions from tidally disrupted globular clusters Gnedin et al. (2014); Brandt and Kocsis (2015); Fragione et al. (2018). However, a fundamental prediction of this model seems to be a spherically symmetric distribution of Galactic bulge MSPs Fragione et al. (2018), which we have demonstrated to be highly disfavored by the data. There is a growing consensus that the origin of the Boxy/Peanut bulge morphology is the result of the dynamical evolution of a disk stars via buckling instabilities Nataf (2017); Bland-Hawthorn and Gerhard (2016). In this picture, the Galactic bulge is predominantly composed of disk stars now on bar orbits. It is possible that by changing the initial conditions of their -body simulations, the disrupted globular cluster scenario could reproduce the Boxy/Peanut bulge morphology that the gamma-ray data requires. Nevertheless, such simulations would still need to explain the nearly isotropic distribution of the surviving population of Milky Way globular clusters. We note that accounting for these two seemingly irreconcilable observations in the disrupted globular cluster scenario Gnedin et al. (2014); Brandt and Kocsis (2015); Fragione et al. (2018) appears to be a difficult task, thus disfavoring current versions of this model for the GCE.

iv.2 Gamma-ray bulge spectrum

Figure 10: Boxy bulge spectrum: Shown are the boxy bulge fluxes (black points) and corresponding statistical (yellow) and systematic (green) errors. Systematic errors were estimated by calculating the dispersion of the best-fit bin fluxes obtained from fits that assumed the alternative 3D IC maps or the standard-2D IC model, respectively. Best-fit spectrum (black line) is given by with MeV cm s, and GeV ( for 15-3 dof, that is, ). This functional form was preferred over a simple power law (PL) with a statistical significance (=5.3 for one extra parameter. We compare the boxy bulge spectrum to the GCE obtained by assuming an NFW template in Ref. Calore et al. (2015) (red data points). Red (black) error bars correspond to their systematic (statistical) uncertainties. The data points were adapted by accounting for the solid angle of our RoI.
Figure 11: Nuclear bulge spectrum: Shown are the nuclear bulge fluxes (black points) and corresponding statistical (yellow) and systematic (green) errors bars. The blue data points correspond to H.E.S.S. Aharonian et al. (2006) data from the Galactic ridge region. Model A: A single power-law (PL) fit to the Fermi gamma-ray data from the nuclear bulge gives an acceptable fit ( for 15-2 dof, that is, ) with best-fit slope of and norm MeV cm s. However, we observe that when this PL fit is extrapolated to higher energies it overshoots the H.E.S.S. measurements from the same region. Model B: We have thus included the H.E.S.S. data points to the test, and obtained for 15+8-2 dof, that is, . A two-component model was preferred to a single power-law with a confidence of (=17.7 for four extra parameters). The best-fit spectrum is given by with MeV cm s, and GeV in addition to with MeV cm s, for GeV. Shaded regions represent the uncertainties on the best-fit model parameters.

The gamma-ray luminosity888In the computation of the luminosities we assumed 8.25 kpc for the distance to the Galactic Center. of the boxy bulge model () was found to be erg s for MeV, as estimated for the whole region of interest. A power-law with an exponential cut-off (PLE) model was preferred to a simple power law (PL) with a confidence of . As the amplitude is restricted to be non-negative, a distribution rather than the distribution is needed. We computed the best-fit spectral parameters with a fitting method to the inferred fluxes at each energy bin. We found GeV and , which as can be seen from Fig. 5 of Ref. Macias and Gordon (2014), are consistent with the resolved MSPs and also the GCE obtained using an NFW template. In Fig. 10, we show the spectra of the boxy bulge component in the inner around the Galactic Center. The yellow (green) error bars show the statistical (systematic) uncertainties. We have estimated the systematic uncertainties by computing the dispersion of the best-fit bin fluxes obtained when all the alternative 3D and standard-2D IC models are included in the maximum-likelihood analyses respectively. Our approach is consistent with earlier analyses by the Fermi-LAT collaboration Ajello et al. (2016) that argued for a large impact in the GCE spectrum due to IC modeling. After rescaling the GCE data from Ref. Calore et al. (2015) to our RoI and comparing to our boxy bulge spectrum (red data points), we find very good agreement with their results. We note that in that work the GCE was obtained by masking the inner region of the Galactic plane and assuming an NFW spatial model.

A fit to the nuclear bulge template preferred a two component spectral model over a single one. Initially, we attempted to fit a PL to only the Fermi-LAT data taken from the nuclear bulge. We found a best-fit slope of and normalization MeV cm s with ( for degrees of freedom). Although this p-value passes the usual acceptance condition , when the resulting PL formula is extrapolated to higher energies, it clearly overshoots the H.E.S.S. measurements from the same region (left-hand side panel of Fig. 11). We also tested a PLE spectrum and obtained ( for degrees of freedom ), which is less than the usual acceptance value.

Next, we included the H.E.S.S. data points Aharonian et al. (2006) to our test so as to make the model consistent with higher energy data from the same region. We therefore opted to add a PLE spectral model to account for data at around GeV and a PL model for GeV. This composite model consisting of a PLE+PL was preferred to a single PL with statistical significance (=17.7 for four extra parameters). Best fitting spectral parameters corresponding to the PLE component were found to be , GeV and erg s calculated in the energy range MeV. As for the PL component we found MeV cm s, for GeV 999This value is motivated by the known typical of the resolved MSPs. . While the H.E.S.S ROI ( and ) is indeed close to that of the nuclear bulge template, the H.E.S.S ROI could be slightly smaller. A dedicated observation of exactly the nuclear bulge region with H.E.S.S or the forthcoming Cherenkov Telescope Array (CTA) could make the consistency across the full gamma-ray energy range much better. Similarly to the boxy bulge component, the best-fit PLE parameters are consistent with spectrum of resolved MSPs. These results are summarized in the right-hand side panel of Fig. 11.

A simple astrophysical model that can match the characteristics of the preferred PLE+PL model is an unresolved population of MSPs in the nuclear bulge in addition to a population of energetic electrons from star formation activity and the same MSPs. Prompt gamma-rays from MSPs could account for the PLE spectrum seen at Fermi-LAT energies and a population of electrons injected from the same region could explain both the high-energy tail observed at energies GeV and the H.E.S.S. Galactic ridge spectrum (see for example Ref. Guépin et al. (2018) for a variation of this picture). Hadronic gamma-rays are also expected to contribute to the Galactic ridge at some level. A combination of hadronic and leptonic scenarios for this region have been explored in Ref. Macias et al. (2015). However, our new results for the nuclear bulge region motivates a much more detailed modeling that self-consistently accounts for hadronic gamma-rays from the CMZ region as well as a new population of electrons from MSPs and star formation activity. This particular topic will be addressed in future work.

It is useful to compare the mass-to-light ratios of our hypothesized boxy bulge and nuclear bulge MSP populations with others residing in different systems101010We caution that this comparison assumes that the ratio of MSP mass to stellar mass is constant across the different regions considered, which can only be approximately correct.. The stellar mass of the Galactic bulge is estimated to be around solar masses () and for the nuclear bulge . Bringing together the bulge mass estimates with our best-fit luminosities of the boxy bulge and nuclear bulge components, the luminosity per stellar mass for the boxy bulge is found to be erg s while for the nuclear bulge erg s . Thus, we find a luminosity per stellar mass for the nuclear bulge that is times higher than for the boxy bulge. This is qualitatively consistent with a picture in which the MSPs of the boxy bulge are systematically older than those in the nuclear bulge, and hence less luminous per stellar mass. The higher luminosity per stellar mass of the nuclear bulge might also be due to additional energy being injected from star formation activity on top of an old MSP population.

V Discussion and conclusions

In this work we have performed a detailed reevaluation of the systematic uncertainties affecting the properties of the Galactic Center Excess (GCE). We have utilized hydrodynamical interstellar gas and dust maps Macias et al. (2018) in order to account for uncertainties associated with gas-correlated gamma-ray emission. In addition, we have employed, for the first time, the new 3D inverse-Compton (IC) maps111111Some of the new 3D IC maps dramatically improve the total log-likelihood of baseline model. This demonstrates the great constraining potential of the Fermi-LAT data from the inner Galaxy for constructing improved CRs and ISRF models. that abandon the assumption of Galactocentric cylindrical symmetry adopted in most previous analyses of the same region. Improved low-latitude Fermi bubbles maps have also been included. By considering this wide range of maps and uncertainties, we have confirmed recent studies Macias et al. (2018); Bartels et al. (2018b) that reported a positive correlation between the stellar distribution of the Galactic bulge and the GCE.

A recent study by the Fermi collaboration Ajello et al. (2016) has raised some concerns with regard to whether the GCE is the result of incomplete modeling of the IC component in the inner Galaxy. However, our results reveal that such concerns are unwarranted. Even though the new 3D IC maps Porter et al. (2017) include the bulge/bar for the CR source and ISRF densities, an additional stellar bulge component in the form of the nuclear bulge + boxy bulge is still required by the Fermi-LAT data. This is consistent with the hypothesis of an unresolved population of MSPs in the Galactic Center region. Since CRs interacting with interstellar gas and the ISRF produce IC morphologies at GeV that are different to both those of the CR source distribution and the target material, statistical analyses of the GCE would still be required to add template model maps to account for the prompt gamma-ray emission from MSPs that closely match their source distribution.

Although the X-shape Ness and Lang (2016) structure of the Galactic bulge is estimated to contain some % of the mass of the bulge stars, these stars spend only a small fraction of their orbital time within the X arms Bland-Hawthorn and Gerhard (2016); Nataf (2017). Therefore, at any given time, about 5–10% of bulge stars are observed to be within the X arms Cao et al. (2013); Nataf (2017); Simion et al. (2017). Under the assumption that gamma-ray sources are kinematically similar to stars, we would therefore expect that the X arms would contribute of the Galactic bulge gamma-ray luminosity. In the maximum-likelihood analyses presented in Table 3, we found that there is statistical evidence for the X-shape bulge in gamma-rays, but that this was weaker than that of the whole boxy bulge structure. In our step-wise statistical procedure, only new sources that increased the log-likelihood the most were allowed to be included in the subsequent baseline model.

In Ref. Macias et al. (2018), using a RoI of size and only the standard 2D IC maps, some of us showed that the GCE was better fit by either the nuclear bulge + X-bulge or nuclear bulge + boxy bulge. In that work both stellar templates were found to have similar statistical significances. Now with our larger RoI and improved IC maps used in this analysis, we are finding that the significance of the nuclear bulge + X-bulge is smaller than that for the nuclear bulge + boxy bulge template. However, we note that it is not at present possible to nest the boxy bulge template with the X-bulge template in the same maximum-likelihood fit, which is necessary for a robust comparison. This is because the boxy bulge model already contains a substantial fraction—if not all—of the X-bulge stars. Our new results strongly encourage a careful revaluation of the luminosities associated to the X-bulge and boxy bulge stellar components using stellar maps that are physically self-consistent.

Figure 12: Fermi Bubbles spectrum: Inferred bin-by-bin fluxes (black points) and corresponding statistical (yellow) and systematic (green) error bars. The systematic errors show the flux dispersion when the maximum-likelihood analyses are repeated using the alternative 3D and 2D IC maps presented in Ref. Porter et al. (2017). For comparison purposes, we show the butterfly plot of the boxy bulge spectrum of our Fig. 10. The red points are the overall high-latitude () flux extracted from Fig. 18 of Ref. Ackermann et al. (2014).

In order to test whether the stellar models can be improved by allowing more flexibility to the templates, we rotated the boxy bulge template around the central position in increments of and computed the log-likelihood at each orientation. The results from this test indicate that the boxy bulge in its original orientation is preferred and that the data is sufficiently constraining to differentiate the East-West and North-South asymmetries present in the stellar bulge model (see Fig. 6).

Figure 13: Best-fit emission generated by MSPs formed through the dynamical and primordial formation channels: See also Fig. 9 for descriptions. Panels (a) and (b) show the relative contribution of the primordial and dynamical MSPs formation channel. This was obtained by simultaneously including in the fit a and map and then estimating their corresponding best-fit gamma-ray luminosities. When the flux fractions calculated this way are superimposed, the total flux profile is approximately described by , in agreement with results shown in Fig. 9. The red solid, gren dotted and black dashed lines show the best-fit nuclear bulge, primordial formation and dynamical formation fractions, respectively. We find that primordially formed MSPs contribute of the Galactic bulge luminosity.

Similarly, with the aim of testing the GCE for any possible translational shifts with respect to the location of the stellar bulge, we have moved the central position of the boxy bulge along latitude and longitude in steps of . Statistically significant shifts were found toward positive latitudes at and negative longitudes at (see Figs. 7 and 8). The former is relatively small and is likely due to a positive excess resulting from the suppression of a cluster of point sources in the 2FIG Ajello et al. (2017) catalog around the region and . However, the later prefers a shifted center at larger angular distance and increased statistical confidence. This result is in line with a similar such analysis performed with the red clump giant stellar data as seen in Table 1 of Ref. Cao et al. (2013), where a westward translation of around 0.3 was found. Although the preferred shifted positions in the two studies differ by , it is worth noting that in the region around and there is a population of young massive stellar objects that might be emitting significant gamma-ray flux which, by construction, the analysis of Ref. Cao et al. (2013) does account for. This is because the red clump giants are systematically older stellar objects. In our view the results of these empirical tests emphasize the validity of the proposal for a stellar origin for the GCE, within the systematic limitations of the current maps. In summary, we find that the gamma-ray data highly prefers a Galactic bulge model that is consistently positioned and aligned within one or two degrees of the known Galactic bulge Freudenreich (1998). All other translational and rotational variations result in significantly poorer fits to the data.

Perhaps the most striking result to emerge from the GCE data is that any model predicting a spherically symmetric morphology is robustly ruled out (see Fig. 5). This, of course, includes the annihilation of weakly-interacting dark matter models but also current versions of the disrupted globular cluster scenario Gnedin et al. (2014); Brandt and Kocsis (2015); Fragione et al. (2018). Any successful disrupted globular cluster model must explain both the spherical distribution of surviving globular clusters in the Galaxy as well as the boxy/peanut distribution of MSPs in the bulge. Taken together, these two observations seem to disfavor this scenario (at least in its current form Gnedin et al. (2014); Brandt and Kocsis (2015); Fragione et al. (2018)) as the main formation mechanism for MSPs in the Galactic bulge.

We have also tested the primordial and dynamical formation scenarios for the Galactic bulge MSPs. In our detailed morphological analysis of the GCE, we have found a significant statistical preference () for an admixture formation scenario in which a fraction of the MSPs are formed through the primordial channel () and the remaining one through the dynamical channel. Due to systematic uncertainties in the Galactic diffuse emission model, our maximum-likelihood method assigns almost identical log-likelihood values to a pure primordial formation as well as to a pure dynamical formation scenario (see Fig. 9). However, models whose spatial morphology start to depart from those predicted by these two scenarios are statistically disfavored by the data. This is an important check that lends further support to the MSP theory for the GCE because we expect both components to exist in the MSP population Voss and Gilfanov (2006).

Reference Eckner et al. (2018) used the results of a correlation between globular cluster gamma-ray luminosity and the stellar encounter rate Hui et al. (2010) to infer the energetics and morphology of the dynamical formation scenario for MSPs in the Galactic bulge. That study predicts a contribution to the bulge luminosity relative to the primordial gamma-ray luminosity at the percent level (). In contrast, using our maximum-likelihood procedure we find that the best-fit contribution for this channel is (see Fig. 13). This was obtained by simultaneously including in the fit a and map and then estimating their corresponding best-fit gamma-ray luminosities121212The fractional luminosity error was obtained by adding in quadrature the inferred energy fluxes at each energy bin. We also note that in this case there are 15 new degrees of freedom (one for each band of the template) while in estimating there was only one new degree of freedom.. Figure 13 shows how the flux fractions computed this way adds up to approximately our best-fit stellar density slope obtained with the slope scan of Fig. 9. Although our detailed morphological analysis does not seem to corroborate the luminosity fraction predicted by Ref. Eckner et al. (2018), we caution that a more direct comparison with that study is not at present possible given that the authors have assumed a spherically symmetric model (see Sec. 3 in Eckner et al. (2018)) for the distribution of the bulge stars. Such stellar distributions are thought to be not realistic Bland-Hawthorn and Gerhard (2016); Nataf (2017) and as showed are robustly ruled out by the gamma-ray data Macias et al. (2018); Bartels et al. (2018b).

In agreement with most previous analyses of the GCE Abazajian and Kaplinghat (2012); Gordon and Macias (2013); Macias and Gordon (2014); Huang et al. (2013); Abazajian et al. (2014, 2015); Calore et al. (2015); Daylan et al. (2016); Ajello et al. (2016), we have found that the spectrum of the boxy bulge is consistent with that expected from an unresolved population of MSPs in the Galactic bulge (see Fig. 10). In particular, we confirm earlier results Linden et al. (2016) that claimed the detection for a high energy component above GeV associated to the GCE. As can be seen in Fig. 12, this high-energy tail has a different spectrum from that of the Fermi Bubbles spectrum. Such a component could be the result of IC emission of electrons and positrons injected by the same Galactic bulge MSP population Linden et al. (2016).

Furthermore, we identified a corresponding high-energy tail also in the nuclear bulge template (Fig. 11). As far as we are aware, this is the first time that a high-energy tail has been detected in this smaller angular region. A fit to the Fermi and H.E.S.S. observations from this site preferred a two component spectral model (with a confidence of ) comprising a power-law with an exponential cutoff model plus a simple power-law for energies GeV. An important implication of this result is that energetic electrons emitted by MSPs, and possibly star forming activity, could potentially contribute a significant fraction of the TeV gamma-rays detected by H.E.S.S.

Our results show that the Galactic Center Excess in gamma rays is correlated with the bulge of the Milky Way. The most likely explanation for this correlation is that the unresolved bulge millisecond pulsars are responsible for the Galactic Center Excess.

We are grateful to Kevork Abazajian, Shin’ichiro Ando, Anthony M. Brown, Henry Freudenreich, Ryan Keeley, Dmitry Malyshev, Harrison Ploeg, Troy Porter, Nicholas Rodd, Tracy Slatyer, and Deheng Song for useful comments and discussions during the project. This work was supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. OM was partially supported by JSPS KAKENHI Grant Numbers JP17H04836, JP18H04340 and JP18H04578. SH thanks support by the U.S. Department of Energy under award number de-sc0018327. MK was supported by NSF Grant No. PHY1620638. DMN was supported by the Allan C. and Dorothy H. Davis Fellowship. OM acknowledges the high performance computer center at Kavli IPMU for providing computational resources and specially to Tatsuhiro Yano for technical support that have contributed to the results reported in this paper.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description