FERMI-LAT OBSERVATIONS OF HIGH-ENERGY γ-RAY EMISSION TOWARD THE GALACTIC CENTER

The diffuse γ-ray emission is produced by the interaction of high-energy CRs with the interstellar gas and radiation fields. The limited angular resolution and statistics that are a characteristic of high-energy γ-ray data, coupled with the relatively intense interstellar emission at low latitudes, make accurate modeling of the latter important for characterizing all but the brightest point sources there (Abdo et al. 2010a; Nolan et al. 2012).

Two analysis approaches for studying the interstellar emission have been used by the Fermi-LAT Collaboration in previous works. Templates tracing γ-ray emission processes were used to determine the γ-ray emissivity of the interstellar gas within several kpc of the Sun (e.g., Abdo et al. 2009a, 2010c; Ackermann et al. 2011). The GALPROP code was used in an extensive study of IEMs constrained by local CR data and their correspondence with the Fermi-LAT data (Ackermann et al. 2012a). There are merits to both approaches. Fitting templates allows for fairly robust extraction of physical quantities, but is a method that is constrained by the assumption that interstellar medium (ISM) densities and other properties (gas-to-dust ratio, X_CO-factor, etc.) and CR spectra remain constant throughout the template, and that a suitable template is available (e.g., the inverse Compton [IC] component of the interstellar emission must be obtained using modeling codes)—see Strong et al. (1988) and Strong & Mattox (1996), and references therein. The GALPROP code can be used to predict the diffuse γ-ray emission throughout the Galaxy, and is capable of reproducing the observations at the ∼20% level. But the predictions of the propagation model based approach are limited by the quality of the inputs to the model calculations, which include the spatial distribution of CR sources and their injection spectra, and the spatial distribution of the interstellar gas density and the interstellar radiation field (ISRF) energy density. In this paper a combination of these methods is used where the GALPROP code is employed to predict templates for the interstellar emission that are fit to the γ-ray data to estimate the foreground and background emission toward the inner Galaxy.

The results of the study by Ackermann et al. (2012a) are used for the baseline IEMs, which are further fit to the Fermi-LAT data. As a reminder, the Ackermann et al. (2012a) study used a grid of IEMs based on diffusion-reacceleration CR propagation models. The spatial distribution of CR sources, the H i spin temperature, H i column density corrections from dust emission, and the size of the CR confinement volume were the fixed parameters in the grid. For each grid point the diffusion coefficient was obtained by adjusting it to reproduce the observed CR secondary/primary ratios iterating with a fit to the γ-ray data for the X_CO distribution for each CR source model. The γ-ray emission for each of the IEMs in the grid was then compared with the Fermi-LAT data in the 200 MeV to 100 GeV energy range. The models in the Ackermann et al. (2012a) study agree at the ∼10%–20% level with the LAT observations over the sky.

A major uncertainty affecting predictions of the interstellar emission toward the inner Galaxy is the spatial distribution of CR sources. The Yusifov & Küçük (2004) pulsar distribution ("Pulsars") and the distribution of OB-stars ("OBstars"; Bronfman et al. 2000) encapsulate this because they represent reasonable extremes for the Galactocentric radial dependence. Figure 1 shows the Galactocentric radial distributions of these CR source models. The Pulsars distribution is non-zero at the GC while the OBstars distribution goes to zero near ∼2 kpc. The models⁷⁰ assume an axisymmetric cylindrical geometry for the CR confinement volume with a halo height z_h = 6 kpc and maximum radial boundary R_h = 30 kpc. This halo height is the closest in the IEM grid⁷¹ to the halo height distribution mean (∼5.5 kpc) determined by Trotta et al. (2011); the exact value of the halo height is not critical for the analysis.

Zoom In Zoom Out Reset image size

For the IEM fitting procedure the GALPROP code is used to calculate all-sky γ-ray intensity maps from 1 to 100 GeV for 10 logarithmically spaced energy bins per decade for the Pulsars and OBstars baseline models, which are normalized to local CR data, using the configuration files for each available from the GALPROP website.⁷² The GALPROP code produces intensity maps in annuli that correspond to ranges in Galactocentric radii; the total intensity map for a given γ-ray production process (π⁰-decay, IC, Bremsstrahlung) is the sum of all the annular intensity maps for that process, and the total predicted γ-ray sky from a GALPROP run is the sum of intensities from all processes. Table 1 lists the Galactocentric annuli and the corresponding longitude ranges for the full extent of each annulus, as well as the "tangent" regions that are used in the fitting procedure for the components interior to the solar circle (see Appendix B of Ackermann et al. 2012a, for a description of the generation of the H i and CO gas annuli).

Table 1.  Galactocentric Annular Boundaries

Annulus	Rmin	Rmax	Longitude	Longitude
#	(kpc)	(kpc)	Range (Full)	Range (Tangent)
1	0	1.5	−10° ≤ l ≤ 10°	⋯
2	1.5	2.5	−17° ≤ l ≤ 17°	$10^\circ \,\leqslant | l| \,\leqslant 17^\circ$
3	2.5	3.5	−24° ≤ l ≤ 24°	$17^\circ \,\leqslant | l| \,\leqslant 24^\circ$
4	3.5	8.0	−70° ≤ l ≤ 70°	$24^\circ \,\leqslant | l| \,\leqslant 70^\circ$
5	8.0	10.0	−180 ≤ l ≤ 180°	⋯
6	10.0	50.0	−180 ≤ l ≤ 180°	⋯

Download table as:  ASCIITypeset image

The annular intensity maps are used as templates together with an isotropic component and a model for γ-ray emission associated with Loop I employing a two-component spatial template from Wolleben (2007) with a power-law spectral model for each, and point sources from the 3FGL source catalog.⁷³ This combined model is fit to the Fermi-LAT data excluding the 15° × 15° region about the GC using a maximum-likelihood method⁷⁴ , but with the point-source normalizations and spectral shapes held constant. Because they make only a small contribution this does not significantly affect the determination of the IEM parameters.

