DEVELOPMENT OF THE MODEL OF GALACTIC INTERSTELLAR EMISSION FOR STANDARD POINT-SOURCE ANALYSIS OF FERMI LARGE AREA TELESCOPE DATA

Approxiamtely 99% of the ISM mass is gas and about 70% of this mass is hydrogen. The hydrogen gas exists in the form of neutral atoms in cold and warm phases, in the form of neutral H₂ molecules, and in an ionized state (diffuse ${{\rm{H}}}^{+}$ and H ii regions; Heiles & Troland 2003). Helium and heavier elements are considered to be uniformly mixed with the hydrogen. The warm H i medium and, to some extent, the cold H i medium are traced by 21 cm line radiation. Most of the cold molecular mass is traced by ¹²CO line emission. At the interface between the atomic and molecular phases, a mixture of dense H i and diffuse H₂ escapes the H i and CO radio surveys because this intermediate medium is optically thick to H i photons and CO molecules are absent or weakly excited. This DNM can be indirectly traced by its dust and CR content (Grenier et al. 2005b; Ade et al. 2014b).

2.1.1. Atomic Hydrogen

We have derived the atomic column density (${N}_{{\rm{H}}{\rm{I}}}$) maps from the 21 cm line radiation temperature under the assumption of a uniform spin (excitation) temperature (T_S). We used the 21 cm all-sky Leiden–Argentine–Bonn (LAB) composite survey of Galactic H i (Kalberla et al. 2005) to determine the all-sky distribution of ${N}_{{\rm{H}}{\rm{I}}}$. The LAB survey merges the Leiden/Dwingeloo Survey (LDS) of the sky north of $\delta =-30^\circ$ with the Instituto Argentino de Radioastronomia Survey (IAR) of the sky south of δ =  −25°. The spatial resolution after regridding is 35' in the case of the IAR survey, and 40' in the case of the LDS. We have derived ${N}_{{\rm{H}}{\rm{I}}}$ from the observed brightness temperature T_B with Equation (1):

$$\begin{eqnarray}&&{N}_{{\rm{H}}{\rm{I}}}=-1.82\times {10}^{18}\;{T}_{S}{\displaystyle \int }_{v}\mathrm{ln}\left[1-\displaystyle \frac{{T}_{B}(v)}{{T}_{S}-{T}_{0}}\right]\;{dv}\;\;[{\mathrm{cm}}^{-2}],\end{eqnarray} \tag{ 1 }$$

where the integral is taken over the velocity range of interest (Ackermann et al. 2012d) and T₀ = 2.66 K represents the background brightness temperature in this frequency range.

Since the cold H i clumps are embedded in the more diffuse warm gas, emission from both cold and warm atomic media can be detected in the same line of sight. Studies of H i absorption against background radio sources have shown that T_S is not uniform in the multi-phase ISM (Heiles & Troland 2003; Kanekar et al. 2011). Heiles & Troland (2003) found that most of the H i in the cold neutral medium has a T_S of less than 100 K and that the warm neutral medium gas lies roughly equally in the thermally unstable region (500−5000 K) and in the stable phase above 5000 K. In the absence of T_S information outside the small samples of background radio sources, we cannot predict the variations of T_S across the H i LAB survey. Instead, we selected the single uniform temperature that provided the best fit to the LAT data (see Section 3) in the anticenter region, at $90^\circ \leqslant l\leqslant 270^\circ$ and $| b| \lt 70^\circ$. Changes in T_S indeed modify the spatial distribution of ${N}_{{\rm{H}}{\rm{I}}}$ as seen from Earth, and these changes can be probed by the γ-ray emission produced by CRs in the atomic gas. The anticenter region was chosen because the uncertain IC emission is dimmer in the outer Galaxy; this is due to the fact that these directions are free of large and bright extended sources unrelated to the gas, such as Loop I and the Fermi bubbles, and because the Galactic warp creates a characteristic signature in the gas maps beyond the solar circle. In that region, the LAT observations of the diffuse emission are accurately reproduced from the gas, dust, and IC distributions. As we describe below, we have traced additional hydrogen in the DNM from dust column densities. For a uniform dust-to-gas ratio and uniform grain radiative emissivity, the DNM map partially corrects for local decreases in T_S. We have investigated seven values of T_S from 90 K to 400 K. For each temperature, we have derived H i and DNM column density maps (see Section 2.1.4). Figure 2 shows the log-likelihood ratio obtained by fitting the seven models at all γ-ray energies. The model obtained with T_S = 140 K gives the best fit to the LAT photon maps. This average temperature is smaller than the 250–400 K values found with H i absorption and emission spectra at 1'–2' resolution beyond the solar circle (Dickey et al. 2009), but it falls well within the temperature distribution found in the inner Galaxy, which peaks between 100 K and 200 K (Dickey et al. 2003), and so we have adopted this temperature of 140 K to derive H i column densities in the whole Galaxy in order to model the LAT data. We postpone dedicated studies of H i complexes in the outer Galaxy to understand why we find a larger fraction of cold H i than the radio studies, even though our model includes additional DNM gas.

Zoom In Zoom Out Reset image size

Because our model assumes a uniform CR density in each template, it is crucial to partition the Galaxy into Galactocentric annuli to account for radial variations in the CR density. The radial velocities measured from the Doppler shifts of the 21 cm line radiation can provide the kinematic Galactocentric distances of H i clouds. We assume that the gas moves in circular orbits around the Galactic Center (GC) and use the rotation curve given by Clemens (1985) with a Galactocentric distance and speed of the Sun of 8.5 kpc and 220 km s⁻¹, respectively. Measuring the Galactic rotation curve from 21 cm and 2.6 mm surveys was a well-understood procedure by the 1980s; for the inner Galaxy, the process primarily involved the measurement of the terminal velocity as a function of longitude. The Milky Way is a barred spiral, and a rotation curve does not strictly apply within Galactocentric radius range of the bar. However, the relatively coarse binning in radius that our model fitting requires mitigates the effect of noncircular motions in a barred potential. Regions located within 10° of the GC and anticenter have poor kinematic distance resolution. We have linearly interpolated the column densities in Galactic longitude from the column density integrated within 5° of longitude from their boundaries and scaled each line of sight to the total ${N}_{{\rm{H}}{\rm{I}}}$. The innermost H i annulus, where this interpolation is not possible, is assumed to contain 60% more gas than its neighboring annulus. We have excised the nearby galaxies LMC, SMC, M33, and M31, which were within the velocity range of the LAB survey. The partitioning of atomic hydrogen and CO into annuli is detailed in Appendix B of Ackermann et al. (2012d). We have generated nine annuli (see Figure 3) with the limits given in Table 1. Annulus number 7, referred to as "local," spans the solar circle. For each annulus, since the scale height of the CR distribution is several times greater than that of the gas, we assumed a uniform density in the axes perpendicular to the Galaxy.

