<figure><img src="image_1.png" /></figure>

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

F. ACERO¹, M. ACKERMANN², M. AJELLO³, A. ALBERT⁴, L. BALDINI⁴,⁵, J. BALLET¹, G. BARBIELLINI⁶,⁷, D. BASTIERI⁸,⁹,
R. BELLAZZINI¹⁰, E. BISSALDI¹¹, E. D. BLOOM⁴, R. BONINO¹²,¹³, E. BOTTACINI⁴, T. J. BRANDT¹⁴, J. BREGEON¹⁵, P. BRUEL¹⁶,
R. BUEHLER², S. BUSON¹⁴,¹⁷,¹⁸, G. A. CALIANDRO⁴,¹⁹, R. A. CAMERON⁴, M. CARAGIULO¹¹, P. A. CARAVEO²⁰, J. M. CASANDJIAN¹,
E. CAVAZZUTI²¹, C. Cecchi²²,²³, E. Charles⁴, A. CHEKHTMAN²⁴, J. CHIANG⁴, G. CHIARO⁹, S. CIPRINI²¹,²²,²⁵, R. CLAUS⁴,
J. COHEN-TANUGI¹⁵, J. CONRAD²⁶,²⁷,⁶², A. CUOCO¹²,¹³, S. CUTINI²¹,²²,²⁵, F. D'AMMANDO²⁸,²⁹, A. DE ANGELIS³⁰, F. DE PALMA¹¹,³¹,
R. DESIANTE¹²,³², S. W. DIGEL⁴, L. Di VENERE³³, P. S. DRELL⁴, C. FAVUZZI¹¹,³³, S. J. FEGAN¹⁶, E. C. FERRARA¹⁴, W. B. FOCKE⁴,
A. FRANCKOWIAK⁴, S. FUNK³⁴, P. FUSCO¹¹,³³, F. GARGANO¹¹, D. GASPARRINI²¹,²²,²⁵, N. GIGLIETTO¹¹,³³, F. GIORDANO¹¹,³³,
M. GIROLETTI²⁸, T. GLANZMAN⁴, G. GODFREY⁴, I. A. GRENIER¹, S. GUIRIEC¹⁴,⁶³, D. HADASCH³⁵, A. K. HARDING¹⁴, K. HAYASHI³⁶,
E. HAYS¹⁴, J. W. HEWITT³⁷, A. B. HILL⁴,³⁸, D. HORAN¹⁶, X. HOU³⁹,⁴⁰, T. JOGLER⁴, G. JÓHANNESSON⁴¹, T. KAMAE⁴², M. KUSS¹⁰,
D. LANDRIU¹, S. LARSSON²⁷,⁴³, L. LATRONICO¹², J. LI⁴⁴, L. Li²⁷,⁴³, F. LONGO⁶,⁷, F. LOPARCO¹¹,³³, M. N. LOVELLETTE⁴⁵,
P. LUBRANO²²,²³, S. MALDERA¹², D. MALYSHEV³⁴, A. MANFREDA¹⁰, P. MARTIN⁴⁶, M. MAYER², M. N. MAZZIOTTA¹¹,
J. E. McENERY¹⁴,⁴⁷, P. F. MICHELSON⁴, N. MIRABAL¹⁴,⁶³, T. MIZUNO⁴⁸, M. E. MONZANI⁴, A. MORSELLI⁴⁹, E. NUSS¹⁵, T. OHSUGI⁴⁸,
N. OMODEI⁴, M. ORIENTI²⁸, E. ORLANDO⁴, J. F. ORMES⁵⁰, D. PANEQUE⁴,⁵¹, M. PESCE-ROLLINS⁴,¹⁰, F. PIRON¹⁵, G. PIVATO¹⁰,
S. RAINÒ¹¹,³³, R. RANDO⁸,⁹, M. RAZZANO¹⁰,⁶⁴, S. RAZZAQUE⁵², A. REIMER⁴,³⁵, O. REIMER⁴,³⁵, Q. REMY¹, N. RENAULT¹,
M. SÁNCHEZ-CONDE²⁷,²⁶, M. SCHAAL⁵³, A. SCHULZ², C. SGRÒ¹⁰, E. J. SISKIND⁵⁴, F. SPADA¹⁰, G. SPANDRE¹⁰, P. SPINELLI¹¹,³³,
A. W. STRONG⁵⁵, D. J. SUSON⁵⁶, H. TAJIMA⁴,⁵⁷, H. TAKAHASHI⁵⁸, J. B. THAYER⁴, D. J. THOMPSON¹⁴, L. TIBALDO⁴, M. TINIVELLA¹⁰,
D. F. TORRES⁴⁴,⁵⁹, G. TOSTI²²,²³, E. TROJA¹⁴,⁷⁰, G. VIANELLO⁴, M. WERNER³⁵, K. S. WOOD⁴⁵, M. WOOD⁴,
G. ZAHARIJAS⁶⁰,⁶¹, AND S. ZIMMER²⁶,²⁷

¹ Laboratoire AIM, CEA-IRFU/CNRS/Université Paris Diderot, Service d'Astrophysique, CEA Saclay, F-91191 Gif sur Yvette, France;
isabelle.grenier@cea.fr, casandjian@cea.fr

² Deutsches Elektronen Synchrotron DESY, D-15738 Zeuthen, Germany

³ Department of Physics and Astronomy, Clemson University, Kinard Lab of Physics, Clemson, SC 29634-0978, USA

⁴ W. W. Hansen Experimental Physics Laboratory, Kavli Institute for Particle Astrophysics and Cosmology, Department of Physics and SLAC National Accelerator Laboratory, Stanford University, Stanford, CA 94305, USA

⁵ Università di Pisa and Istituto Nazionale di Fisica Nucleare, Sezione di Pisa I-56127 Pisa, Italy

⁶ Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, I-34127 Trieste, Italy

⁷ Dipartimento di Fisica, Università di Trieste, I-34127 Trieste, Italy

⁸ Istituto Nazionale di Fisica Nucleare, Sezione di Padova, I-35131 Padova, Italy

⁹ Dipartimento di Fisica e Astronomia “G. Galilei,” Università di Padova, I-35131 Padova, Italy

¹⁰ Istituto Nazionale di Fisica Nucleare, Sezione di Pisa, I-56127 Pisa, Italy

¹¹ Istituto Nazionale di Fisica Nucleare, Sezione di Bari, I-70126 Bari, Italy

¹² Istituto Nazionale di Fisica Nucleare, Sezione di Torino, I-10125 Torino, Italy

¹³ Dipartimento di Fisica Generale “Amadeo Avogadro,” Università degli Studi di Torino, I-10125 Torino, Italy

¹⁴ NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

¹⁵ Laboratoire Univers et Particules de Montpellier, Université Montpellier, CNRS/IN2P3, Montpellier, France

¹⁶ Laboratoire Leprince-Ringuet, École polytechnique, CNRS/IN2P3, Palaiseau, France

¹⁷ Department of Physics and Center for Space Sciences and Technology, University of Maryland Baltimore County, Baltimore, MD 21250, USA

¹⁸ Center for Research and Exploration in Space Science and Technology (CRESST) and NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
¹⁹ Consorzio Interuniversitario per la Fisica Spaziale (CIFS), I-10133 Torino, Italy

²⁰ INAF-Istituto di Astrofisica Spaziale e Fisica Cosmica, I-20133 Milano, Italy

²¹ Agenzia Spaziale Italiana (ASI) Science Data Center, I-00133 Roma, Italy

²² Istituto Nazionale di Fisica Nucleare, Sezione di Perugia, I-06123 Perugia, Italy

²³ Dipartimento di Fisica, Università degli Studi di Perugia, I-06123 Perugia, Italy

²⁴ College of Science, George Mason University, Fairfax, VA 22030, resident at Naval Research Laboratory, Washington, DC 20375, USA

²⁵ INAF Osservatorio Astronomico di Roma, I-00040 Monte Porzio Catone (Roma), Italy

²⁶ Department of Physics, Stockholm University, AlbaNova, SE-106 91 Stockholm, Sweden

²⁷ The Oskar Klein Centre for Cosmoparticle Physics, AlbaNova, SE-106 91 Stockholm, Sweden

²⁸ INAF Istituto di Radioastronomia, I-40129 Bologna, Italy

²⁹ Dipartimento di Astronomia, Università di Bologna, I-40127 Bologna, Italy

³⁰ Dipartimento di Fisica, Università di Udine and Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, Gruppo Collegato di Udine, I-33100 Udine

³¹ Università Telematica Pegaso, Piazza Trieste e Trento, 48, I-80132 Napoli, Italy

³² Università di Udine, I-33100 Udine, Italy

³³ Dipartimento di Fisica “M. Merlin” dell’Università e del Politecnico di Bari, I-70126 Bari, Italy

³⁴ Erlangen Centre for Astroparticle Physics, D-91058 Erlangen, Germany

³⁵ Institut für Astro- und Teilchenphysik und Institut für Theoretische Physik, Leopold-Franzens-Universität Innsbruck, A-6020 Innsbruck, Austria

³⁶ Institute of Space and Astronautical Science, Japan Aerospace Exploration Agency, 3-1-1 Yoshinodai, Chuo-ku, Sagamihara, Kanagawa 252-5210, Japan

³⁷ University of North Florida, Department of Physics, 1 UNF Drive, Jacksonville, FL 32224, USA

³⁸ School of Physics and Astronomy, University of Southampton, Highfield, Southampton, SO17 1BJ, UK

³⁹ Yunnan Observatories, Chinese Academy of Sciences, Kunming 650216, China

⁴⁰ Key Laboratory for the Structure and Evolution of Celestial Objects, Chinese Academy of Sciences, Kunming 650216, China

⁴¹ Science Institute, University of Iceland, IS-107 Reykjavik, Iceland