Two IEMs for each of the Pulsars and OBstars models—4 in total—are constructed. The two variants for each model are termed "intensity-scaled" and "index-scaled." The normalization parameters for the templates are determined in a series of fits to the data, starting at high latitudes for the local components and then working from the outer Galaxy to the inner Galaxy, always fixing the already determined normalization parameters in subsequent fits. For the intensity-scaled variants only the normalizations of the individual intensity maps are allowed to change. For the index-scaled variants the same fitting procedure is followed, but additional degrees of freedom are allowed to the spectrum of the gas-related interstellar emission when fitting to the annuli interior to the solar circle. The details of the procedure for the intensity-scaled variants are given in Appendix A. The motivation for the index-scaled variants is described further below.

Figure 2 upper and center panels show the fractional residuals, (data − model)/model, for 1–10 GeV energies⁷⁵ for the baseline and intensity-scaled Pulsars model. The isotropic component determined for the intensity-scaled IEM has been included in the baseline model for the fractional calculation to show the relative differences from the Galactic components of the IEMs. The regions not used in the fitting procedure are explicitly masked in the figure. They are not used because of localized extended excesses that are most likely unrelated to the large-scale interstellar emission. In particular, the band covering 70° ≤ l ≤ 90° includes the Cygnus region (l ∼ 75°–85°) around the Galactic plane; the corresponding band for negative longitudes is a consequence of the axisymmetric nature of the model being used. The range 90° ≤ l ≤ 270° is used to constrain the IC emission from annulus 5 so the data out of the plane from the 70° ≤ l ≤ 90° region is not required to constrain this component. The $-20^\circ \,\leqslant \,l\,\leqslant 20^\circ ,10^\circ \,\leqslant | b| \,\leqslant 50^\circ$ region where the Fermi haze/bubbles have been detected is also excluded. Including these regions in the fitting procedure would bias the normalization of the IEM components because models for these features are not included in this study.

Zoom In Zoom Out Reset image size

Outside of the Galactic plane the fractional residuals are substantially reduced for the intensity-scaled IEM compared to the baseline. This is due to the scaling of the π⁰-decay interstellar emission for the local annulus, and of the IC component generally. The coefficient for the local gas (annulus 5) interstellar emission is adjusted upward (see Table 5 in Appendix A for a full list of IEM coefficients). Meanwhile, the coefficients of the IC intensity for the local and other annuli interior to the solar circle are increased compared to the baseline IEM. Ackermann et al. (2012a) also found that better fits to the γ-ray data were generally obtained by increasing the IC intensity for large regions of the sky for the baseline IEMs. While the GALPROP version used by Ackermann et al. (2012a) only allowed the calculation of all-sky IC intensity maps, the decomposition of the IC intensity into Galactocentric annuli used here enables mismodeling of the IC emission due to uncertainties of the gradients in the CR electron and ISRF distributions to be more accurately treated (Section 4.1).

Along the Galactic plane the ≳30% under-prediction by the baseline model is reduced to ≲±10% after this scaling, except for scattered regions. The longitude range l ∼ −(15–70)^◦ is the largest such region where the intensity-scaled IEM over-predicts the data by ∼20%–30% in the 1–3.16 GeV band. Within ∼10° of the mid-plane this may indicate that the spectrum of the IEM related to the CR nuclei/gas interaction (π⁰-decay) in this region is too soft. To account for this, additional degrees of freedom to the spectrum⁷⁶ of the π⁰-decay interstellar emission are allowed, and the model is refitted for annuli interior to the solar circle following the same sequence of regions as for the intensity-scaled IEMs—this is the "index-scaled" IEM variant.

Figure 2 lower panels show the fractional residuals for the index-scaled variants of the Pulsars IEM. The fractional difference for the l ∼ −(15–70)^◦ region in the 1–3.16 GeV energy range is reduced to ≲±10% with a slight increase in the residual for the corresponding positive longitude range. At mid-to-high latitudes the residual is reduced because the IC for the annuli interior to the solar circle and Loop-I model are also refit.

The intensity-scaled IEM gives a lower residual for positive longitudes inside the solar circle around 1 GeV, while the index-scaled variant is closer to the data at negative longitudes. The converse appears at higher energies where instead the intensity-scaled IEM gives lower residuals at negative longitudes than the index-scaled variant. The underlying axisymmetric geometry for the IEMs is partly responsible for this: there is not enough freedom in the model parameters, even for the index-scaled variant, to account for differences due to any error in the assumption of an average CR distribution for each Galactocentric annulus (the limitations of the IEM tuning procedure and resulting models are discussed further in Section 5.4). Qualitatively, similar results for the scaled OBstars IEMs (not shown) are obtained.

It is not straightforward to identify a best IEM after fitting because the qualitative improvement for each over the corresponding baseline IEM is similar. Consequently, all 4 (Pulsars/OBstars, intensity-/index-scaled) IEMs are used to estimate the fore-/background toward the 15° × 15° region about the GC below.

3.2.1. Point-source Candidates

3.2.2. Combined Interstellar Emission and Point Source Fit

3.2.3. Residual Maps and Iteration

A potential drawback of using the wavelet detection algorithm is that fainter point sources may be missed with a single iteration. This is remedied in this analysis by iterating the point-source seed detection on the residual maps following the maximum-likelihood fit and rerunning the analysis chain. The iteration is made until there are no new significant excesses in the TS map of the region that satisfy the signal-to-noise criterion adopted in the PGWave detection step. The TS map is determined by moving a putative point source with a PL with spectral index −2 using PointLike through a grid of locations in the region and by maximising the likelihood function at each location. The positions of peaks with TS > 9 are added to the source model.

Figure 3 (left panel) shows the TS map following the gtlike maximum-likelihood fit using these candidates in the first full iteration for the Pulsars intensity-scaled IEM. Some significant excesses remain after the initial pass, in particular around the GC. The PGWave and PointLike seed-finding and optimization steps are repeated (Section 3.2.1), finding 37 additional candidates with a PointLike-assigned TS > 9 for the Pulsars intensity-scaled IEM. The left bottom panel in Figure 3 shows the TS map after the second full iteration of the analysis of the 15° × 15° region for this IEM. Some excesses with TS > 25 remain. However, these correspond to point-source candidates reported by PGWave in the first iteration that were rejected because their gtlike-assigned TS <9, possibly due to structures having angular extensions that are larger than the PSF core.

Zoom In Zoom Out Reset image size