Zoom In Zoom Out Reset image size

Table 1.  Limits of the ${N}_{{\rm{H}}{\rm{I}}}$ Galactocentric Annuli in Kiloparsecs

Annulus	rmin (kpc)	rmax (kpc)
1	0.0	1.5
2	1.5	4.5
3	4.5	5.5
4	5.5	6.5
5	6.5	7.0
6	7.0	8.0
7	8.0	10.0
8	10.0	16.5
9	16.5	50.0

Download table as:  ASCIITypeset image

2.1.2. Molecular Hydrogen

2.1.3. Ionized Hydrogen

2.1.4. Dark Neutral Medium

The ¹²CO J($1\to 0$) line emission is not a perfect tracer of cold H₂. In addition to metallicity variations, the CO molecule is strongly affected by UV photo-dissociation in the outer regions of molecular clouds where H₂ can exist without CO or where CO is only weakly excited (Wolfire et al. 2010; Glover & Low 2011). ${X}_{\mathrm{CO}}$ also depends on the dynamical characteristics of the molecular cloud. Shetty et al. (2011) went beyond the simple "mist" model and calculated the radiative transfer of the CO line in turbulent clouds in order to investigate the dependence of ${X}_{\mathrm{CO}}$ on the physical properties of gases. They predicted that the typical ranges of mean column density, temperature, and velocity dispersion found in molecular clouds lead to a variation in ${X}_{\mathrm{CO}}$ by about a factor of two. Moreover, as we already mentioned, the column density of H i derived from the 21 cm brightness temperature under the hypothesis of a uniform T_S = 140 K is likely to be biased in lines of sight for which T_S is not uniformly 140 K. Those limitations lead to large underestimates of the quantities of gas at the H i–H₂ interface in our Galaxy (Grenier et al. 2005b; Ade et al. 2011; Paradis et al. 2012) as well as in the Magellanic Clouds (Bernard et al. 2008; Dobashi et al. 2008). This transitional region is referred to as the DNM and for any particular region it comprises unknown fractions of cold dense H i and CO-free or CO-quiet H₂.

Neutral gas and large dust grains coexist; observations have shown that the gas-to-dust ratio leads to a mass ratio of ${M}_{\mathrm{gas}}/{M}_{\mathrm{dust}}\sim 100$. The dust column density can therefore provide an alternative template for H i and H₂ in the Milky Way. We used the dust reddening map of Schlegel et al. (1998), which is based on the 6' resolution IRAS/ISSA emission map at 100 μm and on a correction for the dust temperature inferred from the emission ratio measured between 100 and 240 μm at 07 resolution with COBE and DIRBE. The resulting map, once scaled to reddening to match the data from elliptical galaxies, has often been used as an estimator of Galactic extinction. It has been recently superseded by dust optical depth derivations at higher angular resolution and with a broader spectral coverage thanks to the Planck data (Abergel et al. 2013), but those were not available for the development of the present GIEM for Fermi-LAT.

Away from photo-dissociation and hot star-forming regions where dust properties can vary dramatically, most of the dust column density is well correlated with ${N}_{{\rm{H}}{\rm{I}}}$ and W(CO). Regions of dust not correleted with ${N}_{{\rm{H}}{\rm{I}}}$ and W(CO), which spatially correlate with diffuse γ-ray excesses over H i and CO expectations, correspond to DNM-rich clouds (Grenier et al. 2005b). We have derived residual maps by subtracting from the dust reddening map those parts linearly correlated with the ${N}_{{\rm{H}}{\rm{I}}}$ and W(CO) annuli. We first filtered out IR point sources present in the dust map and applied inpainting methods based on wavelet decomposition (Elad et al. 2005) to reconstruct the signal at the point-source locations. We then proceeded in steps in order to minimize the impact of the model uncertainties in the inner Galaxy onto the closer annuli: we first fit the reddening map to the local gas annuli for $| b| \gt 10^\circ ;$ then, we fixed them and fitted outer annuli 8 and 9 in the anticenter region ($90^\circ \lt l\lt 270^\circ$); then, we fixed the previous annuli and fit the inner Galaxy. Subtracting the correlated parts from the total dust has revealed coherent structures across the sky in both the positive and negative residuals (see Figure 4). We note that when we fit the H i maps to derive the T_S that gives the best fit to the Fermi data (Figure 2), we used corresponding T_S values to extract the dust residual maps.

Zoom In Zoom Out Reset image size

The residuals shown in Figure 4 are statistically significant above 0.04 mag. The positive residuals reveal additional gas aside from that traced by ${N}_{{\rm{H}}{\rm{I}}}$ and W(CO). We associate this excess map with the DNM distribution. For a standard ${N}_{{\rm{H}}}/E(B-V)$ ratio (Bohlin et al. 1978; Casandjian 2015b), they translate to column densities in excess of 10²¹ cm⁻² in nearby clouds and in more distant regions of the Galactic disk. There is no velocity information associated with the detection of IR dust emission, and so we could not partition the DNM into annuli. This means that all DNM clouds effectively have the same CR flux in the γ-ray model fitting. Dense H i and diffuse H₂ dominate the DNM column densities, but the residual map also potentially incorporates ionized hydrogen mixed with dust. The ionized mass is, however, small compared to the DNM mass, which is known in nearby clouds to be comparable to or exceed the mass locked in the CO-bright H₂ cores (Grenier et al. 2005b; Ade et al. 2011, 2014b). From the DNM map, we have excised the Magellanic Clouds, M31, as well as those regions with a high density of IRAS sources, including the inner Galaxy for absolute longitudes less than 30° and latitudes less than 2° and part of the Cygnus region. We have also inpainted the excised regions using the same method mentioned above.

The negative residual map in Figure 4 exhibits strong deficits close to the Galactic plane where the H i and CO expectations exceed the data. They are likely related to regions in which an average T_S of 140 K is too low, and thus ${N}_{{\rm{H}}{\rm{I}}}$ is overestimated. The inclusion of these potential corrections to ${N}_{{\rm{H}}{\rm{I}}}$ in the γ-ray model does improve the fits to the interstellar γ-ray emission. In this paper, we call this map the "${N}_{{\rm{H}}{\rm{I}}}$ correction map." We also observe strong negative residuals in regions close to bright OB associations where the dust temperature corrections are too coarse because of the low spatial resolution of DIRBE. Such residuals in Orion have led to the detection of spurious γ-ray point sources in the 2FGL source catalog (flagged as "c" sources). We have zeroed those residuals in the Orion molecular cloud. We also observe small, diffuse residuals that extend to about 20° in latitude. They correspond to regions of warm dust with temperatures of 18–19 K and soft emission spectra (β < 1.5) in the recent spectral studies joining millimeter and IR data from Planck and IRAS (Abergel et al. 2013). They often surround regions of strong PAH emission or free–free emission, thereby suggesting dust exposed to a different ISRF than in the H i and CO clouds and for which the temperature corrections applied by Schlegel et al. (1998) were not precise enough. These residuals are thus more likely to reflect the emissivity changes of the dust grains rather than corrections to the total gas column densities.