⁴² Department of Physics, Graduate School of Science, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

⁴³ Department of Physics, KTH Royal Institute of Technology, AlbaNova, SE-106 91 Stockholm, Sweden

⁴⁴ Institute of Space Sciences (IEEC-CSIC), Campus UAB, E-08193 Barcelona, Spain
<sup>45</sup> Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA

<sup>46</sup> CNRS, IRAP, F-31028 Toulouse cedex 4, France

<sup>47</sup> Department of Physics and Department of Astronomy, University of Maryland, College Park, MD 20742, USA

<sup>48</sup> Hiroshima Astrophysical Science Center, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8526, Japan
<sup>49</sup> Istituto Nazionale di Fisica Nucleare, Sezione di Roma “Tor Vergata,” I-00133 Roma, Italy

<sup>50</sup> Department of Physics and Astronomy, University of Denver, Denver, CO 80208, USA

<sup>51</sup> Max-Planck-Institut für Physik, D-80805 München, Germany

<sup>52</sup> Department of Physics, University of Johannesburg, P.O. Box 524, Auckland Park 2006, South Africa

<sup>53</sup> National Research Council Research Associate, National Academy of Sciences, Washington, DC 20001, resident at Naval Research Laboratory, Washington, DC 20375, USA

<sup>54</sup> NYCB Real-Time Computing Inc., Lattingtown, NY 11560-1025, USA

<sup>55</sup> Max-Planck Institut für extraterrestrische Physik, D-85748 Garching, Germany

<sup>56</sup> Department of Chemistry and Physics, Purdue University Calumet, Hammond, IN 46323-2094, USA

<sup>57</sup> Solar-Terrestrial Environment Laboratory, Nagoya University, Nagoya 464-8601, Japan

<sup>58</sup> Department of Physical Sciences, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8526, Japan

<sup>59</sup> Institució Catalana de Recerca i Estudis Avançats (ICREA), Barcelona, Spain

<sup>60</sup> Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, and Università di Trieste, I-34127 Trieste, Italy

<sup>61</sup> Laboratory for Astroparticle Physics, University of Nova Gorica, Vipavska 13, SI-5000 Nova Gorica, Slovenia

Received 2015 May 8; accepted 2016 February 19; published 2016 April 22

## ABSTRACT

Most of the celestial $\gamma$ rays detected by the Large Area Telescope (LAT) on board the *Fermi Gamma-ray Space Telescope* originate from the interstellar medium when energetic cosmic rays interact with interstellar nucleons and photons. Conventional point-source and extended-source studies rely on the modeling of this diffuse emission for accurate characterization. Here, we describe the development of the Galactic Interstellar Emission Model (GIEM), which is the standard adopted by the LAT Collaboration and is publicly available. This model is based on a linear combination of maps for interstellar gas column density in Galactocentric annuli and for the inverse-Compton emission produced in the Galaxy. In the GIEM, we also include large-scale structures like Loop I and the *Fermi* bubbles. The measured gas emissivity spectra confirm that the cosmic-ray proton density decreases with Galactocentric distance beyond 5 kpc from the Galactic Center. The measurements also suggest a softening of the proton spectrum with Galactocentric distance. We observe that the *Fermi* bubbles have boundaries with a shape similar to a catenary at latitudes below 20° and we observe an enhanced emission toward their base extending in the north and south Galactic directions and located within ~4° of the Galactic Center.

**Key words:** gamma rays: diffuse background – gamma rays: general – gamma rays: ISM

## 1. INTRODUCTION

The hypothesis of interstellar $\gamma$-ray emission, also known as diffuse Galactic emission, dates back to the 1950s when Satio Hayakawa theorized the existence of intense photon production resulting from the decay of the newly discovered neutral pion (Hayakawa 1952). Early estimates of the intensity and distribution of this emission, together with the bremsstrahlung radiation of electrons and positrons, inverse-Compton scattering (IC), as well as other secondary mechanisms responsible for the production of interstellar emission were made by Morrison (1958), Pollack & Fazio (1963), Ginzburg & Syrovatskij (1965), and Stecker (1966).

The *Third Orbiting Solar Observatory OSO-3*, launched in 1967, confirmed the existence of Galactic interstellar emission by observing for the first time a correlation between high-energy $\gamma$ rays and Galactic structures (Kraushaar et al. 1972). Then, the *Small Astronomy Satellite 2* (SAS-2, launched in 1972), which collected 20 times more photons, provided clear evidence of a correlation between the distributions of high-energy $\gamma$ rays and of atomic hydrogen (H I). This evidence was quantified by Lebrun & Paul (1979), who compared the SAS-2 sky intensity with the atomic hydrogen column density ($N_{H\text{I}}$). Later, the comparison between the *Cosmic-Ray Satellite COS-B* data and the dust extinction derived from galaxy counts revealed the contribution from the molecular hydrogen gas to

the interstellar emission (Lebrun et al. 1982; Strong et al. 1982). When the wide-latitude Columbia University CO survey (Dame & Thaddeus 1984) became available, Lebrun et al. (1983) used it as a tracer for molecular hydrogen. In Bloemen et al. (1986) and Strong et al. (1988), high-energy interstellar emission was realistically modeled for the first time for the whole Galaxy with the inclusion of IC and a partitioning of the gas column density into four Galactocentric annuli to account for radial variations in cosmic-ray (CR) density.

The work described in the present paper is in part based on a similar template method. In this approach, we do not calculate the intensity of the model components from assumed CR densities and production cross-sections. Instead, we use the spatial correlation between the $\gamma$-ray data and a linear combination of gas and IC maps in order to (i) model the diffuse background for point-source studies and (ii) estimate the $\gamma$-ray emissivity of the gas in different regions across the Galaxy. The intensity associated with each template is determined from a fit to $\gamma$-ray data. Strong & Mattox (1996) successfully applied this method to observations of the Energetic Gamma-Ray Experiment Telescope (EGRET) on board the *Compton Gamma-Ray Observatory*. Launched in April 1991, EGRET provided a higher angular and energy resolution and collected about four times more $\gamma$ rays than *COS-B* above 100 MeV. This method was extended by

<footnote n="62"><sup>62</sup> Wallenberg Academy Fellow.</footnote>

<footnote n="63"><sup>63</sup> NASA Postdoctoral Program Fellow, USA.</footnote>

<footnote n="64"><sup>64</sup> Funded by contract FIRB-2012-RBFR12PM1F from the Italian Ministry of Education, University and Research (MIUR).</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Mollweide projection in Galactic coordinates of accumulated counts maps for SAS-2, COS-B, EGRET (above 50 MeV), and Fermi-LAT (above 360 MeV, 4 years, Clean class events). Regions with enhanced numbers of counts due to a non-uniform exposure time in observations with pointed observations are apparent in panels corresponding to SAS-2, COS-B, and EGRET.</figcaption></figure>

Casandjian & Grenier (2008) to include dark neutral medium (DNM) gas (Grenier et al. 2005b; Ade et al. 2011, 2014b) and to study the influence of the interstellar emission model on the detection of γ-ray sources.

The Large Area Telescope (LAT; Atwood et al. 2009) is the main γ-ray detector of the *Fermi Gamma-ray Space Telescope* (*Fermi*) launched on 2008 June 11. Its pair-production towers collect γ rays in the energy range from 20 MeV to greater than 300 GeV. *Fermi* was operated in all-sky survey mode for most of its first four years of operation, allowing the LAT, with its wide field of view of about 2.4 sr, to image the entire sky every two orbits (or three hours). The survey mode, together with an on-axis effective area of ~8000 cm² and a 68% containment of the point-spread function (PSF) of 0.8 at 1 GeV, make the LAT data well suited for studies of interstellar emission and large-scale structures.

Figure 1 shows four years of LAT data together with observations<fnref n="65" /> from SAS-2, COS-B, and EGRET. At LAT energies, the diffuse γ-ray emission of the Milky Way dominates the sky. It contributes five times more photons above 50 MeV than point sources, half of them originating from within 6° of the Galactic midplane. The diffuse emission is bright and structured, especially at low Galactic latitudes, and is a celestial background/foreground for detecting and characterizing γ-ray point sources. Standard LAT analyses based on model-fitting techniques to study discrete sources of γ rays require an accurate spatial and spectral model for the

Galactic diffuse emission. The LAT Collaboration has previously released two Galactic Interstellar Emission Model (GIEM) versions based on the template approach corresponding to gll_iem_v02.fit<fnref n="66" /> and gal_2yearp7v6_v0.fits<fnref n="67" /> tuned, respectively, to 10 months and 24 months of observations. The template method was also applied to LAT observations for studying the interstellar emission in several dedicated regions (Abdo et al. 2010a; Ackermann et al. 2011b, 2012c, 2012a; Ade et al. 2014b, 2014b). This paper describes the GIEM recommended for point-source analyses of the LAT Pass 7 reprocessed data (P7REP) where events have been reconstructed using updated calibrations for the subsystems of the LAT (Bregeon & Charles 2013).

An alternative method for modeling interstellar emission consists of a priori calculations of the CR density and folding this density with γ-ray production cross-sections. This method was used to derive the official *EGRET* interstellar model (Bertsch et al. 1993). The model was based on the assumption that the CR density distribution follows the density of matter convolved with a two-dimensional Gaussian distribution whose width, representing the matter and CR coupling scale, was left as an adjustable parameter in a fitting procedure applied to *EGRET* observations. In this model, the molecular-hydrogen-to-CO conversion factor ($X_{CO}$) was also left free to vary. Hunter et al. (1997) found good agreement between the model and *EGRET* observations, except for an excess of γ rays

<footnote n="66">⁶⁶ http://fermi.gsfc.nasa.gov/ssc/data/access/lat/ring_for_FSSC_final4.pdf</footnote>