The procedure for the initial optimization, binned likelihood fit, and iteration is also done using the OBstars intensity-scaled IEM. This yields 85 point source candidates with PointLike-assigned TS > 9. The point-source candidates are combined with their PointLike trial spectra together with the OBstars intensity-scaled IEM in a maximum-likelihood fit using gtlike, following the same procedure as outlined above. The TS map after the gtlike maximum-likelihood fit for this model is shown in the bottom right panel of Figure 3. The results are very similar to the second iteration of the analysis using the Pulsars intensity-scaled model. Consequently, no further iteration is made for the OBstars model.

For the Pulsars and OBstars index-scaled variants, the point-source seeds with TS > 9 from the corresponding intensity-scaled IEM are used in the source model for the maximum-likelihood fit. This procedure yields TS maps that are very similar to the respective intensity-scaled counterparts. The results are summarized in Section 4.

Figure 4 shows the differential spectra of the individual components obtained for the 4 IEMs integrated over the 15° × 15° region about the GC. The figure separates the emission components in terms of the contributions by π⁰-decay and IC for annulus 1, the interstellar emission fore-/background, and point sources over the region. As expected, the fore-/background dominates for each IEM, which is predominantly π⁰-decay in origin. IC scattering is the dominant interstellar emission component over the inner ∼1 kpc. This contrasts with the predictions of the baseline Pulsars and OBstars IEMs, in which the neutral gas π⁰-decay interstellar emission components over the same region have similar fluxes. Combined, the GALPROP-predicted H i and CO-related π⁰-decay annulus 1 templates are brighter by up to an order of magnitude than the IC emission for either model. The fit for the 15° × 15° region preferentially adjusts the annulus 1 IC component while suppressing the H i and CO-related π⁰-decay templates for all IEMs. The scaling factors for the annulus 1 IC templates are ∼6–30, with higher values for the OBstars IEM variants. While the difference in the scaling factors between the IEMs is large, the final flux determined for the IC over annulus 1 over all four IEMs is within a factor ∼1.5.

Zoom In Zoom Out Reset image size

Figure 5 shows the spatial distribution of the individual components over the 15° × 15° region about the GC for energies >1 GeV for each IEM. As noted above, the fore-/background interstellar emission is clearly the brightest component. Although the intensity scale does not directly show a strong variation within ∼±1° of the mid-plane the fore-/background over the four IEMs varies by ∼30%, which is similar to elsewhere in the plane. Outside of the plane the brightness of individual features varies. For example, around l ∼ 0°–2°, b ∼ 3° the Pulsars IEMs are dimmer. Other low-intensity structures appear to change subtly for the different IEMs. This is predominantly due to scaling of the π⁰-decay emission from annuli 2–4, which is discussed further below.

Zoom In Zoom Out Reset image size

The IC intensity is more peaked toward the GC for the OBstars model compared to the Pulsars IEM. This is due to spatial distributions of CR sources employed for the respective IEMs, which result in different spatial distributions for the propagated CR electron intensities over the inner few kpcs. The paucity of sources within a few kpc of the GC for the OBstars IEM gives a constant intensity because the electrons must diffuse from larger Galactocentric radii. On the other hand, the Pulsars source distribution peaks around a few kpc and is non-zero in toward the GC, giving a higher electron intensity with a gradient. Combined with the spatial distribution over the same region of the ISRF intensity, which peaks at the GC, the results are IC templates that are broader (Pulsars) or more peaked (OBstars) for annulus 1.

Table 2 gives the fluxes of the different components integrated over the four energy bins from 1 to 100 GeV. Trends that can explain the low level of emission associated with the annulus 1 π⁰-decay component are not readily apparent. Because the fore-/background is held constant for the maximum-likelihood fit over the 15° × 15° region (Section 3.2.2) there is no correlation matrix with the IEM that can be examined to determine degeneracies for the structured interstellar emission component.

But some understanding of the effect of the fore-/background can be inferred from examining the fluxes per annuli for each of the IEMs (see Table 6 in Appendix A). The IEM scaling procedure (Section 3.1), gives similar contributions by IC emission for annuli >1 to the total flux over the 15° × 15° region for the Pulsars and OBstars IEMs. Also, for each IEM the total emission from the scaled local and outer annulus (annuli 5 and 6) across the 15° × 15° region is very similar.

The major difference is the distribution of the π⁰-decay flux over annuli 2–4. The standout feature is how the H i-related π⁰-decay flux for annuli 3 and 4 (Table 6) is correlated with the essentially complete suppression of the π⁰-decay flux in annulus 1 for the OBstars index-scaled IEM. Similar patterns appear with the other IEMs, but are less pronounced. The scaling procedure results in combinations of the structured fore-/background emission that leave only a small amount of flux for the annulus 1 π⁰-decay templates to be assigned by the likelihood maximization (Section 3.2.2).

While the IC and combined point source fluxes are larger than the formal statistical uncertainties, and the variation across IEMs, this is not true for the annulus 1 π⁰-decay emission. For this component the reliability of the fluxes obtained is uncertain: the intrinsic variation of the structured fore-/background across IEMs is large in comparison, and there is also a potential correlation with some of the flux attributed to the point sources, which is discussed further below.

Table 3 summarizes the properties of the 48 point sources ordered by increasing right ascension over the 15° × 15° region with TS ≥ 25 for the Pulsars intensity-scaled IEM, which was the model used for the point-source positions and localization uncertainties. The list is termed the first Fermi Inner Galaxy point source catalog (1FIG). 27 of the 1FIG sources have a 95% containment localization error ellipse that intersects the 95% containment radius of a 3FGL point source. These associations are also given in the table. The correlation plot (Figure 6) shows that the fluxes of the 1FIG sources compared with their associations in the 3FGL are in good agreement.

Zoom In Zoom Out Reset image size

Fourteen sources in the 3FGL have a multi-wavelength association. For example, the 3FGL source J1701.2-3006 is associated with the globular cluster NGC 6266. For the given TS ≥ 25 detection threshold used for 1FIG source detection 10 counterparts in the 3FGL are obtained that have a multi-wavelength association. The W28 supernova remnant (3FGL J1801.3-2326e) is included in the model of the region as an extended source (see Section 3.2.1), while the 3FGL sources J1716.6-2812 (NGC 6316), J1746.3-2851c (PWN G0.13-0.11), and J1750.2-3704 (Terzan 5) are missing from 1FIG. The latter missing counterparts are discussed in Section 5.2, together with possible multi-wavelength associations for other 1FIG sources.