In the local ISM, Grenier et al. (2005b) have shown that the use of H i, CO, and a dust residual map improves the fits to the interstellar γ-ray emission compared to using the sole dust map for the total gas. While the hydrogen is likely to have the same γ-ray emissivity per atom in the atomic, DNM, and molecular components of the ISM due to the good penetration of CRs at the GeV-TeV energies relevant for the LAT (Skilling & Strong 1976), dust opacity (optical depth per gas nucleon) is known to increase from the diffuse H i to dense H i and H₂ gas (Stepnik et al. 2003; Abergel et al. 2013, 2014; Ade et al. 2014b). Restricting the use of dust to tracing the DNM alleviates the impact of dust opacity gradients because the DNM spans a moderate range of gas volume densities and is rather diffuse in space. Using a unique dust column density map at high latitude instead of H i, CO, and dust-inferred DNM maps would also lead to discontinuity in the model close to the Galactic plane where H i and CO annuli must be used to follow CR variations with distance from the GC.

While the different gas column density maps offer spatial templates for γ-ray photons originating mainly from π⁰-decay and bremsstrahlung emission, there is no direct observational template for the IC emission. Instead, it must be calculated. We have used the prediction from the GALPROP code with GALDEF identification ${}^{S}{Y}^{Z}{6}^{R}{30}^{T}{150}^{C}2$. This model features a radial distribution of CR sources proportional to the distribution of pulsars in the Galaxy given by Yusifov & Küçük (2004), a Galactic halo size and radial boundary equal to 6 kpc and 30 kpc, respectively, and the representative diffusive re-acceleration model described in Ackermann et al. (2012d). In that work, the GALPROP code was run to obtain models for the primary and secondary CR electron and positron intensities and spectra throughout the Galaxy. GALPROP then folds the distributions with ISRF (Porter et al. 2008) to obtain the IC emissivity, which was integrated along the line of sight for each direction and energy range to obtain IC intensity sky maps. We use ${I}_{{\mathrm{IC}}_{p}}$ to denote the intensity of the predicted IC emission.

The sky at GeV energies also features large angular scale structures such as Loop I and the Fermi bubbles. In the radio loops, the γ rays are likely produced by a population of CR electrons trapped in the old supernova remnants (Grenier et al. 2005a; Casandjian & Grenier 2009; Ackermann et al. 2012d). The Fermi bubbles correspond to two lobes of hard emission extending north and south from the direction of the GC (Su et al. 2010; Ackermann et al. 2014). There is no accurate independent template for the γ-ray emission of those large structures at other wavelengths since the synchrotron maps available in the radio fold the CR electron distribution with the complex structure of the magnetic field, and X-ray maps trace the hot gas but not the high-energy particle content. In the first iteration of our fitting, we assumed a proportionality between the γ-ray intensity and the bright radio continuum emission of the north Polar Spur (NPS), which dominates the 408 MHz radio continuum intensity survey of Haslam et al. (1982) at large positive latitudes. To do so, we have subtracted the point sources from the radio map and selected a region around Loop I (first panel of Figure 5). In the absence of external maps of the Fermi bubbles, we assumed a uniform intensity template in the first iteration of our fitting (see Section 4 and Figure 5) and for the final model we extracted intensity maps from the LAT data (see Section 7).

Zoom In Zoom Out Reset image size

We have modeled each of the 2179 point sources derived with a preliminary iteration of the GIEM for the 3FGL catalog (Ballet & Burnett 2013; Acero et al. 2015). We modeled each source with the Science Tool gtsrcmaps,⁶⁹ which takes into account the exposure, the angular resolution, and the source spectrum at each source position and in each energy bin.

Ackermann et al. (2013) analyzed the source population to estimate the contribution of unresolved sources to the diffuse Galactic emission above 10 GeV, and found the contribution to be about 2.5% of the interstellar emission for a reference model with a local source density of 3 kpc${}^{-3}$. For a tenfold increase in the local source density in a "maximum density" model, their contribution rises to  8%. In Acero et al. (2015), it was estimated that their contribution at 1 GeV amounts to 3% in the inner Galaxy. In the present work, the γ rays produced by unresolved sources are likely accounted for by the other templates, in particular, the IC and inner H i templates.

We have built intensity maps for the following 21 extended sources: Centaurus A, Cygnus Loop, Gamma Cygni, HESS J1614−518, HESS J1616−508, HESS J1632−478, HESS J1825−137, HESS J1837−069, IC 443, LMC, MSH 15-52, Puppis A, RX J1713.7−3946, S147, SMC, Vela Junior, Vela X, W28, W30, W44, and W51C. These small extended sources are associated with specific supernova remnants, pulsar wind nebulae, and spatially resolved galaxies. We have modeled each source with a simple disk shape, a ring, a two-dimensional (2D) Gaussian function, or a map derived from other wavelengths, as described in Nolan et al. (2012) and Ballet & Burnett (2013).

The time-integrated γ-ray emissions from the Sun and the moon effectively add a diffuse glow across the sky at low ecliptic latitudes. Their intensities and spatial distributions have been calculated by Abdo et al. (2011, 2012) and used here.

The isotropic emission encompasses the isotropic diffuse γ-ray background (Abdo et al. 2010b; Ackermann et al. 2015) originating from unresolved sources like blazars, star-forming and radio galaxies, as well as contamination from the very small fraction of CRs interacting in the LAT that are misclassified as γ rays and from Earth-limb photons that enter the LAT from the back but are reconstructed to have come from within the field of view. Because of this contamination, the isotropic emission depends on the event class and the conversion type in the LAT (Ackermann et al. 2012b). We have modeled this emission with an isotropic intensity template, with an intensity spectrum to be obtained from the fit to the γ-ray in each energy band.

Due to their proximity, CR protons and electrons interacting with Earth's atmosphere make it by far the brightest γ-ray source in the sky, with an intensity ∼1000 times larger than that of the Galactic plane (Abdo et al. 2009b). The Fermi standard observational strategy is such that Earth is not directly in the field of view of the LAT. However, a large number of limb photons entering the LAT from the side are still detected. We have removed the great majority of those photons with a simple cut in the zenith angle at 100°, but a residual contamination coming from the tails of the PSF is still observed at energies below about 200 MeV. Accurately simulating this component of the Earth limb photons in the far tails of the PSF is challenging, and so we have chosen instead to construct a simple template based on the subtraction of the map derived with a zenith angle cut at 80° from one restricted to angles above 100° for energies between 40 MeV and 80 MeV. We have deconvolved the resulting map to account for PSF broadening. We used the assumption that the spatial distribution of residual limb photons, which over the long time interval analyzed here is largely determined by the inclination and altitude of the orbit, is independent of the energy. The all-sky distribution is displayed in Figure 6 at the γ-ray energy of 100 MeV. The spectrum is soft and well fit by the power law $4.13\times {E}^{-4.25}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}\;{\mathrm{MeV}}^{-1}\;{\mathrm{sr}}^{-1}$, where E denotes the γ-ray energy in MeV. Above 200 MeV, its contribution becomes negligible for our analysis.