<footnote n="67">⁶⁷ http://fermi.gsfc.nasa.gov/ssc/data/access/lat/Model_details/Pass7_galactic.html</footnote>

<footnote n="65">⁶⁵ ftp://heasarc.gsfc.nasa.gov</footnote>
observed above 1 GeV. This excess is likely of instrumental origin (Stecker et al. 2008) and is not seen with the LAT (Abdo et al. 2009a).

The coupling between CRs and matter assumed in Bertsch et al. (1993) may not capture the details of CR propagation in the Galaxy. It predicted a large contrast in γ-ray emissivity inside and outside of spiral arms, at odds with observations (Digel et al. 2001; Ackermann et al. 2011b). An alternative approach for estimating the CR density uses a CR propagation code like GALPROP<fnref n="68" /> (Strong et al. 2007) to model the distribution of CRs across the Galaxy. This code was used in Strong et al. (2000) and Strong et al. (2004a) to predict the Galactic interstellar γ-ray intensity based on the local CR measurements and assumed distribution of CR sources. GALPROP is widely used by the high-energy astrophysics community; recent advances are described in Vladimirov et al. (2011). Extensive comparison between the interstellar emission detected by LAT and GALPROP predictions was performed by Ackermann et al. (2012d). The authors ran GALPROP with various input parameter sets sampled within realistic ranges. They obtained a grid of models associated with different CR source tracers and with various radial and vertical Galactic boundaries. Unlike the present work, the models were not tuned on LAT data, except for $X_{CO}$. They found that the GALPROP models were broadly consistent with LAT data, but noted an under-prediction of the diffuse γ-ray emission in the inner Galaxy above a few GeV and the need for higher CR flux or gas density in the outer Galaxy (Abdo et al. 2010a).

Reasonable agreement has been obtained between propagation code predictions and observations at various wavelengths (Strong et al. 2007, 2011; Bouchet et al. 2011; Gaggero et al. 2013; Orlando & Strong 2013), but our knowledge of the distribution of CR sources, of injection spectra, of CR diffusion properties, and of γ-ray production cross-sections is not sufficiently accurate for these models to be used to precisely describe the diffuse background for point- and extended-source analysis. This is illustrated, for example, in Figure 6 of Ackermann et al. (2012d). The γ-ray emissivity of the gas, which provides the dominant component of the interstellar emission at the energies considered here, is proportional to the CR density folded with γ-ray production cross-sections. The emissivities are more accurately determined with a fit to the data using the template method than from a priori predictions. Nevertheless, propagation codes like GALPROP are essential for the calculation of the spatial distribution of the IC for which no template exists. This emission component is present in every direction of the sky but is brightest in the inner Galactic plane. It is present in the γ-ray data at large angular scales and it dominates the diffuse emission at energies below ~100 MeV. It is therefore important to include a prediction for its spatial distribution in order to properly fit other large-scale components of the GIEM, such as the atomic gas.

Interstellar emission can also be studied directly through observations without using spatial templates. Selig et al. (2015) partitioned the Fermi-LAT counts map above 0.6 GeV into point-like and diffuse contributions. They observed a soft diffuse component tracing the gas content of the interstellar medium (ISM) and a harder component interpreted as IC.

In Section 2, we detail the basic features of the GIEM and the derivation of its templates from observations at other

wavelengths. The model has free parameters that are fit to a selection of γ rays described in Section 3. The model equation and parameters are given in Section 4. The fitting procedure itself and the interpretation of the derived emissivities are detailed in Section 5. The IC intensities predicted with GALPROP for Galactic electrons and positrons need spectral modification as explained in Section 6. We show how we have modeled the part of the diffuse emission that does not correlate with any of the templates in Section 7. We finally describe the construction of the overall GIEM and compare it to LAT observations in Section 8.

## 2. EMISSION COMPONENTS

High-energy interstellar γ-ray emission is produced through the interaction of energetic CRs with interstellar nucleons and photons. The decay of secondary particles produced in hadron collisions, the IC of the interstellar radiation field (ISRF) by electrons, and their bremsstrahlung emission in the interstellar gas are the main contributors to diffuse Galactic emission. Underlying our modeling efforts is the reasonable assumption that energetic CRs uniformly penetrate all of the gas phases in the ISM. CR transport and interactions in magnetized molecular clouds involve complex focusing and magnetic mirror effects and diffusion on small-scale magnetic fluctuations which may lead to an exclusion of CRs from the clouds, or conversely to their concentration in the clouds (Skilling & Strong 1976; Cesarsky & Volk 1978; Gabici et al. 2007; Everett & Zweibel 2011; Padovani & Galli 2013). These processes modify the CR flux at the low energies relevant for gas ionization, but they should leave the CR flux at energies above a few GeV unchanged. We thus do not expect changes in γ-ray emissivity through the rather diffuse gas phases that hold most of the mass, i.e., the HI-bright to CO-bright phases. Experimentally, no evidence of CR screening or re-acceleration in molecular clouds was observed in the molecular complexes or local regions studied with the template method (Digel et al. 1999, 2001; Abdo et al. 2010a; Ade et al. 2014b), except for the Cygnus region where a “cocoon” of freshly accelerated particles was observed with a limited spatial extension of approximately 2° (Ackermann et al. 2011a). We therefore assume that the diffuse γ-ray intensity at any energy can be modeled as a linear combination of templates of hydrogen column density, assuming a uniform CR distribution within each one, an IC intensity map predicted by GALPROP ($I_{IC}$), a template that partially accounts for the emission from Loop I ($I_{Loop}$), an isotropic intensity that accounts for unresolved extragalactic γ-ray sources and for residual CR contamination in the photon data, a map of the solar and lunar emission, and a map for Earth’s limb emission reconstructed in the tails of the LAT point-spread function. To facilitate comparison with the all-sky survey data, the model also includes point-like and extended sources from a preliminary version of the 3FGL catalog (Ballet & Burnett 2013).

### 2.1. Gas Column Densities

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 H⁺ and H II regions; Heiles & Troland 2003). Helium and heavier elements are considered to be uniformly mixed with the

<footnote n="68">⁶⁸ http://galprop.stanford.edu</footnote>
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_{H1}$) 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_{H1}$. 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 $\delta = -25^\circ$. 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_{H1}$ from the observed brightness temperature $T_B$ with Equation (1):

$$N_{H1} = -1.82 \times 10^{18} T_S \int_{v} \ln \left[ 1 - \frac{T_B(v)}{T_S - T_0} \right] dv \ [\text{cm}^{-2}], \quad (1)$$

where the integral is taken over the velocity range of interest (Ackermann et al. 2012d) and $T_0 = 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 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 \le l \le 270^\circ$ and $|b| < 70^\circ$. Changes in $T_S$ indeed modify the spatial distribution of $N_{H1}$ as seen from Earth, and these changes can be probed by the $\gamma$-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 $\gamma$-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^\circ-2^\circ$ 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.

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^\circ$ 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^\circ$ of longitude from their boundaries and scaled each line of sight to the total $N_{H1}$. 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.

### 2.1.2. Molecular Hydrogen

The molecule H₂, which does not have a permanent dipole moment, cannot be observed in direct emission in its dominantly cold phase. The observation of molecular gas relies on other molecules, and especially on the 2.6 mm $J = 1 \rightarrow 0$ line of ¹²C monoxide (CO), which is the second most abundant molecule in the ISM. The millimeter-wave emission of CO can trace H₂ because the molecular hydrogen is its main collisional partner, and collisions excite its rotational transitions (Wilson et al. 1970; Yang et al. 2010). Despite the large optical thickness of the low-level rotational lines of CO, numerous studies suggest that the H₂ column density is
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Evolution, with the H I spin temperature, T<sub>S</sub>, of the log-likelihood obtained for the best fits of the interstellar model to the LAT data in the anticenter region (90° ≤ l ≤ 270° and |b| < 70°) at all energy ranges. The log-likelihood ratio is given with respect to the very best fit obtained at T<sub>S</sub> = 140 K. The temperature T<sub>S</sub> = 10<sup>5</sup> K is equivalent to the optically thin approximation.</figcaption></figure>

proportional to the velocity-integrated CO brightness tempera-
ture W(CO). This relation was experimentally observed by
comparing the virial masses of various molecular clouds to
their CO luminosities and was interpreted by Solomon et al.
(1987), who inferred that molecular clouds have a rich
substructure of small, optically thick regions of distinct
velocity, conceptually analogous to a mist made of discrete
droplets. This molecular “mist” is optically thin at each velocity
and the CO line intensity is proportional to the total number of
molecular “droplets” along the line of sight. The molecular
hydrogen to CO conversion factor (an assumed proportionality
between the integrated column density in the CO line and the
column density of H<sub>2</sub>) is expressed as X<sub>CO</sub> = N<sub>H<sub>2</sub></sub>/W(CO). We
have obtained the W(CO) spatial distribution from the Center
for Astrophysics composite survey (Dame et al. 2001). It is
composed of a Galactic plane survey with a sampling interval
of 0.°125 and surveys covering all of the large local clouds at
higher latitudes with a sampling interval of 0.°25. We have also
used dedicated observations lying outside of the sampling
boundary of the composite CO survey at northern declinations
(T. Dame 2011, private communication). This CO survey
covers the great majority of the Galactic CO emission (Ade
et al. 2014a). We have derived the integrated intensities using
“moment mask” filtering to enhance the signal-to-noise ratio
(Dame 2011). We have derived Galactocentric annuli from
radial CO velocities in the same manner as for H I (Figure 3) in
order to allow for variations in CR density and in X<sub>CO</sub>. The
innermost CO annulus contains all of the high-velocity CO
emission. We have combined outer annulus 8 with annulus 9,
which lacks measurable CO emission to be reliably fitted. This
is equivalent to assuming that the CO is immersed in a uniform
CR flux, with a constant X<sub>CO</sub>, between 10 and 50 kpc. Parts of
the Aquila Rift molecular cloud were incorrectly attributed to
CO annulus 6 using this procedure; we have also merged this
annulus with CO annulus 7.