Figure 7 shows the point sources from the 1FIG and 3FGL overlaid on the total photon counts for the 15° × 15° region about the GC. The 3FGL sources are separated according to whether they have an analysis flag set in the 3FGL catalog: flagged sources indicate their properties depend on the IEM or other details of the analysis in the region. The density of flagged 3FGL sources is higher out of the Galactic plane than that of the 1FIG sources, even if the ${TS}\lt 25$ source candidates are included. This can be partly attributed to differences in the treatment of the IC emission, and its interplay with the gas components, for the IEMs employed for the respective analyses. The 3FGL IEM uses an all-sky IC map based on a GALPROP calculation, with its spatial distribution taken as a fixed template and the spectral parameters adjusted to improve the correspondence with the data (Acero et al. 2015a). The decomposition of the IC intensity map into Galactocentric annulus templates employed here for the first time introduces additional degrees of freedom that can account for Galactocentric radial gradients in the IC emissivities. This allows more flexibility to fit for a spatial distribution of IC emission that is not correctly represented by the baseline GALPROP calculations.

Zoom In Zoom Out Reset image size

The 3FGL source density is higher out of the Galactic plane compared to the 1FIG, while the reverse is the case closer to the plane. The 1FIG sources in the plane do cluster in approximately the same regions as the high-density clusters for the 3FGL: near the W28 supernova remnant and the GC. Outside these regions the density of sources is higher than the 3FGL. That many of these additional sources appear to trace features in at least one of the annulus 1 templates is suggestive that they may be misattributed interstellar emission. This can be seen in Figure 8 where the 1FIG sources are overlaid on the fitted components of the interstellar emission for annulus 1 for the Pulsars intensity-scaled IEM. By way of example, many of the TS > 25 sources appear to trace the edge of the fitted neutral gas π⁰-decay template. Because many of these sources lack a multi-wavelength association it is not straightforward to determine whether they are true point sources.

Zoom In Zoom Out Reset image size

The combined flux of 1FIG point sources and point source candidates across the 15° × 15° region for the Pulsars intensity-scaled IEM > 1 GeV is 44.6 ± 1.4 × 10⁻⁸ ph cm⁻² s⁻¹. Of this total, only 20% is due to point sources with a multi-wavelength association in the 3FGL, while 10% of the total is due to TS < 25 source candidates. Focussing on the region −75 ≤ l ≤ −05, − 15 ≤ b ≤ 15 (the region with the highest density of 1FIG sources without 3FGL counterparts and where they appear to trace the edge of the π⁰-decay template), the combined flux from the TS > 25 1FIG sources is ∼20% of the total point-source flux over the 15° × 15° region. The total flux for the annulus 1 π⁰-decay along the entire plane from the fit is about the same as that from these sources alone.

Although the absolute values differ similar relative contributions to the total point source flux determined for each of the other IEMs are obtained. However, for the Pulsars index-scaled, and both OBstars IEM variants, the annulus 1 π⁰-decay template is even less intense than for the Pulsars intensity-scaled IEM. Over all IEMs, the 17 1FIG sources with TS > 100 have a variation in the combined flux that is ≲1%. For sources without a 3FGL multi-wavelength association and with 1FIG TS in the range 25 ≤ TS < 100 (24 sources) the combined flux over this region is more strongly dependent on the IEM: 6.7–8.3 × 10⁻⁸ ph cm⁻² s⁻¹. The larger of these values is comparable to the variation to the total fore-/background for the IEMs over this region (Table 2).

It is probable that there is some misattribution of interstellar emission to low-flux (e.g., less than ∼few $\times {10}^{-9}$ ph cm⁻² s⁻¹ > 1 GeV) point sources. The low-flux sources are all relatively low-significance sources and modeled using power-law spectra (Section 3.2.1). The distribution of their spectral indices over the 15° × 15° region may provide some information: softer spectral indices (e.g., ≳2.5 in spectral index) can indicate that the low-flux sources are more likely associated with the structured/gas-related interstellar emission, while harder indices can indicate a more "IC-like" distribution. Figure 9 shows all point sources and candidates with a TS < 50 overlaid on the fitted π⁰-decay annulus 1 template for the Pulsars intensity-scaled IEM. The point sources are coded according to the spectral indices: circles show those with indices >2.5, while triangles show those with indices ≤2.5. There is no clear trend of softer spectrum point sources tracing the structured emission, nor one where the harder spectrum point sources have a high density out of the plane. It is difficult to identify the exact fraction of the emission, or to what component (gas-related, IC), the low-flux point sources could be ascribed if they are indeed due to mismodeling of the interstellar emission over the region.

Zoom In Zoom Out Reset image size

Figure 4 shows the fractional residuals below the differential flux spectrum for each IEM integrated over the 15° × 15° region as a function of energy. Some trends are evident: each model over-predicts the data below ∼2 GeV and under-predicts above ∼2 GeV, except for the Pulsars index-scaled IEM that over-predicts the data ≳5 GeV.

Figure 10 shows the longitude and latitude profiles for the energy ranges⁸⁰ 1–1.6, 1.6–10, and >10 GeV for the Pulsars index-scaled model shown in Figure 4, which has the lowest fractional residual across the 1–100 GeV energy range. (The features are mostly the same for the other IEMs with the major difference their magnitude in terms of counts, hence these profiles are not shown because of their similarity.)

Zoom In Zoom Out Reset image size

The lower sub-panel for each figure gives the residual counts (data − model). While there is considerable statistical noise, the total residual counts may be distributed asymmetrically in longitude about the GC below 10 GeV. However, quantifying such an asymmetry using a purely data-driven method, e.g., by forming the ratio $A\,=\,({f}_{+}-{f}_{-})/({f}_{+}+{f}_{-})$ where f₊ and ${f}_{-}$ are the counts for some equally sized regions about some symmetry line, is not useful here because the residuals are of mixed sign.