Zoom In Zoom Out Reset image size

In order to measure the hydrogen emissivity spectra in the various gas phases near the Solar circle and in the outer Galaxy, we have performed a series of successive fits. We took advantage of the broad extent in latitude of the local gas to reduce the influence of the much brighter Galactic ridge.

In addition to the NPS radio template, we have built simple intensity patches to account for the sources of EEE. The patches encompass regions with an excess of photons of at least 20% compared to the best-fit preliminary model (with all parameters in Equation (2) left free). This cut is well above the average level of positive or negative residuals found in the rest of the sky. Figure 5 shows the location and extent of the seven patches. The first four uniform patches fill the region toward the northern part of Loop I, a disk-shaped patch covers the region around the cocoon of hard γ rays observed by Ackermann et al. (2011a) toward Cygnus, and the last two patches cover the Fermi bubbles. Each patch has a spatially uniform intensity spectrum included, with the NPS template, in the fits to the data as additional free components in Equation (2). We note that the GIEM does not include the NPS map or the uniform templates.

The results of the gas emissivity spectra for Galactocentric radii greater than 7 kpc have been obtained by applying the following longitude and latitude cuts and by performing the series of successive fits in each of the 14 energy bins in the following order:

• 1.  
  H i emissivity in the 8–10 kpc annulus about the Solar circle, ${q}_{{\rm{H}}{{\rm{I}}}_{7}}(E)$: all longitudes and $10^\circ \lt | b| \lt 70^\circ$;
• 2.  
  CO emissivities about the Solar circle, ${q}_{{{\rm{CO}}}_{6+7}}(E)$: all longitudes, $4^\circ \lt | b| \lt 70^\circ$, where CO emission is significantly detected in the moment-masked filtered maps;
• 3.  
  gas emissivity in the DNM, ${q}_{\mathrm{DNM}}(E)$, H i emissivity in the 7–8 kpc annulus, ${q}_{{\rm{H}}{{\rm{I}}}_{6}}(E)$, and emissivity associated with the ${N}_{{\rm{H}}{\rm{I}}}$ corrections: all longitudes and $3^\circ \lt | b| \lt 70^\circ$;
• 4.  
  H i and CO emissivities in the outer Galaxy, ${q}_{{\rm{H}}{{\rm{I}}}_{8}}(E)$, ${q}_{{\rm{H}}{{\rm{I}}}_{9}}(E)$, and ${q}_{{{\rm{CO}}}_{8+9}}(E)$: all latitudes and $90^\circ \lt l\lt 270^\circ$.

We have checked that the measured emissivity spectra do not significantly depend on the precise shapes of the patches as long as they approximately follow the edges of the EEE and contain most of its emission. The emissivity spectra in the outer annuli are rather insensitive to the patches as the latter gather in the first and fourth Galactic quadrants. Absolute latitudes have several times been restricted to 70° because the isotropic emission dominates at higher latitudes, as described below. For those independent fits, we left all of the template normalization coefficients of Equation (2) free to vary in each of the 14 energy bins except for the Sun and Moon templates.

The inner Galactic region is particularly difficult to model. The gas column densities are most strongly affected by optical depth corrections and self-absorption, as well as by limited kinetic distance resolution at low longitudes. The determination of dust reddening for the DNM is less precise (Abergel et al. 2013) and we cannot trace CR density variations with distance in this component. Additionally, γ-ray point-like and extended sources are numerous and the IC_p morphology is uncertain. For one or several of these limitations in the interstellar modeling, or because of an over-density of CRs in certain regions, we observe EEE at low latitudes toward the inner Galaxy, with a maximum in the first Galactic quadrant near the base of the NPS (l ∼ 30°). To ensure that the EEE are not taken up at low energies by undue softening of the individual point sources because of the broad PSF, we have fixed the source intensities to those found in the first iteration of the 3FGL catalog, which was calculated with all of the components of the diffuse model and patches set free. The use of uniform intensity patches to account for the EEE would lead to a strong dependence of the inner-annuli emissivities on the patch shapes. Instead, we have used a two-step procedure in which the shape and spectrum of the EEE are iteratively determined from residual emission and incorporated in the γ-ray model in order to measure the gas emissivities. The ${N}_{{\rm{H}}{\rm{I}}}$ and W(CO) annuli in the inner Galaxy have limited morphological differences. With the added presence of the bright, unknown EEE, it became necessary to reduce the number of free components in the inner-Galaxy parts of the model and to proceed with a single hydrogen template for each annulus.

In the first step, we took advantage of the smaller extent and reduced intensity of the EEE in the fourth Galactic quadrant to measure the ${X}_{\mathrm{CO}}={N}_{{{\rm{H}}}_{2}}$/W(CO) conversion factors for the inner annuli 1–5, assuming that the γ-ray emissivity of the H₂ molecule is twice that of H i. Under this assumption the ${X}_{\mathrm{CO}}$ factor is half the ratio between the H i and CO emissivities: ${X}_{{{\rm{CO}}}_{i}}={q}_{{{\rm{CO}}}_{i}}/(2{q}_{{\rm{H}}{{\rm{I}}}_{i}})$. We have fitted Equation (2) to the LAT observations (without patches) in the fourth Galactic quadrant and GC region ($270^\circ \lt l\lt 365^\circ$). We have left all of the components free to vary, including the local annulus, which could be constrained by the latitude extent, $| b| \lt 20^\circ$, of this fit. This was done to optimize the determination of the gas emissivities in the inner annuli. To account for the EEE, we have selected the positive residuals, smoothed them with a 2D Gaussian symmetric kernel of 3° FWHM to weaken the correlation with the gas distributions, and re-injected them as an additional component with a free intensity for the next iteration. We have iterated three times up to the point where the distributions of the positive and negative residual intensities became comparable. We have selected the ${X}_{\mathrm{CO}}$ values obtained for the energy band centered on 2 GeV where high statistics and the narrow PSF reduce the cross-correlation between the H i and CO maps. We have obtained ${X}_{\mathrm{CO}}$ = 0.5 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹ for annuli 1 and 2, and X_CO = 1.5 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹ for annuli 3, 4, and 5. We did not detect a significant signal from the atomic hydrogen in the CMZ, but the gas of this region is largely molecular.