The central molecular zone (CMZ) is a massive complex of
giant molecular clouds located in the central region of the

Galaxy (Serabyn & Morris 1996; Ferrière et al. 2007). It
appears to be pervaded by intense magnetic fields (Morris 2014)
and recently hosted a burst of star formation (Yusef-Zadeh
et al. 2009). Because of these unique conditions, we have cut
out the CMZ from our innermost ring and created a dedicated
CMZ column density map in H I and CO. For that purpose, we
selected a contour corresponding to 20 K km s⁻¹ in W(CO) for
longitude -1.°5 < l < 4.°5 and we have assigned all of the
molecular gas inside the contour to the CMZ. We have used the
same contour to extract from the LAB survey an N<sub>H I</sub> map for
the CMZ.

2.1.3. Ionized Hydrogen

There is no direct observational information on the spatial
distribution of the warm ionized medium. He et al. (2013)
studied 68 radio pulsars detected at X-ray energies and
compared the free electron column density given by dispersion
measures to N<sub>H I</sub> along the line of sight as traced by X-ray
extinction; they obtained a ratio of H<sup>+</sup> to H I column density in
the range 0.07–0.14. The Hα emission, a two-particle process
proportional to the integral of the square of the electron density
along the line of sight, suffers from dust absorption at low
latitudes (Dickinson et al. 2003) and from scattering by
interstellar dust (Witt et al. 2010; Seon & Witt 2012). The free–
free emission, also proportional to the square of the electron
density, is not absorbed in the radio, but is difficult to separate
from the synchrotron emission and has contributions from
numerous H II regions requiring careful temperature correction.
Cordes & Lazio (2002) developed a model (NE2001) of the
density distribution of Galactic free electrons based on 1143
dispersion measures of pulsars with known distances. We have
built column density annuli maps for the warm ionized medium
based on NE2001 predictions, but adding this component to the
γ-ray model did not improve the fit to the LAT data. We
performed several tests, like excluding the Galactic plane from
the fit or removing individual NE2001 clumps of ionized gas
that seemed over predicted, but we were not able to improve
the fit likelihood. Paladini et al. (2007) also found that the fit of
infrared dust emission worsened with an H<sup>+</sup> column density
map (N<sub>H</sub><sup>+</sup>) extracted from NE2001, probably because of its
simplified spatial distribution. Figure 7 of Sun et al. (2008)
shows that NE2001 does not reproduce the structures of the
free-free emission from the Wilkinson Microwave Anisotropy
Probe (WMAP) observations. Based on Gaensler et al. (2008),
we used an exponential scale height of 1 kpc to build a simple
N<sub>H</sub><sup>+</sup> map which also failed to improve the fit to the LAT data.
The H<sup>+</sup> template normalization was set to zero by the fit and no
H<sup>+</sup>-related structure appeared in the final residual γ-ray map.
With the present sensitivity of the LAT survey, we did not
detect γ rays specifically originating from the small mass in the
diffuse H<sup>+</sup> layer and we have dropped the N<sub>H</sub><sup>+</sup> map from the
models.

2.1.4. Dark Neutral Medium