Figure 11 shows in greater detail the spatial distributions of the residual for each of the IEMs and in the three energy bands. Common features across IEMs are present: the model is too bright compared to the data mostly along the Galactic plane for the lowest energy band (1–1.6 GeV) and this behavior is more pronounced for the intensity-scaled IEMs, while the models under-predict the data around the GC in the 1.6–10 GeV energy band.

Zoom In Zoom Out Reset image size

Note that the ecliptic crosses the 15° × 15° region and therefore the Sun and Moon contribute to the observed emission Abdo et al. (2011). The emissions from these objects are not included in the fore-/background models employed in this analysis. But above 1 GeV it is small relative to the observed residual emission.⁸¹

Although the spatial distribution of the residuals is not suggestive of a contribution by the Fermi bubbles, it is also possible that there is some emission from them over the 15° × 15° region. Without a spatial template for the Fermi bubbles over the region their contribution is tested using a model with an isotropic spatial distribution across the 15° × 15° region with intensity and spectrum as determined from analyses at higher latitudes (e.g., Ackermann et al. 2014). The data − model agreement only marginally improves if this contribution is included.

The model over-prediction at the lowest energies is primarily correlated with the Galactic plane, which could be due to mismodeling of the gas component of the IEMs. Some of the positive residual in the few GeV range could be due to an extended component that is more concentrated toward the GC compared to the IEM components. But, the profiles shown in Figure 10 and the spatial distributions shown in Figure 11 represent the situation if the fit is made only for interstellar emission about the GC and point sources. As a consequence it is difficult to establish properties for an additional component not presently included in the model for the region. A spatial and spectral model needs to be assumed and fit to the data together with the interstellar emission and point sources.

A set of templates for the spatial distribution of the additional component is selected with each fit together with the interstellar emission components and point sources using the maximum-likelihood procedure described above. Because the excess emission in the few GeV range is distributed around the GC, templates that peak there are considered. A set of two-dimensional Gaussians with varying HWHM (1°, 2°, 5°, 10°) are used. While these spatial distributions do not have an obvious physical interpretation, they can be used to gauge the radial extent of the positive residuals. Spatial templates to model the predicted distribution for γ-rays produced by DM particles annihilating or decaying in the Galaxy are included. The inner region of the Galaxy is predicted to be the brightest site for a DM signal in γ-rays and could be well within the sensitivity of Fermi-LAT. To model the DM density distribution, the Navarro, Frenk, and White (NFW) (Navarro et al. 1997) profile is employed⁸² with different choices for the slope of the profile in the innermost region, γ = 1, 1.2. The NFW profile is predicted by simulations of cold DM, while the more peaked distribution with γ = 1.2 (NFW-c) is motivated by earlier work (e.g., Hooper & Goodenough 2011; Abazajian & Kaplinghat 2012) and could arise when baryonic effects are included in simulations. The square of the NFW profile is used as a template for DM annihilation, hereafter referred to as simply the "NFW profile." The possibility that an unresolved population of γ-ray point sources such as pulsars distributed along the Galactic plane is contributing to the observed emission is also considered. Predictions of the γ-ray emission from unresolved pulsars exist and they span a range of possibilities (e.g., Story et al. 2007; Faucher-Giguère & Loeb 2010). Here, the spatial distribution of an unresolved pulsar population is modeled using the distribution of the CO gas in annulus 1, because this is a likely tracer for regions of high-mass star formation, and smooth it with a 2° Gaussian to account for the scale height of the associated pulsar population.

For each of the spatial templates listed above, the spectrum is modeled with an exponential cut-off power law. This form has some flexibility to model a pulsar or a DM annihilation spectrum without supposing specific scenarios. For each of the spatial templates listed above and for each of the IEMs, a maximum-likelihood fit is made in the 15° × 15° region as described in Section 3.2.2.

The improvement in likelihood as well as the resulting best-fit parameters for the spectrum of the additional component are summarized in Table 4.⁸³ All templates yield statistically significant improvements compared to the model without the additional component. The largest improvements are observed for the NFW annihilation templates, whereas the unresolved source component yields the smallest improvements.

The new component spectra present harder spectral indices and lower energy cutoffs for the index-scaled IEMs compared to the intensity-scaled variants. This is consistent with the index-scaled models having overall better agreement with the data at higher energy, and therefore attributing the positive residual found for the intensity-scaled IEMs ≳ 10 GeV to gas related emission rather than to the new component. Within the same IEM, the spectrum for the more peaked templates (NFW and NFW-c for DM annihilation, and the 1° Gaussian) present softer indices and higher energy cutoffs. The NFW decay and the 10° Gaussian (the more extended templates) perform similarly to each other for most IEMs.

Among the Gaussian templates, the 2° and 5° gaussians perform better for the Pulsar IEMs, while the 5° and 10° gaussians for the OB stars IEMs. This result is an indication that the Gaussian templates might be compensating for mismodeling of the IC contribution, whose morphology differs for the OB stars and Pulsars IEMs.

By including the NFW profile component the agreement with the data has an overall improvement for all the models up to ∼30 GeV, as shown in Figure 12, with the Pulsars index-scaled variant yielding the best agreement over the full energy range. However, a broad range for the best-fit parameters of the spectral model is found. The variation is not easily ascribed to a covariance with only a single component of the model that is fitted over the 15° × 15° region. For example, the annulus 1 IC and H i-related π⁰-decay normalizations adjust in the fit to compensate for the additional template. But the spectral parameters of the residual template are not solely determined by the fit with the interstellar emission components and point sources over the inner region about the GC; the fore-/background interstellar emission has an effect as well.

Zoom In Zoom Out Reset image size

The intensity-scaled IEMs yield similar spectral parameters for the NFW template, but the results for the index-scaled IEMs have a stronger variation. This can be seen in Figure 13, which shows the flux spectral envelopes from including the uncertainties on the normalization and spectral index obtained for the NFW template for the 4 IEMs. The index-scaled IEMs have the distinction of harder spectra for the π⁰-decay interstellar emission for annuli 2–4 (Table 5), but also modified IC contributions for annuli 2 and 3 compared to their intensity-scaled counterparts. The majority of the π⁰-decay fore-/background interstellar emission is due to annulus 4 and, as already noted in Section 4.1, even small variations in the structured fore-/background interstellar emission can have a follow-on effect on the spatial distribution of the residual emission over the 15° × 15° region. It is difficult to test how small variations in the π⁰-decay fore-/background from this annulus affect the residual model parameters because the annulus 4 fit parameters are determined at an intermediate step in the fitting. But the comparison between the results for the Pulsars and OBstars index-scaled IEMs show that the different spectral parameters obtained for the structured interstellar emission fore-/background can alter the final fitted values for all components over the 15° × 15° region.⁸⁴