For the second step, we fit the whole Galactic disk, $0^\circ \leqslant l\lt 360^\circ$ and $| b| \lt 20^\circ$, with the combined ${N}_{{\rm{H}}}={N}_{{\rm{H}}{\rm{I}}}$ + 2${X}_{\mathrm{CO}}$W(CO) maps with free emissivities for the inner annuli 1 to 5, with the W(CO) CMZ map with a free emissivity, with the H i and CO annuli (not combined) beyond 7 kpc with free emissivities, and with all of the other model components of Equation (2) left free to vary, except for the sources. We used the same iterative procedure to account for the EEE. With this fit, we obtained the γ-ray emissivity spectra per hydrogen atom in the inner annuli.

Figure 7 shows the spectral energy distributions (SEDs) of the emissivities obtained for the nine Galactocentric annuli and for the CMZ region from the differential γ-ray emissivities per hydrogen atom $\tfrac{{dq}}{{dE}}=q/{\rm{\Delta }}E$, where ${\rm{\Delta }}E$ is the energy bin width. We did not correct the emissivities for the LAT energy dispersion. The emissivities plotted in Figure 7 overestimate the true values by approximately 10% at 100 MeV and 50% at 50 MeV; it is negligible above 200 MeV (Casandjian 2015b). As described above, the emissivity per hydrogen atom is derived in the H i phase for annuli 6–9 and in both the atomic and molecular hydrogen phases for annuli 1–5. The emissivity in the CMZ, measured from the CO map, was scaled to emissivity per hydrogen atom assuming ${X}_{\mathrm{CO}}$ = 0.5 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹, which is the same as what we found for annuli 1 and 2. The low ${X}_{\mathrm{CO}}$ value of 0.2–0.7 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹ inferred from infrared observations by Sodroski et al. (1995) in the GC region also supports this choice. We observe that the emissivities follow continuous energy distributions even though the γ-ray fits were performed independently in each energy bin. The error bars plotted in Figure 7 represent only the statistical uncertainties.

Zoom In Zoom Out Reset image size

For the local ISM, in Casandjian (2015b), we show that the systematic uncertainties associated with the measurement of the hydrogen γ-ray emissivity are dominated by the uncertainty in the LAT effective area over the whole energy range. It amounts to 10% below 100 MeV, with a linear decrease in $\mathrm{log}(E)$ to 5% in the range between 316 MeV and 10 GeV, and a linear increase in $\mathrm{log}(E)$ up to 15% at 1 TeV.⁷² Those values represent a lower limit for the hydrogen emissivities in the outer Galaxy where the non-uniformity of T_S cannot be neglected.

For the hydrogen emissivities in the inner annuli ($\lt 7$ kpc), the major source of systematic uncertainties is the model incompleteness. We note that when we restrict the fit to the fourth Galactic quadrant without using any residual template or patch for the EEE, the gas emissivities in the CMZ and in inner annuli 2 and 4 increase globally by up to 40% compared to those given in Figure 7. A second source of uncertainty is the amount of dense, cold H i in the inner spiral arms and the spin temperature correction. The H i mean opacity increases inward from the solar circle and peaks in the molecular ring, suggesting that the cold phase is more abundant and colder there than it is locally (Dickey et al. 2003). However, the median fraction of cold H i is about 20% in column density (Heiles & Troland 2003; Murray et al. 2015), so that the measurements of the γ-ray emissivities per H atom are constrained primarily by the gas content in the warm H i phase which is more reliably quantified. A change in T_S from 140 K to 400 K (Dickey et al. 2009) results in a ∼30% change in ${N}_{{\rm{H}}{\rm{I}}}$ for the inner annuli. A third source of uncertainty is the quality of the dust map at low latitudes, which is hampered by the poor angular resolution of the temperature corrections against a rising spatial density of warm star-forming regions. There exist localized differences by a factor up to 2 or 3 along the Galactic plane between the estimates recently inferred at 5' with Planck (Abergel et al. 2013) and the Schlegel et al. (1998) map available at the start of this work; the outer Galaxy reddening is also significantly underestimated in the Schlegel et al. (1998) map, while the inner Galaxy is too dusty. These differences do not so much affect the gas emissivity in the DNM (because of the presence of massive clouds off the plane to help the fits) as the ${N}_{{\rm{H}}{\rm{I}}}$ correction. We estimate a 30%–40% uncertainty which can partially propagate to the uncertainties of the emissivities in the inner annuli.

We fit the differential emissivity SEDs with a model of bremsstrahlung emission (Gould 1969) and hadronic decay. Most γ rays with energies between 100 MeV and 50 GeV originate from the decay of π⁰ produced in hadronic collisions when CR protons with energies above 0.5 GeV interact with protons from ISM nuclei (Stecker 1970; Aguilar-Benitez 1991). We fit the emissivity SED of each annulus between 200 MeV and 30 GeV using the γ-ray production cross-section of Kamae et al. (2006) and a CR proton flux spectrum parametrized as: $A{\beta }^{{P}_{1}}{R}^{{P}_{2}}$, where $\beta =v/c$ is the proton-to-light velocity ratio, R is the proton rigidity, A is a free normalization, and P₁ and P₂ are free spectral indices (Shikaze et al. 2007). This form tends to a power law with index P₂ at high rigidity (energy) and the P₁ index controls the spectral fall-off at low rigidities, which is indicated by the Voyager 1 data near the heliopause and by γ-ray measurements in the local ISM (Casandjian 2015b; Grenier et al. 2015). We used the results of Mori (2009) to scale the proton-proton cross-section to the nucleus–nucleus one, taking into account the abundance of heavier nuclei in the ISM and in the CRs (Casandjian 2015b). We have accounted for the bremsstrahlung radiation contribution using the following electron spectral form: $B{({E}_{\mathrm{kin}}/{E}_{0}+{(({E}_{\mathrm{kin}}-{E}_{4})/{E}_{0})}^{-0.5})}^{{P}_{3}}$, where E_kin is the kinetic energy of the electrons, ${E}_{0}=1\;{\rm{GeV}}$, B is a free normalization, P₃ is a free spectral index, and E₄ is a free scaling energy. We fit the bremsstrahlung emission together with the hadron decay component to the emissivity spectrum measured in the local annulus 7 (8–10 kpc). We found that the modeled emissivities systematically underestimate the measured ones by about 20% above 2 GeV. However, the measured emissivity spectrum compares well, within ∼10%, with previous measurements obtained with the same T_S in specific regions of the local ISM where there is minimal confusion along the lines of sight and where higher-resolution ISM data were used (see Figure 4 of Grenier et al. 2015). For every annuli, we corrected for this discrepancy by increasing the Kamae et al. production cross-section by 20% for γ-ray energies above 2 GeV, with a smooth transition to lower energies to avoid any discontinuity. Even though we applied this correction to the proton-proton cross-section, it could originate from an incorrect scaling factor from the proton-proton to nucleus–nucleus cross-section, or from the parametrization of the proton spectrum. A detailed interpretation of the emissivity SED, compatible with the one presented here, derived in the local annulus from a similar template fitting method is given in Casandjian (2015b).