The ¹²CO J(1 → 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<sub>CO</sub> also depends on the dynamical characteristics of the molecular cloud. Shetty et al. (2011) went beyond
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Galactocentric annuli of $N_{H1}$ in $10^{20} \text{cm}^{-2}$ (left) and $W(\text{CO})$ in K km s$^{-1}$ (right), displayed in Galactic plate carrée projection with bin size of $0^{\circ}125 \times 0^{\circ}125$. The square root color scaling saturates at $100 \times 10^{20} \text{cm}^{-2}$ for $N_{H1}$ and at $50 \text{K km s}^{-1}$ for $W(\text{CO})$. The Galactocentric boundaries for each annulus are written in each panel.</figcaption></figure>

<table><caption>**Table 1**
Limits of the $N_{H1}$ Galactocentric Annuli in Kiloparsecs</caption><thead><tr><td>Annulus</td><td>r<sub>min</sub> (kpc)</td><td>r<sub>max</sub> (kpc)</td></tr></thead><tbody><tr><td>1</td><td>0.0</td><td>1.5</td></tr><tr><td>2</td><td>1.5</td><td>4.5</td></tr><tr><td>3</td><td>4.5</td><td>5.5</td></tr><tr><td>4</td><td>5.5</td><td>6.5</td></tr><tr><td>5</td><td>6.5</td><td>7.0</td></tr><tr><td>6</td><td>7.0</td><td>8.0</td></tr><tr><td>7</td><td>8.0</td><td>10.0</td></tr><tr><td>8</td><td>10.0</td><td>16.5</td></tr><tr><td>9</td><td>16.5</td><td>50.0</td></tr></tbody></table>

the simple “mist” model and calculated the radiative transfer of the CO line in turbulent clouds in order to investigate the dependence of $X_{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_{CO}$ by about a factor of two. Moreover, as we already mentioned, the column density of H1 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 H1–H2 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 H1 and CO-free or CO-quiet H2.

Neutral gas and large dust grains coexist; observations have shown that the gas-to-dust ratio leads to a mass ratio of $M_{\text{gas}}/M_{\text{dust}} \sim 100$. The dust column density can therefore provide an alternative template for H1 and H2 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 $0.7$ 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_{H1}$ and $W(\text{CO})$. Regions of dust not correleted with $N_{H1}$ and $W(\text{CO})$, which spatially correlate with diffuse γ-ray excesses over H1 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_{H1}$ and $W(\text{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| > 10^{\circ}$; then,
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Mollweide projection in Galactic coordinates of the excess (top) and deficit (bottom) of the dust reddening *E*(B-V) over the best linear combination of H I column density and W(CO) intensity maps. Assuming a uniform dust-to-gas ratio and a uniform dust grain radiative emissivity, the dust reddening map can be used as a template for the total gas column density. The excess map predominantly traces the DNM column density of dense H I and diffuse H$_2$ gas. The strong deficits (<-0.2 mag) along the Galactic plane are interpreted as $N_{H1}$ corrections in regions where $T_S > 140$ K. The smaller deficits that extend to 10° or 20° in latitude are likely due to grain emissivity variations in the warm dust of diffuse cloud envelopes.</figcaption></figure>

we fixed them and fitted outer annuli 8 and 9 in the anticenter region (90° < *l* < 270°); 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.

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_{H1}$ and W(CO). We associate this excess map with the DNM distribution. For a standard $N_{H1}/E(B-V)$ ratio (Bohlin et al. 1978; Casandjian 2015b), they translate to column densities in excess of $10^{21} \text{ cm}^{-2}$ 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$_2$ 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$_2$ 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_{H1}$ is overestimated. The inclusion of these potential corrections to $N_{H1}$ in the γ-ray model does improve the fits to the interstellar γ-ray emission. In this paper, we call this map the “$N_{H1}$ 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 ($\beta$ < 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$_2$ 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.

## 2.2. Galactic Inverse Compton Radiation

While the different gas column density maps offer spatial templates for γ-ray photons originating mainly from $\pi^0$-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 $^5Y^{26}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 reacceleration 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
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Intensity distributions, in Galactic coordinates and arbitrary units, used for the north Polar Spur (from the 408 MHz map) and for uniform patches added to account for regions of Extended Excess Emission while deriving the gas emissivities in the nearby and outer annuli (see Sections 4 and 5.1). The last two patches correspond to the *Fermi* bubbles and the previous one to the region around the cocoon of hard spectrum cosmic rays in Cygnus X. Those patches were not used in the final interstellar model.</figcaption></figure>

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_{IC_p}$ to denote the intensity of the predicted IC emission.

## 2.3. Loop I and Fermi Bubbles

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).

## 2.4. Point Sources

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*,<fnref n="69" /> 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⁻³. 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.

## 2.5. Extended and Moving Sources

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

<footnote n="69">⁶⁹ The *Fermi* Science Tools analysis package, the LAT γ-ray data, and the Instrument Response Functions (IRFs) are available from the *Fermi* Science Support Center, http://fermi.gsfc.nasa.gov/fssc.</footnote>
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.

## 2.6. Isotropic Intensity

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.

## 2.7. Residual Earth Limb Emission

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} \text{ cm}^{-2} \text{ s}^{-1} \text{ MeV}^{-1} \text{ sr}^{-1}$, where E denotes the γ-ray energy in MeV. Above 200 MeV, its contribution becomes negligible for our analysis.

## 3. γ-RAY DATA SELECTION

For the LAT data, we have used the P7REP Clean class events from the first four years of the mission. The Clean selection has a reduced residual background of misclassified charged particles compared to the source selection (Ackermann et al. 2012b). We have excised time intervals around bright GRBs and solar flares. This time selection exactly matches that for the 3FGL catalog analysis. We have generated the exposure

<figure><img src="image_8.png" /><figcaption>Figure 6. All-sky distribution, in Galactic coordinates and photon intensity, of the residual emission originating from the Earth limb at 100 MeV.</figcaption></figure>

and PSF maps using the P7REP_V10 (Bregeon & Charles 2013) IRFs (see also Sections 4 and 8.2). The event selection excluded photons with zenith angles greater than 100° and times when the rocking angle of the spacecraft was greater than 52° in order to limit contamination from photons produced in the Earth limb. We have binned the LAT photon counts into 14 equal logarithmic intervals from 50 MeV to 50 GeV. Below 50 MeV, the combined effects of the worsening energy resolution and the steep dependence of the effective area on energy as well as the increased Earth limb contamination due to the increased breadth of the PSF make the study of the diffuse emission more difficult. Above 50 GeV, the statistics are too low to accurately separate the numerous emission components, especially in the Galactic plane.

## 4. γ-RAY MODEL

Because of the ISM transparency to γ rays and of the uniform CR penetration in the H I, DNM, and CO-bright gas phases, we have modeled the all-sky LAT photon data as a linear combination of (l, b) maps in Galactic coordinates from the emission components presented above, namely: H I column densities, $N_{H\text{I}}(l, b)$, and W(CO) intensities, $W_{\text{CO}_i}(l, b)$, in annuli of different Galactocentric radius $R_i$; the total, dust-inferred, gas column density in the DNM, $N_{\text{H}}^{\text{DNM}}(l, b)$; dust-inferred corrections to the total H I column densities, $N_{\text{H I}}^{\text{corr}}(l, b)$; an inverse-Compton intensity, $I_{\text{IC}_p}(l, b, E)$, predicted in direction and energy (E) by GALPROP; an isotropic intensity, $I_{\text{iso}}$, for the extragalactic and instrumental backgrounds; the accumulated emissions from the Sun and moon, $I_{\text{Sun\_Moon}}(l, b, E)$; the Earth's limb emission, $I_{\text{limb}}(l, b)$; a set of 3FGL point sources in the $(l_j, b_j)$ directions, and a set of extended sources, $I_{\text{ext}_{l_j}}(l, b)$.

For a given map pixel and energy band, we have calculated the predicted number of photon counts, $N_{\text{pred}}(l, b, E)$, detected by the LAT as

$$N_{\text{pred}}(l, b, E) = \sum_{i=\text{annulus}} q_{\text{H}\text{I}}(E) \tilde{N}_{\text{H}\text{I}}(l, b) + \sum_{i=\text{annulus}} q_{\text{CO}_i}(E) \tilde{W}_{\text{CO}_i}(l, b) + q_{\text{DNM}}(E) \tilde{N}_{\text{H}}^{\text{DNM}}(l, b) + q_{\text{H I corr}}(E) \tilde{N}_{\text{H I}}^{\text{corr}}(l, b)$$
$$ 
\begin{aligned}
&+ C_{IC_p}(E) \tilde{I}_{IC_p}(l, b, E) + C_{iso}(E) \tilde{I}_{iso} \\
&+ C_{limb}(E) \tilde{I}_{limb}(l, b) + \tilde{I}_{\text{sun\_Moon}}(l, b, E) \\
&+ \sum_{j=\text{point src}} N_{pt_j}(E) \delta(l - l_j, b - b_j) \\
&+ \sum_{k=\text{extend src}} C_{\text{ext}_k}(E) \tilde{I}_{\text{ext}_k}(l, b).
\end{aligned} 
\quad (2) $$

In Equation (2), the *q* parameters represent the γ-ray emissivity of the hydrogen in the associated column density maps and in the CO-bright phase for a given *X*<sub>CO</sub> conversion factor. Without distance information, the *q*<sub>DNM</sub> and *q*<sub>H<sub>1</sub>corr</sub> parameters correspond to the assumption of uniform CR densities in the whole DNM and for all *N*<sub>H<sub>1</sub></sub> corrections. The *C*<sub>IC<sub>p</sub></sub>, *C*<sub>limb</sub>, *C*<sub>iso</sub>, and *C*<sub>ext<sub>k</sub></sub> terms represent normalization factors for the associated intensity maps. *N*<sub>pt<sub>j</sub></sub> denotes the total number of emitted photons per each energy band for each point source represented by the Dirac δ function. All of the *q*, *C*, and *N*<sub>pt<sub>j</sub></sub> factors are to be determined by fits to the LAT data in each energy band. The free *I*<sub>IC<sub>p</sub></sub> normalization partially allows for possible variations of the CR and ISRF spatial and energy distributions from what was assumed. The solar and lunar intensities were not allowed to vary in the fits.

We use the tilde notation $\tilde{I}$ to denote count maps resulting from the convolution with the LAT PSF of the product of an intensity map *I* and of the instrument exposure and pixel solid angle. We have used the Science Tools *gtpsf* (see footnote 69) and *gtexpcube2* (see footnote 69) to estimate the PSF and the binned exposure with the preliminary set of IRFs P7REP_CLEAN_V10. The final model was ultimately scaled to the publicly available P7REP_CLEAN_V15, see Section 8. Since both exposure and PSF are energy dependent, we have applied Equation (2) in energy bins in which we have averaged the PSF and exposure with a weight corresponding to a power-law spectrum of index 2. This choice has very little impact because the energy bands are narrow.

In order to fit the model to the LAT data, we have used a binned maximum likelihood with Poisson statistics. All of the maps were binned in HEALPix<sup>70</sup> with an *N*<sub>side</sub> value of 256, and so the bin size is about 0.25 × 0.25. We have used the code MINUIT<sup>71</sup> to maximize the logarithm of the likelihood and to calculate the uncertainties on the parameters.

We did not fit Equation (2) to the entire all-sky LAT data set at once. We have applied latitude and longitude cuts to define sub-regions where some templates are prominent. This allows us to account for the increasing level of degeneracy between components with decreasing latitude to optimize the derivation of local versus distant emissivities, and to separate the contributions of structured (gas) versus smooth (e.g., *IC*<sub>p</sub>, *I*<sub>iso</sub>) components at the largest angular scales. We describe them with the results presented in the following sections.

The model includes a comprehensive list of diffuse emission components and of known localized sources, but earlier GIEM versions and the results of a preliminary fit with all of the components left free have revealed extensive regions of highly significant emission in excess of the model. Some excesses exhibit patterns relating them to well-known objects such as the Fermi bubbles and Loop I (along the NPS, but also in lower-latitude parts of the old supernova remnant). Other bright excesses of unknown origin extend along the Galactic plane, in

particular, in the first Galactic quadrant at 10° ≤ l ≤ 50°, in the fourth quadrant around l = 315°, and in the Cygnus region. All of these excesses are due to the lack of suitable templates in the model and are further discussed in Section 7. Hereafter, we generically refer to them as regions of extended excess emissions (EEE). To lower the impact of their existence on the derivation of the modeled components, we have developed specific strategies based either on the addition of uniform patches in the model or on the iterative insertion of residual maps, filtered to remove small angular structures. Because they are statistically highly significant, it is necessary to delineate and account for these diffuse excesses to avoid biasing the spectra of the other components in the model, which would otherwise somewhat compensate for the missing features. We have preferred this approach to masking out large excess zones which could jeopardize the derivation of the other parameters.

## 5. GAS EMISSIVITIES

### 5.1. Beyond 7 kpc in Galactocentric Radius

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_1$ emissivity in the 8–10 kpc annulus about the Solar circle, $q_{H_1}(E)$: all longitudes and 10° < |*b*| < 70°;

2. CO emissivities about the Solar circle, $q_{CO_{6+7}}(E)$: all longitudes, 4° < |*b*| < 70°, where CO emission is significantly detected in the moment-masked filtered maps;

3. gas emissivity in the DNM, $q_{DNM}(E)$, $H_1$ emissivity in the 7–8 kpc annulus, $q_{H_{16}}(E)$, and emissivity associated with the $N_{H_1}$ corrections: all longitudes and 3° < |*b*| < 70°;

4. $H_1$ and CO emissivities in the outer Galaxy, $q_{H_{18}}(E)$, $q_{H_{19}}(E)$, and $q_{CO_{8+9}}(E)$: all latitudes and 90° < *l* < 270°.

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

<footnote n="70"><sup>70</sup> http://healpix.sourceforge.net</footnote>

<footnote n="71"><sup>71</sup> https://wwwasdoc.web.cern.ch/wwwasdoc/minuit/minmain.html</footnote>
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.

5.2. In the Inner Galaxy

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 latitudes. 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 ICp 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 NH1 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 XCO = N<sub>H₂</sub>/W(CO) conversion factors for the
inner annuli 1–5, assuming that the γ-ray emissivity of the H₂
molecule is twice that of H₁. Under this assumption the XCO
factor is half the ratio between the H₁ and CO emissivities:
XCOᵢ = qCOᵢ/(2qH₁). We have fitted Equation (2) to the LAT
observations (without patches) in the fourth Galactic quadrant
and GC region (270° < l < 365°). We have left all of the
components free to vary, including the local annulus, which
could be constrained by the latitude extent, |b| < 20°, 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 XCO 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_{\text{CO}} = 0.5 \times 10^{20} \text{ cm}^{-2} (\text{K km s}^{-1})^{-1}$ for annuli 1 and 2,
and $X_{\text{CO}} = 1.5 \times 10^{20} \text{ cm}^{-2} (\text{K km s}^{-1})^{-1}$ 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° ≤ l < 360° and |b| < 20°, with the combined $N_H = N_{H1} + 2X_{CO}W(\text{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.

5.3. Gas Emissivity Spectra Across the Galaxy

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 $\frac{dq}{dE} = q/\Delta E$, where $\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_{\text{CO}} = 0.5 \times 10^{20} \text{ cm}^{-2} (\text{K km s}^{-1})^{-1}$, which is the same as what we found for annuli 1 and 2. The low X_{\text{CO}} value of $0.2–0.7 \times 10^{20} \text{ cm}^{-2} (\text{K km s}^{-1})^{-1}$ 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.

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 log(E) to 5% in the range between 316 MeV and 10 GeV, and a linear increase in log(E) up to 15% at 1 TeV.<fnref n="72" /> Those values represent a lower limit for the hydrogen emissivities in the outer Galaxy where the non-uniformity of T<sub>S</sub> cannot be neglected.

For the hydrogen emissivities in the inner annuli (<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

<footnote n="72">⁷² http://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a)–(j): spectral energy distributions of the γ-ray emissivity per H atom in the H I and H$_2$ phases for the CMZ and the nine Galactocentric annuli. The solid curve shows the best fit obtained with a combination of pion emission from CR nuclei and bremsstrahlung radiation from CR electrons. The dashed curve shows the best fit for the local annulus. To display the gas SED in the DNM (k) and that associated with the $N_{H1}$ correction map (l), we have used a gas-to-dust reddening ratio of $3.5 \times 10^{21} \text{ cm}^{-2} \text{ mag}^{-1}$. We did not display emissivities below $10^{-25} \text{ MeV}^2 \text{ s}^{-1} \text{ sr}^{-1} \text{ MeV}^{-1}$ or the values for the lowest energy bin for the inner-Galaxy annuli. Those points were not used in the analysis.</figcaption></figure>

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_{H1}$ 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_{H1}$ 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 $\pi^0$ 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_1$ and $P_2$ are free spectral indices (Shikaze et al. 2007). This form tends to a power law with index $P_2$ at high rigidity (energy) and the $P_1$ 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_{\text{kin}}/E_0 + ((E_{\text{kin}} - E_4)/E_0)^{-0.5})^{P_3}$, where $E_{\text{kin}}$ is the kinetic energy of the electrons, $E_0 = 1$ GeV, B is a free normalization, $P_3$ is a free spectral index, and $E_4$ 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_1$, 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_1$ across the Galaxy. As shown in Figure 7, we obtained a reasonable fit to the SEDs. The derived values of the proton index $P_2$ and normalization A are fairly insensitive to the approximations on $P_1$ 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_2$ 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 \times 10^{21} \text{ cm}^{-2} \text{ mag}^{-1}$ (Grenier et al. 2005b). In the case of the $N_{H1}$ correction template, we used the same $N_{H1}/E(B-V)$ ratio and left all of the spectral parameters for electrons and protons free to obtain a better fit to the data.

## 5.4. Gradients of Cosmic-Ray Spectra Across the Galaxy

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).

The proton density profile predicted by the GALPROP model $^5Y^2^6R^{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_{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$_2$ interface, however, is more abundant than CO-bright H$_2$ 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_2$ versus the Galactocentric radius of the annulus in Figure 8(c). The index of $2.81 \pm 0.08$ found in the local annulus is fully consistent with the value of $2.82 \pm 0.003$ (stat.) $\pm 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
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Radial distributions across the Galaxy of (a) the γ-ray emissivity per H atom measured at 2 GeV; (b) the proton flux integrated above 10 GV, with the prediction from the GALPROP model $^{5YZ6R30T150C2}$ (solid curve, Ackermann et al. 2012d); (c) the proton spectral index, $P_2$, with statistical error bars and the prediction for proton rigidities above 1 GV from the same GALPROP model (solid line) and from Gaggero et al. (2015; dashed line). In all plots, the horizontal bars span the radial widths of the gas annuli used for the measurements. The two data points with smallest Galactocentric radii have large systematic uncertainties (see text). Panel (d) shows the proton flux integrated above 10 GV, normalized to its value at the Sun Galactocentric radius, with the star formation rate traced by supernova remnants, H II regions, and pulsars (Stahler & Palla 2005).</figcaption></figure>

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 $\delta$ 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.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Spectral evolution of the normalization factors applied to the GALPROP prediction of the Galactic inverse-Compton intensity. The error bars are those obtained from the fit likelihood maximization. The dashed line shows the interpolation used for the construction of the GIEM.</figcaption></figure>

## 6. NORMALIZATION OF THE INVERSE-COMPTON RADIATION

To study the normalization of the model of the IC intensity, with the spatial and spectral distributions predicted by GALPROP (see Section 2.2), we fixed the emissivities of the highly structured components, namely, the H I, CO, and dust-related components, at the values previously measured (see Sections 5.1 and 5.2). For each band, we fit the whole sky with Equation (2), leaving free all of the smooth components with large angular scales, namely, the IC, isotropic, and Earth limb emission, as well as the source fluxes because of their broad effective PSFs at low energies. As before, we iteratively smoothed the positive residuals in each energy bin and added them to the previous iteration residuals until the IC normalization coefficient remained constant. Figure 9 shows the normalization factors obtained for the 14 energy bins. The values close to the one found near 100 MeV and at the highest energies indicate that no major modification of the GALPROP prediction is required. The prediction is off by about a factor of two at intermediate energies. Given the complexity of predicting the leptonic production and propagation, as well as the calculation of the ISRF in the Galaxy, the agreement can be considered to be satisfactory. A more detailed discussion of the comparison between the GALPROP IC predictions with various initial conditions is given in Ackermann et al. (2012d). They also concluded that a greater IC intensity was needed for all of the models they considered, in particular, in the inner Galaxy, either from an increased ISRF, more CR electron sources in the inner regions, or a larger Galactic halo. The present choice of a 6 kpc halo and of a radial distribution for the CR sources strongly peaking near 3 kpc (inside the inner molecular ring) provides larger IC intensities than for broader source distributions encompassing the molecular ring (Ackermann et al. 2012d). To investigate the spectral dependence of the normalization factor, one would need to separate the Cosmic Microwave Background (CMB) and dust parts of the IC prediction in Equation (2) to test whether it requires an increase in the very-far-infrared ISRF with respect to the model of Porter et al. (2008). Further investigation of the spectrum and source distribution of the CR electrons should await more

precise spatial and spectral models of the γ-ray emissions from Loop I and the *Fermi* bubbles.

## 7. REGIONS OF EEE

To display the spatial distributions of EEE across the sky, we have built a photon count map using Equation (2) with the gas emissivities measured in the different annuli (as in Sections 5.1 and 5.2), the IC normalization factors derived in Section 6, the isotropic and Earth limb intensities derived from the local annulus fit, the Sun and Moon intensities, and the fluxes of the point-like and small extended sources from the preliminary version of the 3FGL catalog. The emission associated with the radio map of the NPS was not included in the calculation since it provides a very limited description of the potential emission originating from Loop I. We refer to this model as the “baseline” model. Figure 10 (left column) shows the positive difference between the LAT count map and the count map obtained with this model in three energy bands: 50 MeV–1 GeV, 1–11 GeV, and 11–50 GeV. We did not observe strong negative residuals, except in the direction of the Carina arm tangent where the model greatly over-predicts the observations.

Figure 10 exhibits coherent emission features across the sky. They include the NPS and broader Loop I, which dominate at medium-to-high latitudes at low energy, and the *Fermi* bubbles along with the strong emission toward their base, which are both conspicuous above 1 GeV. Spurs of emission visible at medium northern latitudes toward “the interior” of the NPS (roughly within $-10^{\circ} \le l \le 30^{\circ}$ and $10^{\circ} \le b \le 30^{\circ}$) spatially relate to a structure in the local DNM gas distribution, thereby indicating the need for more gas than described by the dust reddening used in this work or a possible enhancement in the CR flux.

We also observe extended sources broadly distributed along the plane at longitudes less than 50° and to a lesser extent at longitudes around 315°. The origin of these excesses is not known. Part of the excess may be caused at low energy by Loop I in the foreground of the Galactic disk. Strong radio recombination line emission has been detected for longitudes around 30° and 330° by Alves et al. (2012) and Alves et al. (2015), and so the excess of γ radiation could also partly relate to ionized gas. The asymmetry observed between the first and fourth Galactic quadrants below 10 GeV could also have a Galactic IC origin since the GALPROP calculation uses a cylindrical geometry instead of the tilted bar and spiral arms of our Galaxy. Other extended excesses are present at low latitudes along the Galactic plane.

Figure 11 shows a close-up view of the fractional excesses found above the baseline model toward the GC region between 1.7 GeV and 50 GeV. Since the inner Galaxy hosts many point sources, we show the excesses with and without the point-source contribution in the model. The *Fermi* bubbles are clearly visible in Figures 11(a)–(c). We have fit the edges of the bubbles within 20° of the GC using various mathematical curves. The edges of the *Fermi* bubbles are well reproduced by two catenary curves: $10^{2.5} \times (\cosh((l-1^{\circ})/10^{2.5}) - 1)$ for the Northern bubble and $-8^{2.7} \times (\cosh((l+1^{\circ})/8^{2.7}) - 1)$ for the Southern one. In Casandjian (2015a), we noted that those catenary curves also correctly reproduce the structures observed close to the GC in the ROSAT X-ray observations (Bland-Hawthorn & Cohen 2003).

We also observe an extended excess of photons close to the GC at the base of the *Fermi* bubbles. This feature is oriented
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Left column: mollweide projection in Galactic coordinates of the EEE found in the energy bands 50 MeV–1 GeV (top), 1–11 GeV (middle), and 11–50 GeV (bottom): the *Fermi*-LAT count maps have been obtained after subtraction of the baseline interstellar model described in Section 7 and they have been smoothed with a two-dimensional symmetric Gaussian of 3° FWHM. Right column: photon specific intensity, at energies 204 MeV (top), 3.4 GeV (middle), and 22 GeV (bottom), of *I*<sub>EEE</sub> that has been developed to describe the EEE at angular scales larger than 2°. All the maps are displayed with a square root scaling and a pixel size of 0.25.</figcaption></figure>

nearly perpendicular to the Galactic plane and extends ~4° from the GC. In this energy range, several studies have reported emission from the GC region in excess of the expectations from standard emission components (Vitale & Morselli 2009; Hooper & Linden 2011; Abazajian & Kaplinghat 2012; Calore et al. 2015). We note that the shape of this excess is very uncertain for latitudes less that 1°.5 due to systematics in the subtraction of foreground emissions.

## 8. THE GIEM FOR THE LAT

To characterize point sources detected by the LAT, a detailed spatial and spectral model of the total diffuse emission visible in their direction is required. We have used the component decomposition of Equation (2) and the results of

its fits to the LAT data to build the GIEM that is publicly available at the *Fermi* Science Support Center (FSSC) website.<fnref n="73" /> The GIEM is meant to describe the total emission originating from the Galactic ISM prior to its detection by the LAT, in other words, not convolved by the instrument response functions. The isotropic spectrum and the intensities due to the Sun, Moon, and Earth limb are not part of the GIEM, but distributed separately on the website because they depend on the photon selection and time span of observations.

### 8.1. Modeling the EEE

In order to build a 3D intensity cube in position and energy, *I*<sub>EEE</sub>(*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

<footnote n="73">⁷³ http://fermi.gsfc.nasa.gov/ssc</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Close-up view of a region within 20° of the GC showing the *Fermi*-LAT count maps integrated between 1.7 GeV and 50 GeV after subtracting the baseline interstellar model described in Section 7, excluding (top row) and including (bottom row) the point-and extended sources from a preliminary 3FGL list in the model. To reduce the emission contrast in latitude, we display the residuals in fractional units (a), dividing the residuals by the model, and in units of standard deviation (b), dividing the residuals by the square root of the model. In (d) we show the residual map after the further subtraction of $I_{\text{EEE}}$; it contains structures smaller than the angular scale included $I_{\text{EEE}}$. The red dashed lines correspond to the catenary functions that reproduce approximately the edge of the *Fermi* bubbles for latitudes below 20° (see text for details). We have smoothed the four maps with a Gaussian of 1° FWHM.</figcaption></figure>
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_{\text{EEE}}$, we rebuilt the LAT count maps in each of the 14 energy bins to a 1°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_{\text{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_{\text{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_{\text{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_{\text{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_{\text{EEE}}$.

## 8.2. GIEM Construction

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_{\text{H}}$), from the Galactic IC emission, and from $I_{\text{EEE}}$ (see Equation (3)). The gaseous components have been scaled according to the emissivity spectra $\frac{dq_{\text{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_{\text{EEE}}$, and so they cannot be studied with this diffuse background model.

$$ I(l, b, E) = \sum_{i=\text{H},\text{I},\text{H}_2,\text{DNM}} \frac{dq_{\text{fit},i}}{dE} (E) N_{\text{H}_i}(l, b) + N_{\text{IC}}(E) I_{\text{IC}_p}(l, b, E) + I_{\text{EEE}}(l, b, E) \quad (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_2$. 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
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Mollweide projection, displayed in log scaling, of the photon specific intensity in the Galactic Interstellar Emission Model at 1 GeV.</figcaption></figure>

compare the LAT count map integrated between 360 MeV and 50 GeV to that predicted by the GIEL, 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σ.

In Figure 14, we present the SEDs of various components of the GIEL, averaged over regions covering the Galactic disk ($|b| < 10°$) 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\alpha$ contribution as an intensity by taking the negative of the photon intensity associated with the $N_{H\alpha}$ 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.

The GIEL 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.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 GIEL in photons $sr^{-1} s^{-1} cm^{-2} MeV^{-1}$. 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 GIEL for the final (V15) IRFs, we have rescaled its intensity by the exposure ratio evaluated for each of the 30 energy planes of the GIEL. The model is then intended for use with the instrument response functions versions

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Top: All-sky Mollweide projection for 4 years of Fermi-LAT γ-ray counts in the 0.36–50 GeV energy band. Middle: counts prediction in the same energy range based on the Galactic Interstellar Emission Model combined with modeled point and extended sources (including the Sun and the moon), the residual Earth limb emission and the isotropic emission. Both maps are displayed with square root scaling to enhance emission away from the plane. Bottom: residual map in units of standard deviations after smoothing with a Gaussian of 2° FWHM. The pixel size for the three maps is 0.25°.</figcaption></figure>

$P7REP_SOURCE_V15$, $P7REP_CLEAN_V15$, and $P7REP_Ultraclean_V15$.

### 8.3. GIEL Accuracy

The GIEL 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.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Spectra of interstellar emission model components for $|b| > 10°$ (upper panel) and $|b| < 10°$ (lower panel). We have decomposed the total intensity (solid line) into emission originating from hydrogen in its different phases: H I (long-dashed), CO (dashed-dotted), DNM (dotted). The emission from IC assuming an axisymmetric ISRF and electron distribution is shown as short dashed lines and the large-scale structures like Loop I and the *Fermi* bubbles are shown as dashed-double-dotted. We show also the negative of the intensity associated with $N_{H\text{I}}$ correction from the negative dust residual as dashed-triple-dotted.</figcaption></figure>

At high latitudes and outside the region covered by $I_{\text{EEE}}$, the uncertainties are likely dominated by the determination of the gas emissivities. The thinness of the local H I is such that $N_{H\text{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_{\text{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_{H\text{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_{H\text{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_{\text{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_{\text{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 $\pm 2\sigma$;

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).

## 9. CONCLUSION

We have constructed a model for Galactic interstellar emission to allow the characterization of γ-ray point and small extended sources in the LAT data with the best precision possible. The model is based on linear combinations of templates spatially correlated with production sites of γ rays. We used H I, CO, and dust reddening maps to describe the γ rays resulting from collisions between CRs and the ISM through hadrons decay and bremsstrahlung emission. The spatial distribution of γ rays resulting from the IC of CRs on the ISRF was calculated by the CR propagation code GALPROP. We determined the intensity associated with each template with a fit to LAT observations in several energy bands. In the first stage of the fit, extended emission like Loop I and the *Fermi* bubbles were accounted for by patches or through iterative procedures. This extended emission was included in the final model by re-injecting LAT residual counts above a baseline model. Those counts were filtered to only include structures with angular scales larger than $2°$. The model is publicly available at the FSSC website.

From this study, we derived the γ-ray emissivity spectrum at various Galactocentric distances. We interpreted those emissivities and observed that the spectrum of CR protons measured in the inner Galaxy is harder than in the outer Galaxy. We derived the radial distribution of the density of CR protons in the Galaxy and found that it shows similarities to the distribution of tracers of massive star formation. In this work, we characterized the shape of the *Fermi* bubbles within $20°$ from the plane and observed a non-centrosymmetric excess of γ rays in the Galactic center above 1 GeV. We also observed strong soft emission in the first and fourth quadrants from an unknown origin.

The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the
LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the KA Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden.

Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France.

Some of the results in this paper have been derived using the HEALPix (Górski et al. 2005) package.

REFERENCES

Abazajian, K. N., & Kaplinghat, M. 2012, *PhRvD*, 86, 083511

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009a, *PhRvL*, 103, 251101

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009b, *PhRvD*, 80, 122004

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010a, *ApJ*, 710, 133

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010b, *PhRvL*, 104, 101101

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011, *ApJ*, 734, 116

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2012, *ApJ*, 758, 140

Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, *ApJS*, 208, 17

Abergel, A., Ade, P. A. R., Aghanim, N., et al. 2011, *A&A*, 536, 21

Abergel, A., Ade, P. A. R., Aghanim, N., et al. 2013, arXiv:1312.1300

Abergel, A., Ade, P. A. R., Aghanim, N., et al. 2014, *A&A*, 566, A55

Acero, F., Ackermann, M., Ajello, M., et al. 2015, *ApJS*, 218, 23

Ackermann, M., Ajello, M., Allafort, A., et al. 2011a, *Sci*, 334, 1103

Ackermann, M., Ajello, M., Baldini, L., et al. 2011b, *ApJ*, 726, 81

Ackermann, M., Ajello, M., Allafort, A., et al. 2012a, *A&A*, 538, 71

Ackermann, M., Ajello, M., Albert, A., et al. 2012b, *ApJS*, 203, 4

Ackermann, M., Ajello, M., Allafort, A., et al. 2012c, *ApJ*, 755, 22

Ackermann, M., Ajello, M., Atwood, W. B., et al. 2012d, *ApJ*, 750, 3

Ackermann, M., Ajello, M., Allafort, A., et al. 2013, *ApJS*, 209, 34

Ackermann, M., Ajello, M., Albert, A., et al. 2015, *ApJ*, 799, 86

Ackermann, M., Albert, A., Atwood, W. B., et al. 2014, *ApJ*, 793, 64

Ade, P. A. R., Aghanim, N., Alves, M. I. R., et al. 2014a, *A&A*, 571, A13

Ade, P. A. R., Aghanim, N., Aniano, G., et al. 2014b, arXiv:1409.3268

Ade, P. A. R., Aghanim, N., Arnaud, M., et al. 2011, *A&A*, 536, 19

Ade, P. A. R., Aghanim, N., Arnaud, M., et al. 2013, *A&A*, 554, A139

Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2011, *Sci*, 332, 69

Aguilar, M., Aisa, D., Alpat, B., et al. 2015, *PhRvL*, 114, 171103

Aguilar-Benitez, M. 1991, *ZPhyC*, 50, 405

Aloisio, R., Blasi, P., & Serpico, P. 2015, *A&A*, 583, A95

Alves, M. I. R., Calabretta, M., Davies, R. D., et al. 2015, *MNRAS*, 450, 2025

Alves, M. I. R., Davies, R. D., Dickinson, C., et al. 2012, *MNRAS*, 422, 2429

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, *ApJ*, 697, 1071

Ballet, J., & Burnett, T. H. 2013, Proc. 33rd ICRC, (arXiv:1501.02003)

Bernard, J.-P., Reach, W. T., Paradis, D., et al. 2008, *ApJ*, 136, 919

Bertsch, D. L., Dame, T. M., Fichtel, C. E., et al. 1993, *ApJ*, 416, 587

Bland-Hawthorn, J., & Cohen, M. 2003, *ApJ*, 582, 246

Bloemen, J. B. G. M., Dogiel, V. A., Dorman, V. L., & Ptuskin, V. S. 1993, *A&A*, 267, 372

Bloemen, J. B. G. M., Strong, A. W., Mayer-Hasselwander, H. A., et al. 1986, *A&A*, 154, 25

Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, *ApJ*, 224, 132

Bouchet, L., Strong, A. W., Porter, T. A., et al. 2011, *ApJ*, 739, 29

Bregeon, J., & Charles, E. 2013, arXiv:1304.5456

Breitschwerdt, D., Dogiel, V. A., & Völk, H. J. 2002, *A&A*, 385, 216

Calore, F., Cholis, I., & Weniger, C. 2015, *JCAP*, 03, 038

Casandjian, J. M. 2015a, arXiv:1502.07210

Casandjian, J.-M. 2015b, *ApJ*, 806, 240

Casandjian, J.-M., & Grenier, I. 2009, arXiv:0912.3478

Casandjian, J.-M., & Grenier, I. A. 2008, *A&A*, 489, 849

Cesarsky, C. J., & Volk, H. J. 1978, *A&A*, 70, 367

Clemens, D. P. 1985, *ApJ*, 295, 422

Cordes, J. M., & Lazio, T. J. W. 2002, arXiv:astro-ph/0301598

Dame, T. M. 2011, arXiv:1101.1499D

Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, *ApJ*, 547, 792

Dame, T. M., & Thaddeus, P. 1984, NASCP, 81, 235

Dickey, J. M., McClure-Griffiths, N. M., Gaensler, B. M., & Green, A. J. 2003, *ApJ*, 585, 801

Dickey, J. M., Strasser, S., Gaensler, B. M., et al. 2009, *ApJ*, 693, 1250

Dickinson, C., Davies, R. D., & Davis, R. J. 2003, *MNRAS*, 341, 369

Digel, S. W., Aprile, E., Hunter, S. D., Mukherjee, R., & Xu, F. 1999, *ApJ*, 520, 196

Digel, S. W., Grenier, I. A., Hunter, S. D., Dame, T. M., & Thaddeus, P. 2001, *ApJ*, 555, 12

Dobashi, K., Bernard, J.-P., Hughes, A., et al. 2008, *A&A*, 484, 205

Elad, M., Starck, J.-L., Querre, P., & Donoho, D. 2005, *Appl. Comput. Harmon. Anal.*, 19, 340

Everett, J. E., & Zweibel, E. G. 2011, *ApJ*, 739, 60

Everett, J. E., Zweibel, E. G., Benjamin, R. A., et al. 2008, *ApJ*, 674, 258

Ferrière, K., Gillard, W., & Jean, P. 2007, *A&A*, 467, 611

Gabici, S., Aharonian, F. A., & Blasi, P. 2007, *Ap&SS*, 309, 365

Gaensler, B. M., Madsen, G. J., Chatterjee, S., & Mao, S. A. 2008, *PASAu*, 25, 184

Gaggero, D., Maccione, L., Bernardo, G. D., Evoli, C., & Grasso, D. 2013, *PhRvL*, 111, 021102

Gaggero, D., Urbano, A., Valli, M., & Ullio, P. 2015, *PhRvD*, 91, 083012

Ginzburg, V. L., & Syrovatskij, S. I. 1965, *SSRv*, 4, 267

Glover, S. C. O., & Low, M.-M. M. 2011, *MNRAS*, 412, 337

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, *ApJ*, 622, 759

Gould, R. J. 1969, *PhRv*, 185, 72

Grenier, I. A., Black, J. H., & Strong, A. W. 2015, *ARA&A*, 53, 199

Grenier, I. A., Casandjian, J. M., & Terrier, R. 2005a, Proc. 29th ICRC (Pune), 4, 13

Grenier, I. A., Casandjian, J.-M., & Terrier, R. 2005b, *Sci*, 307, 1292

Haslam, C. G. T., Salter, C. J., Stoffel, H., & Wilson, W. E. 1982, *A&AS*, 47, 1

Hayakawa, S. 1952, *PThPh*, 8, 571

He, C., Ng, C.-Y., & Kaspi, V. M. 2013, *ApJ*, 768, 64

Heiles, C., & Troland, T. H. 2003, *ApJ*, 586, 1067

Hooper, D., & Linden, T. 2011, *PhRvD*, 84, 123005

Hunter, S. D., Bertsch, D. L., Catelli, J. R., et al. 1997, *ApJ*, 481, 205

Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, *A&A*, 440, 775

Kamae, T., Karlsson, N., Mizuno, T., Abe, T., & Koi, T. 2006, *ApJ*, 647, 692

Kanekar, N., Braun, R., & Roy, N. 2011, *ApJL*, 737, L33

Kraushaar, W. L., Clark, G. W., Garmire, G. P., et al. 1972, *ApJ*, 177, 341

Lebrun, F., & Paul, J. A. 1979, Proc. 16th ICRC (Kyoto), 12, 13

Lebrun, F., Paul, J. A., Bignami, G. F., et al. 1982, *A&A*, 107, 390

Lebrun, F., Bennett, K., Bignami, G. F., et al. 1983, *ApJ*, 274, 231

McClure-Griffiths, N. M., Green, J. A., Hill, A. S., et al. 2013, *ApJL*, 770, L4

Mori, M. 2009, *APh*, 31, 341

Morris, M. R. 2014, arXiv:1406.7859

Morrison, P. 1958, *Il Nuovo Cimento*, 7, 858

Murray, C. E., Stanimirović, S., Goss, W. M., et al. 2015, *ApJ*, 804, 89

Nolan, P. L., Abdo, A. A., Ackermann, M., et al. 2012, *ApJS*, 199, 31

Orlando, E., & Strong, A. 2013, *MNRAS*, 436, 2127

Padovani, M., & Galli, D. 2013, Cosmic Rays in Star-Forming Environments, Vol. 34 (Berlin: Springer)

Paladini, R., Montier, L., Giard, M., et al. 2007, *A&A*, 465, 839

Panov, A. D., Adams, J. H., Ahn, H. S., et al. 2006, arXiv:astro-ph/0612377

Paradis, D., Dobashi, K., Shimoikura, T., et al. 2012, *A&A*, 543, 103

Pineda, J. L., Langer, W. D., Velusamy, T., & Goldsmith, P. F. 2013, *A&A*, 554, 103

Pollack, J. B., & Fazio, G. G. 1963, *PhRv*, 131, 2684

Porter, T. A., Moskalenko, I. V., Strong, A. W., Orlando, E., & Bouchet, L. 2008, *ApJ*, 682, 400

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, *ApJ*, 500, 525

Selig, M., Vacca, V., Oppermann, N., & Enßlin, T. A. 2015, *A&A*, 581, A126

Seon, K.-I., & Witt, A. N. 2012, *ApJ*, 758, 109

Serabyn, E., & Morris, M. 1996, *Natur*, 382, 602

Shetty, R., Glover, S. C., Dullemond, C. P., et al. 2011, *MNRAS*, 415, 3253

Shibata, T., Honda, N., & Watanabe, J. 2007, *APh*, 27, 411

Shikaze, Y., Haino, S., Abe, K., et al. 2007, *APh*, 28, 154

Skilling, J., & Strong, A. W. 1976, *A&A*, 53, 253

Sodroski, T. J., Odegard, N., Dwek, E., et al. 1995, *ApJ*, 452, 262

Solomon, P. M., Rivolo, A. R., Barrett, J., & Yahil, A. 1987, *ApJ*, 319, 730

Stahler, S. W., & Palla, F. 2005, The Formation of Stars (Weinheim: Wiley-VCH)

Stecker, F. W. 1966, *SAOSR*, 220

Stecker, F. W. 1970, *Ap&SS*, 6, 377
Stecker, F. W., Hunter, S. D., & Kniffen, D. A. 2008, *APh*, **29**, 25

Stepnik, B., Abergel, A., Bernard, J.-P., et al. 2003, *A&A*, **398**, 551

Strong, A. W., Bignami, G. F., Caraveo, P. A., et al. 1982, *A&A*, **115**, 404

Strong, A. W., Bloemen, J. B. G. M., Dame, T. M., et al. 1988, *A&A*, **207**, 1

Strong, A. W., & Mattox, J. R. 1996, *A&A*, **308**, L21

Strong, A. W., Moskalenko, I. V., & Ptuskin, V. S. 2007, *ARNPS*, **57**, 285

Strong, A. W., Moskalenko, I. V., & Reimer, O. 2000, *ApJ*, **537**, 763

Strong, A. W., Moskalenko, I. V., & Reimer, O. 2004a, *ApJ*, **613**, 956

Strong, A. W., Moskalenko, I. V., Reimer, O., Digel, S., & Diehl, R. 2004b, *A&A*, **422**, L47

Strong, A. W., Orlando, E., & Jaffe, T. R. 2011, *A&A*, **534**, 54

Su, M., Slatyer, T. R., & Finkbeiner, D. P. 2010, *ApJ*, **724**, 1044

Sun, X. H., Reich, W., Waelkens, A., & Enßlin, T. A. 2008, *A&A*, **477**, 573

Sun, Y., Xu, Y., Yang, J., et al. 2015, *ApJL*, **798**, L27

Uhlig, M., Pfrommer, C., Sharma, M., et al. 2012, *MNRAS*, **423**, 2374

Vitale, V., & Morselli, A. 2009, arXiv:0912.3828

Vladimirov, A. E., Digel, S. W., Jóhannesson, G., et al. 2011, *CoPhC*, **182**, 1156

Wilson, R. W., Jefferts, K. B., & Penzias, A. A. 1970, *ApJL*, **161**, L43

Witt, A. N., Gold, B., Barnes, F. S., et al. 2010, *ApJ*, **724**, 1551

Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010, *ApJ*, **716**, 1191

Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, *ApJ*, **718**, 1062

Yoon, Y. S., Ahn, H. S., Allison, P. S., et al. 2011, *ApJ*, **728**, 122

Yusef-Zadeh, F., Hewitt, J. W., Arendt, R. G., et al. 2009, *ApJ*, **702**, 178

Yusifov, I., & Küçük, I. 2004, *A&A*, **422**, 545