Zoom In Zoom Out Reset image size

With the interstellar emission fore-/background held constant for each IEM, the interplay between the centrally peaked positive residual template and the interstellar emission components is not surprising. Because the IC component is maximally peaked toward the GC for all IEMs an additional template that is also peaked there will also be attributed some flux when fit. Over all IEMs the effect of including the NFW model for the residual results in an IC annulus 1 contribution that is up to three times smaller and H i annulus 1 contribution that is up to three times larger.

Note that even if a centrally peaked template is included as a model for the positive residual, it does not account for all of the emission. This can be seen in Figure 14, which shows the residual counts for the NFW template and IEM with the best spectral residuals (Pulsars index-scaled). Qualitatively, the remainder does not appear distributed symmetrically about the GC below 10 GeV, and still has extended positive residuals even at higher energies along and about the plane.

Zoom In Zoom Out Reset image size

This study is the first using the Fermi-LAT data that has made a separation between the large-scale interstellar emission of the Galaxy and that from the inner ∼1 kpc about the GC. The IC emission from annulus 1 is found to dominate the interstellar emission from the innermost region, and represents the majority of the IC brightness from this component along and through the line of sight toward the GC. The contribution by the IC from annulus 1 to the total flux depends on the IEM and whether the residual is fitted (Section 4.3). For the latter case the IC from annulus 1 is still up-scaled compared to the GALPROP predictions, but by a factor ∼2 lower than if fitted solely for the interstellar emission components and point sources. The remainder is distributed across the H i-related π⁰-decay annulus 1 component and the template used to fit the residual centered on the GC. For either case (residual template used/not-used), the fitted fluxes attributed to the IC annulus 1 component across all IEMs are within a factor ∼2—the flux and its range is the important quantity, instead of the individual (model-dependent) scaling factors.

The Pulsars intensity-scaled IEM with the residual template gives the minimal "enhanced" flux for IC annulus 1. The average CR electron intensity ≳5 GeV in the Galactic plane is estimated for this model within ∼1 kpc of the GC as ∼2.8 ± 0.1 × 10⁻⁴ cm⁻² s⁻¹ sr⁻¹, where the uncertainty is statistical only. This energy range is used because its lower bound corresponds to the CR electron energies producing ∼1 GeV IC γ-rays. This is ∼ a factor of two higher than the local total CR electron density for this same energy range for the Pulsars baseline model. On the other hand, the OBstars intensity-scaled IEM fitted without the residual component gives the maximal "enhanced" flux for IC annulus 1. The average CR electron intensity ≳5 GeV in the Galactic plane within ∼1 kpc of the GC for this IEM is ∼9.4 ± 0.1 × 10⁻⁴ cm⁻² s⁻¹ sr⁻¹.

Measurements of the interstellar emission at hard X-ray energies to MeV γ-rays by INTEGRAL/SPI (Bouchet et al. 2011) show that the majority is due to IC scattering by ∼GeV energy CR electrons off the infrared component of the ISRF.⁸⁵ The GALPROP calculations, which follow the same "conventional" model normalization condition to local CR measurements as used in this paper, made to interpret the SPI measurements indicate that IEMs with at least factor of two higher CR densities toward the inner Galaxy are a plausible explanation for the data. Another possible explanation is a higher intensity for the radiation field energy density in the inner Galaxy than used in the standard ISRF model of Porter et al. (2008); these possibilities are not tested here because they require detailed investigations that are beyond the scope of the current work. The higher CR electron densities obtained from this analysis are plausible given the same electrons are IC scattering different components of the ISRF to produce the interstellar emission ≳1 GeV and at SPI energies.

The purpose for fitting the baseline IEMs to the data was to obtain estimates for the interstellar emission fore-/background. However, the fit results for the individual rings for each IEM potentially give some information on the large-scale distribution of CRs througout the Galaxy. Tables 5 and 6 in Appendix A.1 give the fit coefficients and fluxes for the scaled IEMs, while Figure 15 shows the integrated fluxes for the 1–10 (top) and 10–100 GeV (bottom) energy ranges, respectively, over the 15° × 15° region for the GALPROP-predicted and scaled version of each IEM for the Pulsars (left) and OBstars (right) source distributions.

Zoom In Zoom Out Reset image size

The fitting procedure generally increases the intensity of each annulus relative to the nominal model. The coefficients for the intensity-scaled Pulsars and OBstars IEMs are mostly higher than the GALPROP predictions toward the inner Galaxy (annuli 2–3). Those for the OBstars IEM are higher than the Pulsars, which reflects the fact that the spatial distribution for the CR sources in this model cuts off within ∼2 kpc of the GC. The cut off in the OBstars source spatial distribution produces a predicted CR intensity that is lower compared to the Pulsars IEM over this region. The fitting procedure adjusts the OBstars predictions upward more than the Pulsars to compensate. This indicates that a Pulsars-style spatial source distribution is closer to the real spatial distribution of sources within ∼2 kpc of the GC. But, even the Pulsars spatial source distribution is scaled up by the fit over this region, indicating that even more "peaked" source models, or some modification to the propagation model, is required to describe the distribution of CRs toward the inner Galaxy. Meanwhile, there is more similarity in the scaling coefficients for annuli 4–6. This reflects that the CR source distributions and propagation conditions for both IEMs are not significantly different in their Galactocentric radial distributions in these annuli.