From the CR spectra fit to the emissivity SED in the local annulus, we derived the bremsstrahlung contribution and the proton functional parameter P₁, which describes the low-rigidity turn-over of the proton spectrum. We assigned those same values to the other annuli. The bremsstrahlung radiation contributes only 14% at 200 MeV to the total emissivity spectrum in the local ISM (Casandjian 2015b). We interpreted the pion-decay emission above 1 GeV to study potential variations in flux and spectrum of the bulk of the CRs pervading the Galaxy. The corresponding CRs have energies per nucleon well above the low-energy turn-over (see Figure 1 of Grenier et al. 2015). At GeV energies and above, the improved LAT performance enables a better separation of the emissions originating from the different annuli. We fit the measured SEDs above 280 MeV assuming a uniform CR electron flux and P₁ across the Galaxy. As shown in Figure 7, we obtained a reasonable fit to the SEDs. The derived values of the proton index P₂ and normalization A are fairly insensitive to the approximations on P₁ and the bremsstrahlung contributions. We tested these approximations by also fitting the SEDs above 100 MeV with the electron normalization and index parameters allowed to be free and found no significant variation in P₂ and A. In Figure 7, in order to compare the gas emissivity spectrum in the DNM and in the local annulus, we used a gas-to-dust reddening ratio of 3.5 × 10²¹ cm⁻² mag⁻¹ (Grenier et al. 2005b). In the case of the ${N}_{{\rm{H}}{\rm{I}}}$ correction template, we used the same ${N}_{{\rm{H}}}/E(B-V)$ ratio and left all of the spectral parameters for electrons and protons free to obtain a better fit to the data.

Figure 8(a) shows the radial distributions of the γ-ray emissivity measured at 2 GeV across the Galaxy and Figure 8(b) shows the radial distribution of proton density integrated above 10 GV. We observe a marked increase in CR density around 3 kpc from the GC. The Fermi-LAT counts associated with the second gas annulus dominate the region within ±30° in longitude. The number of counts integrated above 2 GeV and associated with this annulus is twice that of 3FGL sources and four times that of IC in a region defined by the annulus contours. The EEE represents a few percent of the total in this region. This increase around 3 kpc might be associated with an enhanced CR production in the molecular ring. The steep increase relative to the next annulus is reminiscent of the marked increase in star formation rate indicated by massive stars (H ii regions, supernova remnants, pulsars; see Stahler & Palla 2005), as shown in Figure 8(d).

Zoom In Zoom Out Reset image size

The proton density profile predicted by the GALPROP model ${}^{S}{Y}^{Z}{6}^{R}{30}^{T}{150}^{C}2$ (Ackermann et al. 2012d) reproduces the trend of deviations from measurements by a factor of two at the maximum in the molecular ring region. For Galactocentric distances greater than 5 kpc, the predicted proton density gradient is steeper than the observed one. This discrepancy, referred to as "the CR gradient problem" has been known since the γ-ray surveys made by COS-B (Bloemen et al. 1986; Strong et al. 1988) and EGRET (Strong & Mattox 1996; Hunter et al. 1997). Bloemen et al. (1993) suggested that the radial distribution of CR sources may be flatter than inferred from pulsar and supernova remnant observations, or that the diffusion parameters derived from the local CR measurements are not the same throughout the Galaxy. A solution to this issue in terms of CR-driven Galactic winds and anisotropic diffusion has been proposed by Breitschwerdt et al. (2002). Uhlig et al. (2012) note that CR-driven winds could also suppress the star formation rate by a significant factor. Shibata et al. (2007) proposed a non-uniform diffusion coefficient that increases with Galactocentric radius and distance from the Galactic plane. Indeed, models with a large Galactic halo, and thus faster CR diffusion, are able to better reproduce the γ-ray emissivity in the outer Galaxy (Ackermann et al. 2012d). Such a position-dependent CR diffusion coefficient, linked to the ambient power in the magnetic turbulence induced by stellar and supernova activity, allows for a good reproduction of both Fermi-LAT and local CR observables (Gaggero et al. 2015 and references therein). On the other hand, more molecular gas in the outer Galaxy is still being discovered (Sun et al. 2015), implying that the star formation rate, and thus the density of CR sources, may be underestimated at large distances. An increase in ${X}_{\mathrm{CO}}$ with distance beyond the solar circle is expected because of the metallicity gradient (Pineda et al. 2013), but it cannot explain the large γ-ray emissivity values found in the outer Galaxy in correlation with the H i gas (Abdo et al. 2010a), contrary to what has been proposed by Strong et al. (2004b). The DNM gas at the H i–H₂ interface, however, is more abundant than CO-bright H₂ beyond the solar circle (Abergel et al. 2011; Pineda et al. 2013) and it should better correlate with H i spatially. It can offer an alternative or complementary solution to the CR gradient problem that can be tested using forthcoming radio line and dust extinction surveys.

We have plotted the proton spectral index P₂ versus the Galactocentric radius of the annulus in Figure 8(c). The index of 2.81 ± 0.08 found in the local annulus is fully consistent with the value of 2.820 ± 0.003 (stat.) ±0.005 (syst.) measured by the PAMELA experiment in the range 30 GV to 1.2 TV, which is beyond the influence of solar modulation (Adriani et al. 2011). This agreement indicates that the same CR population pervades the local spiral arm and the immediate solar neighborhood. Measurements of γ-ray emissivity in nearby clouds had indicated such a spectral uniformity out to a few hundred parsecs from the Sun (for a review see Section 4 of Grenier et al. 2015). The present result extends this uniformity to the Local Arm, which dominates the local annulus.

We observe a hardening of the proton spectra when moving from the outskirts of the Galaxy to the inner molecular ring. This spectral hardening cannot be due to contamination by the EEE at low latitudes since their emission contributes on average only 10% of the gas emission. Contamination by unresolved sources such as pulsars is also unlikely. They generally have much harder spectra than the interstellar emission (Abdo et al. 2013), but they contribute at best a few percent of the diffuse emission (see Section 2.4). One would need a huge increase in source density in the inner Galaxy, at variance with pulsar and supernova remnant observations. The proton spectral indices extracted in the CMZ and the first annulus (first two radial bins in Figure 8) are evaluated for a region extending ±10° from the GC which is dominated by the emission of point sources listed in the 3FGL catalog, together with IC; in this region, the confusion with gas emission is maximal. Moreover, as we discussed in Section 2.1.1, in this region, the gas column density is calculated in part using interpolations of adjacent regions and may therefore be inaccurate.

CR transport models such as GALPROP do not predict spectral variations across the Galaxy because they assume uniform diffusion properties and a uniform injection spectrum from CR sources. Gaggero et al. (2015) have also recently noted gradual CR hardening toward the inner regions under the cruder assumption that the γ radiation originating from the pion decays in the gas dominates all of the other emission components, so that the diffuse γ-ray spectrum above 5 GeV directly maps the CR spectrum. They propose to explain this hardening by varying the CR diffusion properties through the Galaxy, specifically by linearly decreasing the rigidity index δ of the diffusion coefficient toward the GC. This decrease allows for harder CR spectra at small Galactic radii and can also explain the emission deficit noted above a few GeV in the (inner) disk with uniform CR transport models. Their model also includes strong convection within 6.5 kpc from the GC. Energy-independent CR transport by Galactic winds keeps the spectrum near the hard distribution injected by the sources, and so an increasing wind toward the inner Galaxy would eventually dominate over energy-dependent diffusion. Convection in the form of Galactic winds is supported by X-ray (Everett et al. 2008) and radio (McClure-Griffiths et al. 2013) observations. Figure 8(c) compares the spectral variation across the Galaxy observed by Fermi-LAT data with the one predicted by Gaggero et al. (2015). We observe a reasonable agreement.

In order to build a 3D intensity cube in position and energy, ${I}_{\mathrm{EEE}}(l,b,E)$, of the EEE at angular scales relevant for point-source analyses, we derived their spectral distributions in each sky pixel and we used wavelet decomposition to retain structures on angular scales broader than 2°. We parametrized the spectral distributions by assuming IC interactions of a population of CR electrons with the CMB radiation. This parametrization is motivated by the presence of radio-emitting CR leptons in Loop I and the Fermi bubbles (Ade et al. 2013), but its goal is to provide a flexible parametric form to describe the spectrum without Poisson noise in each pixel of the sky, rather than to interpret the spectra. We discuss the specifics of this procedure below.

To construct the I_EEE, we rebinned the LAT count maps in each of the 14 energy bins to a $1\_\_AMP\_\_fdg;8$ grid. In each pixel of this grid, we fitted the electron power-law spectra so that the sum of the baseline model (see Section 7) and the supplementary electron IC emission reproduces the total photon count spectrum of the pixel. We need two independent power-law electron spectra in order to match the data over the whole γ-ray energy range from 50 MeV to 50 GeV. The power laws are, respectively, constrained by the LAT data below and above a γ-ray energy of 965 MeV. The resulting IC intensity maps are found to be in good agreement with those of the EEE presented in Section 7 between 50 MeV and 50 GeV. We have built sky maps of the power-law indices and normalizations for each electron population. To filter out point-like and small extended sources not present in our source list, as well as large Poisson fluctuations, we have transformed the spatial distributions of the electron spectral parameters into wavelets and filtered out scales smaller than 2°. Then, by applying the inverse transform, we derived separate IC intensity maps from the low-energy and high-energy filtered electron spectra. We applied a smooth spectral transition to merge the two IC distributions in each sky pixel. Finally, in order to reduce the amount of LAT data to be reintroduced into the GIEM, we restricted I_EEE to those regions where the EEE are bright. We further verified at energies above 50 GeV, in five energy bins spanning from 50 GeV to 600 GeV, that the sum of I_EEE and the baseline model agrees with the LAT observations. We note that we did not perform any deconvolution of the LAT counts map in order to avoid introducing structures with angular scales less than 2° into the model.

In the right-hand column of Figure 10, we show I_EEE obtained at energies close to the geometric averages of the energy intervals used to display the count maps of the left column. We observe good agreement between the large-scale structures in the LAT count maps and in the filtered component that we developed in order to account for these bright structures in the GIEM. In Figure 11(d), we show a close-up view of the residuals integrated between 1.7 GeV and 50 GeV in the direction of the GC when I_EEE is added to the baseline model. The residual map is flat apart from small-scale residuals toward the GC and the Fermi bubbles. They correspond to the small angular scales not retained in the wavelet decomposition for the construction of I_EEE.

We built the GIEM by summing the emission components originating from the gas annuli and the DNM (allowing for the dust-related correction in ${N}_{{\rm{H}}{\rm{I}}}$), from the Galactic IC emission, and from I_EEE (see Equation (3)). The gaseous components have been scaled according to the emissivity spectra $\tfrac{{{dq}}_{\mathrm{fit}}}{{dE}}$ obtained in Sections 5.1 and 5.2, and parametrized with the pion decay and bremsstrahlung decomposition described in Section 5.3 to provide continuous functions. They are represented by solid lines in Figure 7. The GALPROP IC distribution has been scaled with the normalization factors derived in Section 6 and interpolated in energy as necessary (see Figure 9). We stress that extended γ-ray sources with sizes larger than 2° are partially incorporated into the GIEM through the addition of I_EEE, and so they cannot be studied with this diffuse background model.

$$\begin{eqnarray}I(l,b,E) & = & \displaystyle \sum _{i=\mathrm{HI},{{\rm{H}}}_{2},\mathrm{DNM}}\displaystyle \frac{{{dq}}_{{\mathrm{fit}}_{i}}}{{dE}}(E){N}_{{H}_{i}}(l,b)\\ & & +{N}_{\mathrm{IC}}(E){I}_{{\mathrm{IC}}_{p}}(l,b,E)+{I}_{\mathrm{EEE}}(l,b,E)\end{eqnarray} \tag{ 3 }$$

The gas emissivities and IC normalization factors have been derived in the energy range 50 MeV to 50 GeV. Above this range, the low photon statistics do not allow for reliable component separation, especially toward the Galactic plane. To build the GIEM for the range 50–600 GeV, we extrapolated these factors as follows.

• 1.  
  The spectral form chosen for the CR protons pervading the gas is equivalent at large rigidity to a simple power law of index P₂. We checked the validity of this form for producing the right intensities beyond 50 GeV by comparing the LAT count maps recorded in 5 equal logarithmic bands between 50 GeV and 590 GeV to a model based on the extrapolation of the power-law proton spectra to several TV in rigidity. We observed an excess of high-energy γ rays in the LAT sky maps, including in the Galactic plane. To approximately account for this hardening, we globally scaled the gas emissivity spectra above 50 GeV by the coefficient $0.8+(5.9\times {10}^{-3}\times E)-(2.8\times {10}^{-6}\times {E}^{2})$ with the γ-ray energy E in GeV. This scaling would correspond to a gradual hardening for proton rigidities between approximately 200 GV and 2 TV. Near Earth, the compilation of PAMELA (Adriani et al. 2011), AMS-02 (Aguilar et al. 2015), CREAM (Yoon et al. 2011), and ATIC-2 (Panov et al. 2006) data indicates an up-turn in the proton spectrum starting around 500 GV. This break may signal a change in CR transport from diffusion on self-generated waves at low rigidities to diffusion on pre-existing magnetic turbulence at high rigidities (Aloisio et al. 2015). The present need to scale-up the local gas emissivity above 50 GeV to explain the LAT data at medium latitudes supports this CR hardening, but we defer a quantitative comparison to a later dedicated study. The marked need for a larger high-energy intensity at low latitudes may be due to a comparable break in CR spectra and transport properties further out in the Galactic disk, but it can also stem from a population of unresolved hard sources, such as pulsar wind nebulae that are abundantly detected as TeV sources.