The spectral parameters for the annuli interior to the solar circle for the index-scaled variants give results that are strongly dependent on the IEM being fit. For the Pulsars IEM the spectrum of the CR nuclei/gas interstellar emission is consistently harder across annuli 2–4 for both CO and H i components than the intensity-scaled IEMs. For the OBstars IEM only the H i component has a hardening to the spectrum across annuli 2–4. For this IEM the fits for annuli 2–3 were unstable when fitting both CO and H i components. Because the size of the regions are small, the low flux of the annuli 2–3 components in comparison to those that are already-determined from fitting to the outer longitude ranges means that the data are insufficiently constraining. However, a convergent fit is obtained if the CO-related π⁰-decay templates is set to the GALPROP prediction. The motivation for allowing the additional freedom to fit the spectrum for the gas-related interstellar emission interior to the solar circle is solely to improve the fit residuals. But, the harder index for the H i and CO component when fitting the Pulsars IEM can be an indication that the assumption of a uniform CR source spectrum across the Galaxy is insufficient, or that the diffusive propagation of CRs is non-uniform.

Generally, the fitting results can be interpreted as a reconfirmation that the CR gradient in the Galaxy is flatter than expected based on current knowledge of the Galactocentric radial distribution of CR sources, which has been known since the SAS-2 (Stecker & Jones 1977), COS-B (Bloemen et al. 1986; Strong et al. 1988), and EGRET (Hunter et al. 1997; Digel et al. 2001) all-sky surveys. The explanation is not clear. Bloemen et al. (1993) suggested that the radial distribution of CR sources derived from observations may be biased and their real distribution is flatter or the diffusion parameters derived from the local CR measurements are not the same throughout the Galaxy. Solutions to this issue in terms of CR propagation phenomenology have been proposed: CR–driven Galactic winds and anisotropic diffusion (Breitschwerdt et al. 2002), or non-uniform diffusion coefficient that increases with Galactocentric radius and the distance from the Galactic plane (Shibata et al. 2007).

The current analysis has focussed on finding IEMs to estimate the fore-/background toward the inner Galaxy. The broader implications of the scaled IEMs developed in this study for the large-scale distribution of CRs in the Galaxy are deferred to future work.

Figure 16 shows the 1FIG sources and source candidates overlaid on the Fermi-LAT data used in this paper, and 3FGL multi-wavelength associated sources, together with SNRs from Green's SNR catalog⁸⁶ (Green 2014) and pulsars from the ATNF catalog⁸⁷ (Manchester et al. 2005), respectively, that are within 95% of the 1FIG source/source candidate error ellipse. The 3FGL sources that have likely counterparts at other wavelengths that are listed in the catalog not detected in the 1FIG are either due to a too low TS (3FGL J1716.6-2812—NGC 6316), or are more than the 95% containment radius of the error ellipse from a potential 1FIG counterpart (3FGL J1750.2-3704—Terzan 5 and 3FGL J1746.3-2851c—PWN G0.13-0.11).

Zoom In Zoom Out Reset image size

There are 14 1FIG sources and source candidates with overlaps with the above mentioned SNR and pulsar catalogs. Multiple overlaps occur across and within the catalogs, e.g., SNR 354.1+00.1 and PSR J1701–3006A,D,E overlap with 1FIG J1701.1–3004.

The 1FIG source J1801.4–2330 overlaps with SNR 006.5–00.4, which has been detected in the first LAT catalog of SNRs (Ackermann et al. 2016). 1FIG J1740.1–3057 overlaps with SNR 357.7–00.1 (MSH 17–39), and has been detected previously in Fermi-LAT data (Castro et al. 2013). The source 1FIG J1745.5–2859 overlaps with SNR 000.0+00.0 (Sgr A East), but this is in a strongly confused region and other counterparts may be possible. The 3 other SNRs (SNR 354.1+00.1, SNR 355.4+00.7, and SNR 000.3+00.0 corresponding to 1FIG J1730.2-3351, 1FIG J1731.3–3235, and 1FIG J1746.4–2843, respectively) are new detections in high-energy γ-rays at Fermi energies. Follow-on studies are required to better characterize their spatial and spectral properties.

The comparison with the ATNF catalog yields 9 1FIG sources overlapping with known pulsars. The 1FIG source J1750.2–3705 is the counterpart of the globular cluster NGC 6441, which has been detected in high-energy γ-rays (Tam et al. 2011). Four of the remaining eight overlap with nearby pulsars (≲0.2 kpc) and are listed in the LAT Second Catalog of Gamma-ray Pulsars⁸⁸ (Abdo et al. 2013). The remaining four sources have been identified previously and searches for pulsed emission have been made but with no detections (Abdo et al. 2013).

Obviously the comparison made here between the 1FIG sources and γ-ray source classes is not exhaustive. However, more than two thirds of the 1FIG sources do not have associations with sources in known classes of γ-ray emitters. The unassociated 1FIG sources tend to be close to the Galactic plane. It remains a significant possiblity that a majority of the point sources found over the 15° × 15° region can be attributed to mis-identified interstellar emission, as already discussed in Section 4.2.

A number of studies of the residual emission toward the inner Galaxy have been performed, as described in Section 1. Figure 13 compares the results from this analysis with selected results from the literature. The comparison is useful because a similar spectral model to other authors is used when fitting the residual emission associated with the centrally peaked spatial templates. The IEMs that are used in this paper are scaled to the data outside of the 15° × 15° region to reduce the discrepancies, particularly along the Galactic plane, by the a-priori GALPROP-generated IEMs from Ackermann et al. (2012a). The developments also made as part of this work have allowed additional degrees of freedom to be included for the scaling of the IEMs—treating the IC template as the sum of individual templates with the same Galactocentric radial boundaries as the π⁰-decay templates, and allowing for spectral variations in the π⁰-decay templates from those predicted interior to the solar circle—that go beyond the modeling of the interstellar emission employed by other analyses (e.g., Calore et al. 2015b). In addition, a catalog of point sources is derived for each IEM that is used in the analysis of the inner region about the GC. The prescriptive method of determining the fore-/background interstellar emission, together with the self-consistent treatment of the point sourcess⁸⁹ and interstellar emission for the inner ∼kpc about the GC allows the least biased estimate to-date to be made of the positive residual emission about the GC. This work finds that for individual IEMs the spectral parameters for a spatial template that peaks at the GC, such as the NFW profile, can be relatively tightly constrained. However, over all IEMs considered in this work the variation of, for example, the cut-off energy for an exponential power-law spectral model is much wider than that for any individual model.

Although the spectral residuals are generally improved by an additional template, discrepancies remain that are more pronounced for the intensity-scaled variants of the IEMs. It should be emphasized that despite this observation the intensity-scaled IEMs cannot be excluded on the basis of fits made to the 15° × 15° region about the GC. All four of the IEMs are tuned to data outside this relatively small region, providing similar improvements to the all-sky residuals, and hence are equivalent representations of the fore-/background toward and through the GC. Because of the limitations in modeling the interstellar emission, the higher energy cutoff spectra for the NFW profile component with the intensity-scaled IEMs cannot be ruled out. With the limited freedom for the interstellar emission components (only the normalization for the IC, H i and CO-related π⁰-decay intensity maps are allowed to vary—three parameters) and in the specification of the spectral model for the positive residual (normalization, spectral index, and cut-off energy—three parameters) the spread in the positive residual template parameters is considerable (see above.) If more freedom is allowed for the spectrum of the positive residual then the spectral residuals for all four IEMs over the 15° × 15° region can be very small. Figures 17 and 18 show the results for all IEMs if more degrees of freedom are allowed to model the spectrum of the NFW profile⁹⁰ (note that the feature at ∼40 GeV is not significant when the fit uncertainties are considered, as shown in Figure 18). For this choice of spectral model indeed the residuals are very good for all IEMs. It is therefore premature, because of the variations in the IEMs and their limitations, to favor a specific IEM among those we considered and to attribute the high energy residual to a particular origin.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although a large formal statistical significance may be indicated for the detection of a new component, note that fitting a centrally peaked profile does not account for all of the positive residual over the 15° × 15° region. Ascribing a singular origin to such a residual component is premature given the limited constraints on the other emission components over the 15° × 15° region. A complete assessment of the uncertainties (see Section 5.4) is required to understand the nature of its spatial and spectral parameters. The current work demonstrates that even in the optimistic scenario where the presence of a DM component in the data might be established based on the spatial distribution of the associated γ-ray emission, important information on the DM particle such as its mass and annihilation spectrum is strongly dependent on the IEM. This was first demonstrated by Agrawal et al. (2015) using preliminary results based on this work.

The IEMs used in this analysis are cylindrically axisymmetric averaging the CR densities and other details azimuthally about the GC. Some shortcomings of these models were noted earlier (Section 3.1), particularly when fitting for the large-scale interstellar emission interior to the solar circle.⁹¹ The density of CR sources and other ISM components is highest inside the solar circle and is most likely radially and azimuthally dependent, e.g., being associated with spiral arms or the Galactic bar/bulge. However, in the absence of detailed three-dimensional models for the interstellar gas, radiation field, and CR sources the axisymmetric models are the only viable method for estimating the fore-/background toward and through the line of sight to the GC.

The IEM fitting interior to the solar circle uses the tangent ranges for positive and negative longitudes to obtain parameters for the annuli 2–4 (Table 5). To examine the effect of the azimuthal averaging, fits to the tangent ranges were made for positive and negative longitudes to gauge the difference in the parameters for the IEMs obtained when considering each separately. The scaling factors for annulus 4 obtained when fitting negative and positive longitude ranges were statistically consistent⁹² with those found when fitting both ranges combined. For annuli 2 and 3 the fits to the positive and negative tangent longitude ranges result in scaling parameters that differ by factors up to ∼2 from each other, which is well beyond the statistical uncertainty; the average value obtained by fitting both tangent ranges together is approximately in-between for the intensity-scaled IEMs over annuli 2 and 3. For the index-scaled IEMs the spectral parameters are harder or softer than the average when using the positive/negative tangent ranges individually for annuli 2–4. However, there is no clear trend and the over/under-prediction is not confined to a particular energy interval.

The uncertainty for the IEM fore-/background flux toward the GC due to the azimuthally averaged IEMs is difficult to quantify precisely. A minimal estimate can be made from the statistical uncertainty for the annulus 4 π⁰-decay flux for each IEM, because the fit results for the combined tangent ranges are within these uncertainties when fitted to the positive and negative ranges individually. Above 1 GeV this is ∼4 × 10⁻⁸ ph cm⁻² s⁻¹ for the 15° × 15° region about the GC across all IEMs. This is comparable to the fitted flux from annulus 1 π⁰-decay or the TS < 25 point sources over the same region.

Any analysis employing the Galactocentric annulus decomposition for the gas column densities is subject to the loss of kinematic resolution for sight lines within l ∼ ±12° of the GC/anti-GC. Appendix B of Ackermann et al. (2012a) details the transformation of H i and CO gas-survey data into the column density distributions over Galactocentric annuli used in this analysis, and employed by many others. The assumptions made in the transformation for the site lines over the 15° × 15° region about the GC have an impact on the interstellar emission and point sources in the maximum-likelihood fitting and consequently the spatial distribution of residuals. Approximations made interpolating the gas column density across the l ± 10° range can result in an incorrect gas density distribution along the line of sight. Spurious point sources in the analysis and structure in residuals can result from this because a higher/lower CR intensity compared to where the gas should be placed is used in creating the interstellar emission templates. The scaling procedure for the IEM then adjusts the individual annuli potentially producing low-level artifacts due to a combination of the effects described above.

To obtain an estimate of the uncertainties associated with misplacement of the gas new maps of the column density per annuli are created. 10% of the H i gas column density is randomly displaced over the annuli and recombined with the π⁰-decay emissivity⁹³ in each annulus to create modified intensity maps for this process, which are summed to produce new fore-/background intensity maps. The 68% fractional change per pixel from 100 such realizations for each IEM is compared with the fore-/background resulting from the scaling procedure (Section 3.1). Depending on the IEM and energy range, variations from 1% to 15% in the intensity per pixel for the fore-/background from the structured interstellar emission across the 15° × 15° region are obtained, with the largest for OBstars index-scaled and smallest for the Pulsar intensity-scaled IEM, respectively. Because of the somewhat arbitrary choice of the precise fraction of H i column density⁹⁴ that is redistributed over the annuli these variations are illustrative rather than providing a true "systematic uncertainty" associated with the gas misplacement. Note that the uncertainty is maximized toward the GC because it is furthest away from the gas column density interpolation base points at l ∼ ±12°.