• 2.  
  The Galactic IC intensity varies very smoothly across the sky. With increasing γ-ray energies, most of the IC emission occurs at low latitude toward the inner Galaxy. The lack of spatial variation at high latitude together with the low γ-ray statistics are such that the fit fails to reliably differentiate the IC from the isotropic emission above 50 GeV. Since the normalization of the GALPROP IC distribution is close to 1 at 50 GeV, we relied on GALPROP predictions, without rescaling, for the extrapolation to higher energies.

In Figure 12, we show a map of the interstellar γ-ray emission in the GIEM at about 1 GeV. In Figure 13, we compare the LAT count map integrated between 360 MeV and 50 GeV to that predicted by the GIEM, combined with point and extended sources from a preliminary version of the 3FGL catalog, the isotropic emission, and the emission from the Sun, the Moon, and from the residual Earth limb emission. The overall agreement between observations and the model is very good, partly because some of the detected excesses have been modeled and re-injected into the interstellar model. We still observe discrepancies along the Galactic plane at a level of less than 2σ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 14, we present the SEDs of various components of the GIEM, averaged over regions covering the Galactic disk ($| b| \lt 10^\circ$) and higher latitudes. We have decomposed the total SEDs into contributions originating from the interstellar gas (atomic, molecular, and DNM), from an axisymmetric Galactic ISRF for the IC radiation, and from the EEE. We have represented the correction to the H i contribution as an intensity by taking the negative of the photon intensity associated with the ${N}_{{\rm{H}}{\rm{I}}}$ correction map. Above 100 MeV, hadronic interactions with the atomic hydrogen dominate the interstellar γ-ray emission. The hardening of this emission above 50 GeV, obtained by scaling the gas emissivities as described above, is visible in both panels.

Zoom In Zoom Out Reset image size

The GIEM is available at the FSSC website as a FITS file named gll_iem_v05_rev1.fit. For this file, we have resampled all of the maps to a $0\_\_AMP\_\_fdg;125$ grid in CAR projection, which is the format required by the LAT Science Tools analysis software. The FITS file comprises 30 logarithmically spaced energies between 50 MeV and 600 GeV. It gives the photon specific intensity of the GIEM in photons sr⁻¹ s⁻¹ cm⁻² MeV⁻¹. This model, tuned to LAT data, is not corrected for the energy dispersion; therefore, it can be used directly with LAT data. Version 15 (V15) of the P7REP IRFs is the recommended version for Pass 7 reprocessed data and for this model. The difference with the P7REP_V10 set of IRFs used for the derivation of the gas emissivity spectra and IC normalization factors in this work mainly resides in an improved Monte Carlo PSF and in an updated fitting procedure to determine the parameters representing the LAT effective area. Those minor differences modify the exposure, but not the reconstructed LAT events. The minimum ratio in exposure (V15/V10) is 0.98 at 50 MeV and the maximum is 1.05 at 1 GeV. To correct the GIEM for the final (V15) IRFs, we have rescaled its intensity by the exposure ratio evaluated for each of the 30 energy planes of the GIEM. The model is then intended for use with the instrument response functions versions P7REP_SOURCE_V15, P7REP_CLEAN_V15, and P7REP_ULTRACLEAN_V15.

The GIEM aims to provide a representation of the interstellar emission that closely reproduces the LAT observations. It combines three methods: the robust template fitting method, which uses no assumption on CR transport and spectra, but is sensitive to cross-correlations between the diffuse components and to source confusion at the lowest energies; the prediction of a propagation model (GALPROP) for the Galactic IC emission, with a scalable intensity in energy, but a fixed spatial distribution at each energy; and an iterative detection and data-based modeling of the EEE for which we have no external information. The interplay between these methods is such that deriving the uncertainties of our model is very challenging.

At high latitudes and outside the region covered by I_EEE, the uncertainties are likely dominated by the determination of the gas emissivities. The thinness of the local H i is such that ${N}_{{\rm{H}}{\rm{I}}}$ is not very sensitive to variations in T_S, and so the uncertainties in the absolute determination of the LAT effective area dominate (Casandjian 2015b).

Toward the outer Galaxy and outside the region covered by I_EEE, uncertainties in the LAT effective area and in the uniformity of T_S both must be accounted for. We note that the column density in a line of sight can vary by up to a factor of 2 when assuming optically thin H i or a T_S of 95 K. The impact of a non-uniform T_S is partially reduced by the ${N}_{{\rm{H}}{\rm{I}}}$ correction applied with the dust-derived template. However, the dust optical depths are easily biased by temperature confusion at low latitudes (see our discussion in Section 5.3). The use of the "negative" ${N}_{{\rm{H}}{\rm{I}}}$ correction map in the model improves the global fits to the LAT data, but can artificially lower the interstellar emission in a specific region. We caution the reader concerning potentially spurious features in the direction of hot dust or of steep gradients in dust temperature.

In the other directions, the largest uncertainty in our model is its degree of incompleteness (large-scale sources off the plane, like Loop I and the Fermi bubbles, optically thick H i and CO, poor gas distribution with distance, and missing DNM mass in the inner Galaxy). We mitigated this incompleteness by including a component I_EEE extracted from the data at angular scales broader than 2°. However, in Figure 13 (bottom), we observe that our model is still not complete since small deviations remain visible above the statistical fluctuations. This is a consequence of the interplay between the different components in the fits, which converge to the least-bad solution rather than the best one (a great improvement in one zone can be reached at the expense of minor worsening in others). It is also a consequence of adding a filtered map issued only from the positive residuals. In Section 5.3, we also discussed the uncertainties on the gas emissivities in the inner annuli.

Although we cannot quantify the systematic uncertainties or the degree of incompleteness of our model, we have assessed several indicators for the quality of this work:

• 1.  
  coherent spectral distributions for the gas emissivities and IC normalization factors, in agreement with the gas emissivity spectra previously obtained in dedicated studies of less confused regions;
• 2.  
  coherent spatial structures of the EEE and I_EEE, strongly reminiscent of well-known features at other wavelengths;
• 3.  
  a coherent and continuous shape for the edges of the Fermi bubbles at low latitudes;
• 4.  
  the detection of the extended apparent path of the Sun and of the Moon across the sky in residual maps when they are not added to the model;
• 5.  
  a flat final residual map within ±2σ;
• 6.  
  the need for less than 5% corrections to the GIEM specific intensities when fitting the data in the large majority of the 840 regions of interest used in the generation of the 3FGL catalog (see Figure 25 of Acero et al. 2015).