FERMI-LAT OBSERVATIONS OF THE DIFFUSE γ-RAY EMISSION: IMPLICATIONS FOR COSMIC RAYS AND THE INTERSTELLAR MEDIUM

We use the recently released version 54 (ver. 54) of the GALPROP code (Vladimirov et al. 2011a) to create models of the DGE. We limit ourselves to models using diffusive-reacceleration with no convection for a Kolmogorov spectrum of interstellar turbulence as has been successfully used to explain CR data and the EGRET γ-ray sky (Strong et al. 2004b) as well as INTEGRAL data (Porter et al. 2008). For a detailed description of the GALPROP code and the improvements in ver. 54 with respect to earlier versions we refer the reader to the dedicated Web site.⁶¹

The parameter files of the GALPROP models used in this paper are available in the supplementary material to this paper available in the online journal. These give a precise definition of the models used which can be reproduced as required. Note that only one scaling factor for the ISRF is given in each file which is the average of the local and inner scaling factors found from the fit (see Sections 3.4 and 3.5.3).

The collection of models used in this analysis is created using the CR source distributions and propagation parameters as described in Section 3.2 for different sizes for the CR confinement region: for an assumed cylindrical geometry where the Sun is located in the Galactic plane 8.5 kpc from the Galactic center, we use radial boundaries, R_h, of 20 kpc and 30 kpc, and vertical boundaries (halo size), z_h, of 4, 6, 8, and 10 kpc, respectively. In addition, we use two assumptions for the optical depth correction of the H i component (see Section 3.3.1) and also two values for the cut at which dust emission is no longer used to correct the total column density (see Section 3.3.4). This results in a total of 128 models.

Supernova remnants (SNRs) are believed to be the principal sources of CRs. However, their Galactic distribution is not well determined (Case & Bhattacharya 1998; Green 2005). Therefore we consider, in addition to the measured SNR distribution from Case & Bhattacharya (1998)⁶² (hereafter SNR distribution), other tracers of supernovae (SNe) explosions. The pulsar distribution is a prime candidate as a proxy distribution, because pulsars are an SN explosion end state. The distribution of pulsars is also better determined than SNRs. Still, it suffers from observational biases and for that reason we test two different pulsar distributions, one from Yusifov & Küçük (2004, hereafter Yusifov distribution) and another from Lorimer et al. (2006, hereafter Lorimer distribution). One of the main differences between the two distributions is the functional form used to fit the observational points. Both have a maximum between R = 0 and R = R_☉ and fall off exponentially in the outer Galaxy. However, Lorimer et al. (2006) force the source spatial distribution to zero at R = 0, whereas it is non-zero in Yusifov & Küçük (2004). As an additional proxy for the distribution of CR sources, we consider the distribution of OB stars from Bronfman et al. (2000, hereafter OB stars distribution). OB associations are putative CR acceleration regions and these stars are also the progenitors of core collapse SNe that can leave compact objects, such as pulsars. The four CR source distributions used in this paper are plotted in Figure 1.

Zoom In Zoom Out Reset image size

To determine the CR injection spectra and propagation parameters we perform a χ² fit to CR nuclei, electron, and positron data, using the Minuit2⁶³ minimizer. To reduce computation time, the fit is done in two parts. The propagation parameters and CR nuclei injection spectrum are found from a fit to the CR nuclei data first. The electron injection spectrum is then found by fitting to the total CR electron and positron spectrum, including the contribution by secondary electrons and positrons from CR protons and He interacting with the interstellar gas. Fitting the propagation parameters in the first step decreases the computation time because nuclei up to ₁₄Si must be included in the propagation calculation for an accurate B/C ratio determination, while for the secondary electrons and positrons only protons and He are important. Not calculating the secondary electrons and positrons in the fit to the propagation parameters also saves a considerable amount of time. We use the CR database created by Strong & Moskalenko (2009) and use the data sets and parameters as discussed below. When comparing the models to the CR data we account for solar modulation using the force-field approximation (Gleeson & Axford 1968). In addition to the propagation parameters, the modulation potential for each experiment that has data below a few GeV (AMS, BESS, and ACE) is a free parameter during the fit. For Fermi-LAT and HEAO-3 we fixed the modulation at 300 MeV and 600 MeV, respectively, appropriate for the low solar activity during observations by each instrument as observed with the ACE satellite (Wiedenbeck et al. 2005). Solar modulation is unimportant for JACEE data.

3.2.1. Protons and Heavier Nuclei

We assume that the injection spectra for all CR nuclei species are described by the same rigidity-dependent function

$$\begin{equation} n_p(\rho) \propto \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\rho ^{\gamma _{p,1}} & \rho < \rho _p, \\ \rho ^{\gamma _{p,2}} & \rho _p \le \rho, \end{array}\right. \end{equation} \tag{ 1 }$$

where the indices and break rigidities are obtained by tuning the model to the observed spectrum of CR protons as well as the He, C, and O nuclei spectra. The low-energy intensity and spectrum are affected by solar modulation so we use data taken at low periods in the solar activity cycle. The inclusion of nuclei up to O is to ensure the major contributor species to the production of the secondaries B and Be are properly included. For the proton and He spectra we use low-energy data from BESS (Sanuki et al. 2000) and high-energy data from JACEE (Asakimori et al. 1998). For the C and O spectra we use low-energy data from ACE (ACE Team 2005) and high-energy data from HEAO-3 (Engelmann et al. 1990). To determine the diffusion coefficient, D₀, and Alfvén speed, v_A, for an assumed halo size, we use the B/C ratio because it is the most accurately observed secondary-to-primary ratio. Low-energy ACE (Davis et al. 2000) and high-energy HEAO-3 (Engelmann et al. 1990) data are used for this ratio.

The break in the CR proton and He spectrum observed by ATIC-2, CREAM, and PAMELA (Wefel et al. 2008; Ahn et al. 2010; Yoon et al. 2011; Adriani et al. 2011b) is not taken into account in this modeling and neither are the different spectral indices for protons and He. Vladimirov et al. (2011b) explore different scenarios for the break and different indices and find that the γ-ray intensities and spectra for their models are smaller than the systematic uncertainty of the Fermi-LAT effective area.

Because we derive constraints on the halo size from the γ-ray data, the radioactive secondary ratios are not directly used to fix the propagation conditions, as is usually the case in propagation model studies. However, the models are compared to the ¹⁰Be/⁹Be ratio to check the consistency of the constraints derived from the γ-rays. The ¹⁰Be/⁹Be data uncertainties are large enough to allow a halo size range from 4 kpc to 10 kpc (Strong et al. 2007), depending on the assumed propagation model. We keep the source abundances of nuclei fixed to values determined for ACE data (Moskalenko et al. 2008).

3.2.2. Electrons and Positrons

We assume the injection spectrum of primary CR electrons is described by the rigidity-dependent function

$$\begin{equation} n_e(\rho) \propto \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\rho ^{\gamma _{e,1}} & \rho < \rho _{e,1}, \\ \rho ^{\gamma _{e,2}} & \rho _{e,1} \le \rho < \rho _{e,2}, \\ \rho ^{\gamma _{e,3}} & \rho _{e,2} \le \rho, \end{array}\right. \end{equation} \tag{ 2 }$$

and use data from AMS (Aguilar et al. 2002), Fermi-LAT (Abdo et al. 2009c; Ackermann et al. 2010) and H.E.S.S. (Aharonian et al. 2008, 2009) to determine the spectral indices and break rigidities.⁶⁴ Unfortunately, the data are insufficient to constrain the entire parameter set and unphysical values are obtained if all are freely fit. We therefore fix the values of γ_{e, 1} = 1.6 and γ_{e, 3} = 4, allowing only for freedom in ρ_{e, 1}, ρ_{e, 2}, and γ_{e, 2}. The γ_{e, 1} used in this paper is consistent with that employed by Strong et al. (2011) for modeling the Galactic synchrotron emission, while the exact value of the high-energy index does not significantly affect our results. There is a strong correlation between γ_{e, 1} and the modulation potential for the AMS data in the fits. We therefore caution that the derived values for the modulation potential are biased, although they are in reasonable agreement with the values derived by ACE for the same period (Wiedenbeck et al. 2005). The primary electron injection spectrum is dependent on that obtained from the CR nuclei fits through the size of the confinement volume and corresponding propagation parameters, as well as the contribution of secondary CR electrons and positrons produced by the CR nuclei interacting with the interstellar gas.

The overall fit is made to the data as described above. The normalization is essentially determined from the Fermi-LAT total electron and positron data. No attempt is made to fit the increasing positron fraction reported by Adriani et al. (2009). These can be neglected because the contribution from the excess positrons at high energies to the γ-ray emission is small.

The DGE in the energy range considered in this paper has a strong contribution from π⁰-decay emission.⁶⁵ The treatment in GALPROP is described in detail by Moskalenko & Strong (1998). For proton–proton interactions we use the formulation described in Kamae et al. (2006) for the calculation of the production cross sections. The production of pions from interactions of He nuclei with the interstellar hydrogen, as well as from collisions of CRs with interstellar He, are explicitly included, where we assume in this paper an interstellar He/H ratio of 0.11 by number (Strong & Moskalenko 1998). This is slightly higher than the canonical value of 0.1 found by observations of H ii regions (Deharveng et al. 2000), but it is within systematic uncertainty of those observations. We also ignore production of pions from CR and interstellar gas nuclei heavier than He while their contribution could be as high as ∼10% (Mori 2009). It is assumed that the distribution of interstellar He follows that of interstellar hydrogen, detailed below.

3.3.1. Atomic Hydrogen

3.3.2. Molecular Hydrogen

3.3.3. Ionized Hydrogen

3.3.4. Dust as a Tracer of Gas

The use of dust as a tracer of gas for γ-ray studies goes back to Strong et al. (1982) and Strong & Lebrun (1982). Infrared emission from cold interstellar dust is an alternative to surveys of H i and CO emission lines, which may not trace all of the neutral gas due to various reasons (cold/optically thick H i, variations in X_CO, H₂ not traced by CO). An extensive study of this topic with application to EGRET data has been performed by Grenier et al. (2005), where the total hydrogen column density was derived for each pixel using the E(B − V) reddening maps given by Schlegel et al. (1998). This procedure significantly reduced the residual in the DGE modeling of EGRET data (Grenier et al. 2005). The addition of dust as a tracer of gas has also been used successfully in analysis of Fermi-LAT data (Abdo et al. 2010d; Ackermann et al. 2011b).

In this work, we apply a similar procedure as Grenier et al. (2005) and create a map of "excess" dust column density, E(B − V)_res. We obtain a gas-to-dust ratio for both H i and CO using a linear fit of the N(H i) map and W(CO) map to the E(B − V) reddening map of Schlegel et al. (1998). For simplicity, we assume a constant gas-to-dust ratio throughout the Galaxy. To minimize errors in the fitting, we first determine the gas-to-dust ratio for H i (H i ratio) in regions where no CO is observed and then use that to determine the gas-to-dust ratio for CO (CO ratio) in regions rich in CO. Because the quantity of dust traced by E(B − V) cannot be reliably determined in regions with high extinction, we apply a magnitude cut to the E(B − V) map. To gauge the uncertainty involved with this procedure, we consider two values: magnitude cuts at 2 and 5, respectively. Figure 2 shows that the region affected by these cuts is only a narrow strip around the Galactic plane for both values. The gas-to-dust ratio obtained from our procedure depends on the assumed value of the spin temperature T_S and the E(B − V) magnitude cut. Our derived ratios are given in Table 2. The X_CO factor in the table is determined by assuming a constant proton-to-dust ratio as X_CO = H i ratio/(2 × CO ratio).

Zoom In Zoom Out Reset image size

Table 2. The Gas-to-dust Ratio Determined from a Linear Fit to the H i and CO Component

TSa	E(B − V) Cutb	H i Ratioc	CO Ratiod	XCOe
150	2	74.37	19.52	1.91
150	5	73.00	21.87	1.67
100,000	2	61.39	21.13	1.45
100,000	5	59.99	23.78	1.26

Notes. The X_CO is determined from the gas-to-dust ratios under the assumption that the dust-to-proton ratio is the same for both H i and H₂. See the text for details. ^aIn units of K. T_S = 10⁵ K is equivalent to the optically thin assumption. ^bIn units of mag. ^cIn units of 10²⁰ cm⁻² mag⁻¹. ^dIn units of K (km s⁻¹) mag⁻¹. ^eIn units of 10²⁰ cm⁻² (K (km s⁻¹))⁻¹.

Download table as:  ASCIITypeset image

For simplicity, we use the H i gas-to-dust ratio to turn E(B − V)_res into a column density map, N(E(B − V)_res). This should not cause a significant systematic effect because the X_CO values in Table 2 are similar to those found from the γ-ray fits (see Section 4.3), but the dust-reddening map does not contain distance information. Because the gas-related γ-ray emissivity varies throughout the Galaxy, and depends on the incident CR intensity together with the gas density, correct placement of the residual gas that is traced by the reddening map is essential. While the E(B − V)_res map corrects for uncertainties in N(H i) and W(CO) and its density distribution along each LOS should be similar to N(H i) and W(CO), unique assignment from dust to N(H i) or N(H₂) is not possible. It is difficult to use the density distribution of W(CO) for this placement for two reasons: the sky coverage of W(CO) is smaller than E(B − V)_res, and the X_CO factor is susceptible to variations. N(H i) is not limited in these ways. Therefore, we choose to distribute N(E(B − V)_res) proportionally to the density distribution of N(H i) along each LOS.

N(H i) in the optically thin limit provides a robust lower limit on the H i column density. To account for spurious negative residuals in the reddening map we limit the residual such that the sum of E(B − V)_res and the equivalent reddening of N(H i) and W(CO) is never less than the equivalent reddening of N(H i) in the optically thin limit. The equivalent reddening of W(CO) and N(H i) is evaluated using the determined gas-to-dust ratios, implicitly using the X_CO ratio given in Table 2. W(CO) is included in the sum to account for possible variations in the X_CO ratio in the Galaxy, i.e., N(H i) − N(E(B − V)_res) might be less than N(H i) in the optically thin limit because we overestimate X_CO. We further limit the absolute value of the negative residuals to be less than the H i column density for each LOS so no pixels in the reddening-corrected annular column density maps are negative. This last requirement is needed because our method for calculating the expected model counts assumes no negative pixels. The number of pixels affected by these two cuts is a small fraction of the total and does not affect our results significantly.

Note that this method effectively replaces the N(H i) estimate with N(H i)+N(E(B − V)_res) in the regions not affected by the E(B − V) magnitude cut (see Figure 2). As described earlier, this changes the meaning of T_S because it now acts only as a proxy for the gas-to-dust ratio for a large part of the sky.

The Galactic ISRF is the result of emission by stars, and the scattering, absorption, and re-emission of absorbed starlight by dust in the ISM. Because the ISM is not optically thin for the stellar emission due to the interstellar dust, a radiation transport code must be used to model the distribution of low-energy photons throughout the Galaxy. We calculate the ISRF using the FRaNKIE⁶⁶ code (Porter et al. 2008; see Appendix C for more details). The ISRF model we use in this paper (the "maximum metallicity gradient" model from Porter et al. 2008) has an input bolometric stellar luminosity ∼4 × 10¹⁰ L_☉. This is distributed across the stellar components boxy bulge/thin disc/thick disc/halo with fractions ∼0.1/0.7/0.1/0.1. Approximately 20% of the input stellar luminosity is reprocessed by dust and emitted in the infrared.

A major uncertainty with the ISRF model is the overall input stellar luminosity and how it is distributed among the components of the model. Higher input stellar luminosities for a particular component, e.g., the bulge, will increase the CR electron/positron losses via inverse Compton (IC) scattering and hence the overall output in γ-rays approximately over the spatial region where the stellar model component dominates. Estimates available in the literature illustrating the range for the total Galactic stellar luminosity are, e.g., 6.7 × 10¹⁰ L_☉ (Kent et al. 1991) and 2.3 × 10¹⁰ L_☉ (Freudenreich 1998), with different distributions of the total luminosity across the stellar components used in the models of these authors. Also, the metallicity gradient is important for determining the distribution of interstellar dust (see Porter et al. 2008 for the variation due to the range of Galactic metallicity gradients).

Because of these details, the uncertainty in the ISRF can be considerable in regions like the inner Galaxy. A full exploration of the model parameters for the ISRF is beyond the scope of the current work, so we account for the uncertainty in the ISRF by allowing freedom in the IC emission associated with the optical and infrared (IR) components. This is done by separately calculating with GALPROP the contributions to the IC intensity by optical, IR, and cosmic microwave background photons. Because the optical and IR are physically related, we use a common scaling parameter for both components.

Once the parameters of the propagation model have been determined, the predicted γ-ray maps are compared to the Fermi-LAT data. The comparison is non-trivial due to the uncertainties in some of the DGE parameters described above, along with other γ-ray sources emitting in the Fermi-LAT energy range. To account for the uncertainties we perform a maximum-likelihood fit to the data using the GaRDiAn tool described in Appendix A including in the model the detected point sources and an isotropic component described below.

3.5.1. Detected Sources

3.5.2. Instrumental and Extragalactic Backgrounds

3.5.3. Fitting Region Subdivision

Figure 3 shows the CO annuli used in this analysis. There is very little CO emission in the outer Galaxy. To minimize the effects that the bright and complex inner Galaxy has on the determination of the CO scaling parameters in the outer Galaxy, we split the maximum-likelihood fit into regions, separating the inner and outer Galaxy. In addition, we also minimize the effect of the bright Galactic plane when determining the isotropic background by fitting low- and high-latitude regions separately. This subdivision results in three regions: |b| > 8° (local), |b| ⩽ 8° and 80° < l < 280° (outer Galaxy), and |b| ⩽ 8° and l < 80° or l > 280° (inner Galaxy). A latitude cut of b = 8° was chosen because all CO with |b| > 4° is considered to be in the local annulus, with the extra 4° accounting for the extension of the Fermi-LAT PSF, and also to reduce effects of the bright plane for the determination of the isotropic spectrum. We first fit in the local region and determine the scaling parameter for the local CO annulus (8–10 kpc) and the spectrum of the isotropic emission. We also allow freedom in the ISRF scaling parameter because there is significant IC emission in this region, a fraction of which originates in the inner Galaxy where the ISRF is most uncertain. These parameters are then fixed and the fit is performed for the outer Galaxy region to determine the CO scaling parameters there. Finally, we fit the remaining CO scaling parameters in the inner Galaxy region and allow the ISRF scaling parameter to be free in the fit because the fraction of IC emission originating from the inner Galaxy is much higher in this region than the local region. The fluxes and spectra of point sources in each region are fit as described above.

Zoom In Zoom Out Reset image size

To account for the overlap between regions caused by the Fermi-LAT PSF we create a model of the whole sky for each fit, setting non-fitted scaling parameters to their nominal values of 1 and the spectra of the point sources to their values in the 1FGL catalog. This does not affect our results significantly because the overlapping area is a small fraction of the total area for each region, the scaling parameters do not vary significantly from 1, and the point-source fluxes and spectra are close to the 1FGL catalog values.

To account for the effect of the radial variation of X_CO on the CR propagation and the LOS integration when generating the γ-ray sky maps with GALPROP, the above process is iterated using the X_CO distribution found from the γ-ray fit back into the propagation parameter determination/transport calculation. The iteration is done in two steps. First, we calculate for each annulus the average X_CO value, weighted with both the parameterized CO gas distribution used in GALPROP and the integrated CO intensity from the annular maps. These values are then scaled with the values found from the γ-ray fit. To have a smoothly varying X_CO(R), we use power-law interpolation between the scaled values.

We use X_CO(R) = 2 × 10^{19 + 0.1R/(1 kpc)} cm² (K km s⁻¹)⁻¹ for the initial radial variation of X_CO, compatible with the results from Strong et al. (2004c). There is no formal criterion for stopping the process, but we have found that it converges after a few iterations. The results we report in this paper are obtained after four iterations.

The best-fit DGE models to the γ-ray data are determined by comparing their maximum likelihoods (see Appendix A) where higher values are a qualitatively better fit. Figure 4 shows the logarithm of the maximum likelihoods of all the models for the three different fit regions: local, outer Galaxy, and inner Galaxy. No single model stands out as providing the best fit in all three regions simultaneously. The largest difference between models occurs in the outer Galaxy. Because the difference is about three times larger in the outer Galaxy than other regions, the outer Galaxy would dominate in an all-sky likelihood ratio test.

Zoom In Zoom Out Reset image size

While none of the models provides a best fit for all three regions simultaneously, there are patterns in the likelihood results that are similar between regions. The most general trend is that increasing z_h improves the likelihood in all regions, though the effect is strongest for the outer Galaxy. It is also in the outer Galaxy that the difference between models employing different CR source distributions is most pronounced, with the flat SNR distribution favored over the distributions of pulsars and OB stars, which are more peaked in the inner Galaxy. However, this is strongly dependent on z_h and the effect nearly disappears for z_h = 10 kpc. The outer Galaxy also shows an increase in likelihood with larger values of R_h, especially combined with high values of z_h. The models giving the highest CR flux in the outer Galaxy therefore give the largest likelihood. This need for an increased flux in the outer Galaxy compared to standard propagation models has been shown in other Fermi-LAT analyses (Abdo et al. 2010d; Ackermann et al. 2011b).

The value of T_S also has a significant impact on the likelihood values of the models, although the effect differs from region to region. A value of T_S = 150 K is preferred in the outer Galaxy, which is consistent with requiring an increased flux in this region. Lower values of T_S give higher column densities of H i that increase the γ-ray intensity of the models. The effect of T_S is different in the local region, where the optically thin assumption for H i is preferred. As discussed in Section 3.3.4, the H i column density is replaced by that estimated from dust and T_S becomes a proxy for a certain gas-to-dust ratio given in Table 2. A similar consideration applies in the inner Galaxy region, where optically thin H i gives both the maximum and minimum likelihood depending on the value of the adopted E(B − V) cut. The higher cut of 5 mag gives the best fit and thus the E(B − V) column density estimator seems to be preferred even in the inner Galaxy region. The lower gas-to-dust ratio from the optically thin H i assumption is also preferred while the large difference in the likelihood for different cuts of E(B − V) indicate that the optically thin assumption for H i is not appropriate in the Galactic plane as is generally known (see, e.g., Taylor et al. 2003). An E(B − V) cut of 5 mag is also preferred in the outer Galaxy for both values of T_S, showing that E(B − V) is a better total column density tracer than H i and CO combined in the Galactic plane.

While the likelihood ratio test allows comparison between different models, it is not an absolute measure. As we show in Section 4.2, there are large-scale residuals remaining after model subtraction, which indicate missing components in the models that might bias the comparison. However, because the residuals exhibit a spatial structure that is different from the DGE, we do not think there is a strong bias.

While the likelihood ratio test is effective for comparing different models, it is not able to describe the accuracy of each model separately. Examining residual maps and spectra for different sky regions, along with the longitude and latitude profiles, is a direct method for comparison of models with data. Figure 5 shows the counts observed with the Fermi-LAT in the energy range 200 MeV–100 GeV considered in this paper and also the predicted counts from model ^SS^Z4^R20^T150^C5, which we take as our reference model (the use of this as the reference model is not arbitrary because its parameters are similar to the "conventional" model employed in earlier work (Strong et al. 2004b)). This illustrates the general good agreement across the sky between model and data. However, looking in detail reveals discrepancies in particular regions. We discuss these in the following sections. Due to space constraints, we will not show figures for all of the models considered in the paper. A few models are chosen for display, selected to show the range of results, emphasizing the differences between the models. The figures for all of the models are available in the online supplementary material. Note that the comparison models incorporate the factors found from the fit to the γ-ray data so directly comparing the GALPROP output using the GALDEF files provided in the online supplementary material will not give identical results.

Zoom In Zoom Out Reset image size

4.2.1. Residual Sky Maps

Figure 6 shows the residual sky maps in units of standard deviations⁶⁹ for models ^SS^Z4^R20^T150^C5 and ^SL^Z6^R20^T∞^C5. All models display large-scale residuals with similar, but not identical, features. A more physical way of comparing the models to the data are fractional residual maps, (data-model)/data, shown in Figure 7 for the same models. The Galactic plane shows significant (greater than 4σ) positive and negative structure in the inner Galaxy, but mainly positive in the outer Galaxy. While the residuals are statistically significant, Figure 7 shows that the fractional difference in the inner Galactic plane is less than 10%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

All of the models considered have large positive residuals at intermediate and high latitudes about the Galactic center, most notably features coincident with those described by Su et al. (2010) and Dobler et al. (2010), and a feature that is similar to the radio-detected Loop I (Casandjian et al. 2009). The negative residual of the Magellanic stream is also visible in the southern hemisphere. It was not subtracted from the H i annular column density maps because its contribution to the column density was incorrectly assumed to be negligible. However, this does not affect our model comparison because the models all include this same extra column density. Due to the limited freedom in our fits of the DGE to the γ-ray sky, no attempt will be made here to characterize these residual structures but we do note that their shapes depend on the assumed DGE model.

Point sources are also evident in the large-scale residuals, indicating that the point-source fluxes determined by the fit are biased in these areas. However, their PSF-like spatial extent prevents them from affecting the DGE modeling significantly. Only in areas with many overlapping point sources, such as in the Galactic ridge, can they mimic the structure of the DGE. Our tests have shown that inaccurate source modeling causes less than 20% variations in the derived X_CO factors, less than the variation caused by the CR source distribution and gas properties (see Section 4.3).

The track of the Sun along the ecliptic can also be seen (particularly in the north), although it is not very prominent. The quiet Sun is a source of high-energy γ-rays from CR nuclei interacting in its atmosphere (Seckel et al. 1991) and CR electrons and positrons IC scattering of the heliospheric photon field (Moskalenko et al. 2006a; Orlando & Strong 2007, 2008; Abdo et al. 2011). However, when averaged over a year the overall intensity of this component is very small, being less than 5% of the isotropic background over most of the sky (Abdo et al. 2010g), and does not affect the large-scale DGE modeling significantly. The Moon also contributes to the emission from the ecliptic, being nearly as bright as the Sun around 100 MeV. But, the equivalent diffuse intensity from the moving Moon is less due to the inclination of its orbit relative to the ecliptic. The γ-ray spectrum from the Moon is also steeper than that of the Sun (Moskalenko & Porter 2007) and does not contribute at a detectable level above 10 GeV.

For comparisons between models, we calculate the differences between the absolute values of the fractional residual maps for the models. These maps show directly which model better fits the data while the difference between the models might be larger than shown in these plots. To study the effects of individual parameters, we compare models where only a single model grid parameter is varied. In Figure 8 we show the difference residual maps for variation of only the CR source distribution, changes in the size of the CR confinement volume in Figure 9, and variations of the gas properties in Figure 10. While only a single model grid parameter is changed between the models, there are related changes in the propagation parameters, CR source injection spectrum, ISRF scale factor, X_CO factors, isotropic spectrum, and point-source spectra resulting from the CR and γ-ray fits that also affect the results. The variation of the gas-related parameters has the largest and most distributed effect across the sky. Varying the E(B − V) magnitude cut produces differences that are mostly confined to the Galactic plane, but the small change in the gas-to-dust ratio (see Table 2) has an effect at higher latitudes. Changing T_S gives large positive and negative differences over the sky. The most striking feature is toward the outer Galaxy, where changing from T_S = 150 K to the optically thin approximation mostly improves the agreement at intermediate latitudes, but generally worsens it in the outer Galaxy plane. The improvement at intermediate latitudes seems to correlate at least somewhat with the distribution of CO at intermediate latitudes seen in Figure 3, indicating that the gamma-ray signal is sensitive to the ratio of the gas-to-dust ratios for H i and CO. Another explanation might be a Galactocentric gradient of the gas-to-dust ratios, it being higher in the outer Galaxy in agreement with the increased metallicity in the outer Galaxy. Other variations at high latitude are not as strong but still significant. Due to the diffusive propagation of CRs throughout the Galaxy, the steady state CR distribution should not show strong variations on the scales that are needed to account for the differences shown in the figure. This indicates rather that the assumption of a single T_S value, and hence gas-to-dust ratio, should be reconsidered in subsequent work.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While variation of the assumed CR source distribution does not show as strong of an effect as for the gas properties, significant differences are still seen over the sky. No single CR source distribution is best for all regions of the sky. A strong asymmetric feature in the direction of the inner Galaxy can be seen in the top panel of Figure 8, having opposite signs above and below the plane, indicating a missing asymmetry in the model, either in terms of gas properties or CR flux. The outer Galaxy shows similar features, where the intermediate latitudes and the plane have opposite signs. This is most easily seen in the middle panel of Figure 8 where we have an improvement in the fit in the plane but worsening at intermediate latitudes. Because the SNR distribution provides more CR flux in the outer Galaxy, this indicates that the gas distribution could be more closely confined to the plane in the outer Galaxy than estimated in our modeling. This possibility has been studied by Kalberla et al. (2007) who found a reduced extension in z of the gas distribution in the outer Galaxy when assuming the gas in the halo rotated more slowly than gas in the plane.

Variation of the halo size parameters produces low-level residuals, both positive and negative, in different regions of the sky. The halo size, z_h, has the strongest influence in the outer Galaxy and in the region above and below the Galactic plane in the direction of the Galactic center. The former can be explained by increased CR flux in the outer Galaxy when z_h is increased, while the latter is due to increased IC emission in the direction of the Galactic center caused by a longer integration path length along the LOS. The increase in IC emission is suppressed somewhat because the normalization of the ISRF is anti-correlated with z_h (see Section 4.4). Increasing R_h affects only the outer Galaxy significantly, where the models with larger R_h better agree with the data.

The effect of varying a single parameter on the derived residuals can be strongly interdependent on the other adopted input parameters. This is illustrated in Figure 10 where the difference in residuals when varying T_S for two different sets of the other input parameters is shown. The changes are clearly different depending on the input parameters and the resulting parameters found from the CR and γ-ray fits.

Finally, to illustrate the differences between models where more than one parameter is changed, we compare in Figure 11 models ^SL^Z6^R20^T∞^C5, ^SY^Z10^R30^T150^C2, and ^SO^Z8^R30^T∞^C2 to the reference model ^SS^Z4^R20^T150^C5. The models were chosen to illustrate the range of our parameter scan. While some models fare better than others in the likelihood ratio tests (see Figure 4), there is no model that is uniformly better than the other models. There are large-scale positive and negative differences between the models that are distributed over the entire sky, although the greatest differences are at low latitudes. We emphasize that differences between models can be nonlinear, and caution should therefore be exercised when interpreting low-level residual structures because of subtle interplay between the different model parameters.

Zoom In Zoom Out Reset image size

4.2.2. Plots of Spectra

We plot the spectra of models and data for several selected regions. It is evident from Figure 12 that the models give on average a good description of the Fermi-LAT data at high and intermediate latitudes. Even though the likelihoods of the models differ significantly, the model predictions for the total intensity fall within the systematic error of the Fermi-LAT effective area, deviating less than 10% from the data over the entire energy range. This is partly due to the freedom we have when fitting for the isotropic background (see Section 4.5). Figure 13 shows that we overpredict the data in the south polar cap, an indication that the isotropic component is too large and compensating for inaccuracies in the DGE models. But, even for the intermediate-latitude region shown in Figure 14, where the DGE dominates the isotropic component, the agreement is very good.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The models in the current paper agree better with the intermediate-latitude data (Figure 14) than the model presented in Abdo et al. (2009a) for two main reasons. First, we use dust as an additional tracer for gas densities that has been shown to give better results than using only H i and CO tracers (Grenier et al. 2005). This is especially true for intermediate latitudes in the direction toward the inner Galaxy, which is the brightest part of the low intermediate-latitude region. Second, we allow for freedom in both the ISRF scale factor and X_CO to tune the model to the data, which is well motivated given the uncertainty in those input parameters.

The models in general do not fare as well in the Galactic plane where they systematically underpredict the data above a few GeV but overpredict it at energies below a GeV. This is most pronounced in the inner Galaxy (Figure 15), but can also be seen in the outer Galaxy (Figure 16), with even a small excess at intermediate latitudes (Figure 14). Possible explanations for this discrepancy are deferred to the discussion section. We note that the dip in the data visible between 10 and 20 GeV is due to the IRFs used in the present analysis. Figure 17 shows a comparison of model ^SS^Z4^R20^T150^C5 to the data in the outer Galaxy using the Pass 7 clean photons. The dip between 10 and 20 GeV is greatly reduced because of the improved effective area of the new photon class. Because our results do not depend strongly on the exact shapes of the spectra of the data, these improvements in the effective area do not affect our conclusions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The maximum-likelihood trend of preferring models with larger z_h, lower T_S, and flatter CR source distribution (see Figure 4) is illustrated in Figure 16. Going from the SNR distribution to the OB star distribution has very similar effects as changing from T_S = 150 K to the optically thin approximation and also increasing z_h from 4 kpc to 10 kpc. We note that changing z_h with the SNR distribution has a much smaller effect. It is also evident that all of the models underpredict the data in this region above ∼800 MeV, and some even for the entire energy range.

4.2.3. Longitude and Latitude Profile Plots

We compare longitude and latitude profiles of representative models and data for selected regions. For the profile plots, we use three energy bands (200 MeV–1.6 GeV, 1.6–13 GeV, and 13–100 GeV) to increase the statistics in the profiles. Our discussion below is mainly focused on the lowest energy band, because this has the highest statistics and even though the PSF is broader than at higher energies the profiles are wide enough to be relatively unaffected. In general, the models agree well with the data, deviating less than ∼10% from the data over a large fraction of the sky while covering almost two decades of dynamic range in the latitude profiles. From the profile figures, the component associated with CRs interacting with the H i dominates the DGE in most sky regions and for most of the energy range of the Fermi-LAT. The IC component approaches a similar intensity to the H i for high latitudes, and dominates only in the 13–100 GeV energy band. The H₂ component extends only a few degrees from the Galactic plane and is dominant only in the inner Galaxy.

Despite the overall good agreement, the profile residuals do show structure on scales from few degrees to tens of degrees. For the latitude profile in the outer Galaxy shown in Figure 18, it is evident that the models underpredict the data in the Galactic plane, but overpredict it at intermediate latitudes. The exact shape and magnitude of this residual depend on the model. The underprediction in the plane is mostly dependent on the CR flux in the outer Galaxy (CR source distribution and halo height), while the overprediction at intermediate latitudes depends mostly on the assumed T_S value and therefore gas-to-dust ratio (see Section 3.3.4). These effects can be seen also toward the inner Galaxy (Figure 19), but the effect is mostly absent toward the Galactic center (Figure 20). The residual map differences in Figures 8 and 10 also illustrate this.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The dip around the Galactic plane in the residual in Figure 18 is caused by unreasonably large X_CO factors found from the fits (see Section 4.3), artificially increasing the H₂ component. A residual structure coincident with the H₂ component is not seen in any of the other latitude profiles. The underprediction in the outer Galaxy can also be seen in the longitude profiles in the Galactic plane (Figure 21) where peaks in the H₂ component corresponds with dips in the residual. The contribution from detected point sources is also strongest in the plane with a similar profile as the H₂ component, which can also compensate for a lack of freedom in the DGE model during the fitting procedure. The longitude profile in the Galactic plane does not show a correlation of peaks in the source intensity and dips in the residual indicating that sources from the 1FGL catalog are not able to compensate for large-scale inaccuracies in the diffuse emission.

Zoom In Zoom Out Reset image size

All of the latitude profiles display a north–south asymmetry in the residuals, as was shown in the spectra of the polar cap regions in Figure 13. The effect is most noticeable in Figure 19, which is caused mostly by the gas from the Magellanic stream (Mathewson et al. 1974) that was not removed from the H i annular column density maps as mentioned earlier. As the north–south asymmetry is also visible in the outer Galaxy profile where the Magellanic stream has very little effect, there must be some underlying asymmetry. The origin of this asymmetry is not currently known. It is more likely associated with an asymmetry in the CR flux rather than the ISM because the ISM is more observationally constrained.

The model underprediction above a few GeV seen in Figures 15 and 16 is confined to the Galactic plane, as can be seen in Figure 22. The model systematically underpredicts the data in the plane in the 1.6–13 GeV and 13–100 GeV energy bands, but very little excess emission is seen at higher latitudes. This is not seen as clearly in the Galactic center profile (Figure 20) because that region also includes other large-scale residuals, most notably due to features coincident with those described by Su et al. (2010) and Dobler et al. (2010). Note that while these are prominent above 1.6 GeV, they can also be seen at lower energies, but the details of the residual features depend on the DGE model.

Zoom In Zoom Out Reset image size

Figure 21 shows the longitude profile about the Galactic plane for a few different models. It shows how the H i component is affected by different assumptions for T_S, the magnitude cut in the dust map, and the different CR source distributions. The difference in the CR source distribution is also seen in the IC component that is more peaked for the Lorimer source distribution than the SNR distribution. This can be better seen at intermediate latitudes in Figure 23. The effect is noticeable both at intermediate latitudes as well as in the outer Galaxy where CO from the local annulus dominates.

Zoom In Zoom Out Reset image size

The residuals in the plane show signs of small-scale features that are not compatible with statistical fluctuations. Similar residual structure is also seen at intermediate latitudes in Figures 23 and 24, where the most significant structures in the residuals are correlated with peaks in the H i distribution. Note that some peaks in the H i distribution are not associated with residual structure. It is unlikely that the small angular scale fluctuations are due to small-scale CR intensity variations because the bulk of the CR nuclei producing the DGE for the energy range shown are smoothly distributed. The variations are then mostly caused by features in the annular gas maps that introduce artifacts on small angular scales. This suggests that the gas-to-dust ratio is not constant over the sky and can fluctuate by at least 10%. However, comparing the panels in Figure 24, the residual structure can be seen to be energy dependent. The largest variation is toward the inner Galaxy that can be associated with structure coincident with those identified by Su et al. (2010) and Dobler et al. (2010) but smaller variations around l = 100° indicate spectral variations in the CR flux. See, e.g., Bykov & Fleishman (1992) for how OB associations and super-bubbles might have an effect on the CR flux on smaller spatial scales.

Zoom In Zoom Out Reset image size

Figure 25 shows the radial dependence of X_CO for a few selected models. X_CO for all models can be found in the online supplementary material. Our analysis finds that X_CO(R) depends both on the assumed CR source distribution and the gas properties. This is illustrated in Figure 26, which shows X_CO derived for the local annulus for all models. The local X_CO varies by up to a factor of two depending on the value of T_S and the E(B − V) magnitude cut but is nearly independent of the other input parameters. Because the emissivity of the local gas is well determined by observations of CRs, this shows that an accurate determination of the H i column density is important for determining X_CO values from γ-ray data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The scatter in the X_CO is not surprising. The limited freedom in the γ-ray fit can bias the derived values as has already been mentioned. For X_CO the bias can be twofold. For an accurate determination of X_CO we need to know the γ-ray intensity associated with CO as well as the emissivity per H₂ atom. For our case, we calculate the emissivity assuming a CR distribution and estimate the intensity from a fit to the γ-ray data. If the emissivity is incorrectly estimated, the intensity associated with CO will be biased in the opposite direction, resulting in the twofold bias. Methods that determine the emissivity of the gas simultaneously with the intensity associated with CO are independent of this bias, as long as the emissivity is accurately determined from the data (e.g., Abdo et al. 2010d; Ackermann et al. 2011b). Note that the effect of variations in N(H i) applies in all cases. However, the scatter in the X_CO we determine does not significantly affect our comparison between the models except possibly in the inner Galaxy where the molecular gas content is greatest.

Despite the large variations of the X_CO factors between the models there are several features that are consistent. The X_CO factors in the inner Galaxy are systematically higher than the estimate from Strong et al. (2004c), even when using the same CR source distribution. Only in the innermost annulus is there agreement between our X_CO values and those of Strong et al. (2004c). The strong decrease of X_CO in the innermost annulus is consistent for all our models as has also been seen in other analyses (see, e.g., Ferrière et al. 2007). This has been attributed to the breakdown of X_CO as a tracer of H₂ because the ¹²CO line becomes less optically thick in the Galactic center region (Dahmen et al. 1998). There also seems to be a dip in X_CO for the local annulus that was not in Strong et al. (2004c). Our values for the local region X_CO agree very well with the value found for the nearby Gould Belt by Abdo et al. (2010d). Because the X_CO values for this annulus are determined from fits to the local region, the value is associated with high-latitude clouds. This indicates that molecular clouds in the vicinity of the solar system may have different properties than clouds at a similar Galactocentric distance. High-latitude translucent clouds have also been shown before to have lower X_CO values (de Vries et al. 1987; Heithausen & Thaddeus 1990), but more recent work based on other tracers of molecular hydrogen shows that CO intensities might not be linearly related to H₂ column densities in those clouds (Magnani et al. 2003).

There is an exponential increase in the outer Galaxy that is strongly dependent on the assumed CR source distribution, halo size, and T_S. Figure 27 shows the fractional residual (see Section 4.2) for model ^SS^Z4^R20^T150^C5 with the integrated CO emission in the outer Galaxy overlaid. A correlation between the CO emission and negative residuals is evident. On this basis, we conclude that the X_CO values in the outer Galaxy derived in our analysis are strongly biased and we do not show them in Figure 25. This bias is caused by the underprediction of γ-ray intensity in the outer Galaxy by all of the models considered here. Because there is very little CO in the outer Galaxy (see Figure 3) this bias will not strongly affect our results, but only slightly reduce the scatter when comparing the models in the outer Galaxy.

Zoom In Zoom Out Reset image size

The scaling factor of the ISRF is shown in Figure 28 for each model in both the local and inner region. The derived ISRF scaling factor is model dependent and varies with CR source distribution, gas densities, and halo size. In general, the scaling factor is smaller for the pulsar CR source distributions that are more peaked toward the inner Galaxy than both the OB stars and SNR distributions. More CRs will be injected in the inner Galaxy from the pulsar-like source distributions. This produces a corresponding increase in the IC emission in this region because of the larger number of CR electrons and positrons for these source distributions.

Zoom In Zoom Out Reset image size

The ISRF scaling factor is also dependent on T_S and the E(B − V) magnitude cut, indicating that the normalization of the ISRF (and of the IC intensities) serves to compensate for different gas densities in the fits. This is despite the IC component being both spectrally and also spatially different from the H i component. The latitude dependence of the IC component is similar enough to the H i component (see Figure 20) to allow for the correlation between the H i and IC components in the fit. Coupling that with the trend between increased likelihood in the inner region for the optically thin case seems to indicate that a greater intensity of IC is needed for all models, either from an increased ISRF, more CR electron sources, or a larger z_h. The longer confinement time for the CR electrons/positrons for the larger halo sizes, together with the approximately 1/z decrease of the ISRF perpendicular to the Galactic plane (Porter et al. 2008), produces more IC emission for cases with larger z_h (see also Strong et al. 2010). This effect can be seen as a decrease in the scaling factor in the local region with increasing z_h, although the magnitude of the decrease depends on the assumed CR source distribution and gas densities. The ISRF scaling factors follow a trend in that they are larger for the inner Galaxy region for larger z_h. This indicates that further enhancement is required in addition to that provided by the large z_h, which could be due to additional CR electrons/positrons, or an increase in the ISRF in the inner Galaxy region. Disentangling these effects is difficult. One way could be to consider the IC from CMB photons along with the bremsstrahlung emission. Unfortunately, their overall intensity is not high enough to allow them to be used to constrain the CR electron component. The ISRF scaling factors therefore only allow us to infer that the IC modeling requires modifications but a detailed investigation is needed to unravel their origin.

The spectrum of the derived isotropic background is shown in Figure 29 for selected models. Isotropic spectra for all models can be found in the online supplementary material. As for the ISRF normalization and X_CO values, the isotropic background is model dependent. However, the variation in the overall normalization is not very large because it is spatially very different from the other fitted components. The derived spectral shape is also very similar amongst the various models. Both of these indicate that the isotropic background is not strongly biased by variations between the different models.

Zoom In Zoom Out Reset image size

Figure 29 also shows a comparison with the isotropic spectrum from Abdo et al. (2010g) after combining the derived EGB and RIB components from that work. Our estimates of the isotropic background are consistent with Abdo et al. (2010g) below 1 GeV, and systematically higher above this energy. Despite our efforts to minimize the contribution of the Galactic ridge to the determination of the isotropic background by fitting for the local region, Figure 14 shows that the model slightly underestimates the data at low intermediate latitudes above 1 GeV while still being within the systematic uncertainty of the Fermi-LAT data. The fit will compensate for this with the isotropic component as can be seen in Figure 13 where the model overestimates the data in the polar cap regions above 1 GeV. This is especially true in the south polar cap. The difference in the two estimates of the isotropic background can therefore be attributed to the additional freedom allowed in the diffuse Galactic emission model in the analysis of Abdo et al. (2010g) where the intensity spectrum of the local H i annulus and the IC component was allowed to vary freely. The motivations for this additional freedom were the uncertainties associated with the observed CR intensities that are around few percent at energies above 100 GeV reaching more than 10% below 10 GeV caused by the uncertainty in solar modulation, and the size of the CR halo. The combination of these can produce variations both in the H i-related and IC emission. Because our models predict the Fermi-LAT data within the systematic error we do not try to account for this uncertainty in this analysis. The isotropic spectrum from Abdo et al. (2010g) is therefore a better measurement. This does not affect the comparison between the models considered in this paper because they are all treated identically.

Because the main purpose of the fit to CR data is to obtain a propagation model consistent with CR observations, we defer most of the discussion of the results to Appendix D and only summarize the few key points here. We emphasize that all the models give a good representation of the CR data as can be seen in Figures 30 and 31. This has been shown earlier for similar diffusive-reacceleration models (Strong & Moskalenko 1998). But models with z_h = 10 kpc are at the limit of consistency with the observed ¹⁰Be/⁹Be ratio and therefore considering larger values for z_h is not warranted.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The propagation parameters from the fits shown in Figures 32 and 33 are generally in agreement with values found from similar analyses (Strong et al. 2004a; Trotta et al. 2011). The values are also stable between the different models, with notable exceptions for D₀ and v_A. These parameters are strongly correlated with z_h and the properties of the gas distribution, that is, the X_CO profile.⁷⁰ The values of D₀ and v_A are also slightly affected by the assumed CR source distribution. Therefore, the values for these parameters should be interpreted in the context of the entire model that they were determined from because they depend not only on the propagation set up, but also on the distributions of gas and assumed CR sources.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Gamma Ray Diffuse Analysis (GaRDiAn) package is a full sky binned maximum-likelihood analysis tool. It was designed for fitting DGE models to the Fermi-LAT data, although it is general enough to accommodate other instruments. While the framework is capable of fitting nonlinear models, the main emphasis has been on DGE models consisting of a linear combination of template sky maps. At this point, the GaRDiAn package is not publicly available.

The photon data and model are spatially binned on a HEALPix grid (Górski et al. 2005) that is hierarchical, equal-area, and isolatitude allowing for fast spherical harmonics decomposition and nearest neighbor search. As is appropriate for photon-limited data we use a forward folding method for the analysis, turning the model into counts using the IRFs. The GaRDiAn package requires knowledge of the exposure as a function of energy for each direction on the sky and a sky average representation of the PSF as a function of energy. This information needs to be provided in tabulated form, which is standard for the Fermi-LAT IRFs.⁷¹

Denoting the model flux by $f(\mathbf {\theta }, E)$ and the exposure as ${\rm Exp}(\mathbf {\theta }, E)$, we can write the model counts for energy bin i as

$$\begin{equation} F_i(\mathbf {\theta }) = \int _{E_{{\rm min},i}}^{E_{{\rm max},i}} {d}E\, f(\mathbf {\theta }, E) \, {\rm Exp}(\mathbf {\theta }, E). \end{equation} \tag{ A1 }$$

Here, $\mathbf {\theta }$ is a position in the sky, E is the photon energy, and we have assumed that the photon data have been binned in energy bins defined by E_{min, i} ⩽ E ⩽ E_{max, i}. To account for the energy dependence of the PSF, ψ(α, E), we calculate an effective PSF for each energy bin as the spectrally weighted average with the spectra of the model

$$\begin{equation} \Psi _i(\alpha) = \frac{\int _{E_{{\rm min},i}}^{E_{{\rm max},i}} {d}E\, \psi (\alpha, E)\, \langle F\rangle (E)}{\int _{E_{{\rm min},i}}^{E_{{\rm max},i}} {d}E\, \langle F\rangle (E)}. \end{equation} \tag{ A2 }$$

Here, α is the angle between the true photon direction and the reconstructed direction and $\langle F\rangle (E) = \int {d}\Omega \,F(\mathbf {\theta },E)/4\pi$ is the sky average of the model counts as a function of energy. Using the spherical harmonics $Y_{{\rm lm}}(\mathbf {\theta })$, we can write

$$\begin{equation} F_i(\mathbf {\theta }) = \sum _{l=0}^{l_{\rm max}} \sum _{m=-l}^l a_{{\rm lm},i} Y_{{\rm lm}}(\mathbf {\theta }), \end{equation} \tag{ A3 }$$

$$\begin{equation} \Psi _i(\alpha) = \sum _{l=0}^{l_{\rm max}} c_{l0,i} Y_{l0}(\alpha). \end{equation} \tag{ A4 }$$

When using the spherical harmonic decomposition of the PSF, we assume that α is the angle between a point in the sky and the north pole. We also assume that the PSF is azimuthally symmetric so all c_lm = 0 for |m| ⩾ 1. The convolved model is now calculated as

$$\begin{equation} \tilde{F_i}(\mathbf {\theta }) = \sum _{l=0}^{l_{\rm max}} \sum _{m=-l}^l \sqrt{\frac{4\pi }{2l+1}} c_{l0,i} a_{{\rm lm},i} Y_{{\rm lm}}(\mathbf {\theta }). \end{equation} \tag{ A5 }$$

To allow for arbitrary energy binning of the photon data while still handling strong energy dependence of the IRFs, we integrate Equations (A1) and (A2) semi-analytically. We use power-law interpolation of the tabulated input values of the IRFs and model. For the PSF weighting in Equation (A2), we use a single effective power-law index for the entire bin because fine structure within the energy bin is lost in the conversion to counts.

While employing a spherical harmonic decomposition for the convolution of the PSF with the sky maps is extremely efficient, it has limitations. We are limited to using an azimuthally symmetric PSF and must assume the PSF is the same over the entire sky. Fortunately, the tabulated Fermi-LAT PSF is azimuthally symmetric and its variations over the sky are minimal due to both the uniform exposure of the Fermi-LAT in its nominal survey mode operation, and the small variations of the PSF with incident angle.⁷²

Having the model converted to counts and properly convolved, we calculate the likelihood using

$$\begin{equation} L(\mathbf {X}) = \sum _{i,j} D_i(\mathbf {\theta _j})\log (F_i(\mathbf {\theta _j}, X)) - F_i(\mathbf {\theta _j}, X) - \log (D_i(\mathbf {\theta _j})!), \end{equation} \tag{ A6 }$$

where $D_i(\mathbf {\theta _j})$ are the binned photons for energy bin i and HEALPix pixel j, and X are the parameters of the model. The best-fit parameters are found by maximizing the likelihood using Minuit2.⁷³

Under the assumption of uniform circular motion around the Galactic center with rotation curve V(R), the velocity with respect to the local standard of rest of a region with Galactocentric distance R viewed toward direction l, b (in Galactic coordinates) is

$$\begin{equation} v_{\rm LSR} = R_\odot \left(\frac{V(R)}{R} - \frac{V_\odot }{R_\odot }\right) \sin (l) \cos (b). \end{equation} \tag{ B1 }$$

This relation provides a one-to-one relationship between v_LSR and R for any given LOS. We use the parameterized rotation curve of Clemens (1985) using the IAU-recommended values R_☉ = 8.5 kpc for the distance from the Galactic center to the Sun and V_☉ = 220 km s⁻¹ for the velocity of the Sun around the Galactic center (Kerr & Lynden-Bell 1986).⁷⁴ We applied this relation to the 21 cm LAB survey of H i (Kalberla et al. 2005) and the 115 GHz Center for Astrophysics survey of CO (Dame et al. 2001) to transform the spectral measurements into maps of the emission for a range of Galactocentric annuli. The boundaries of the annuli are given in Table 1. The ∼1 kpc width of the annuli is set by the finite non-circular (random and systematic) motions of the gas traced by these surveys as well as internal velocity dispersions of molecular clouds. These non-circular and internal motions limit the practical linear resolution of the velocity-to-distance relation. The outer annuli are broader because the gradient of v_LSR with Galactocentric distance decreases approximately as 1/R beyond the solar circle.

Due to non-circular motion of gas in the Galaxy, a small fraction of the emission has forbidden velocities. This can be due to the v_LSR being greater than the terminal velocity or having an incorrect sign. In our procedure, for the former case the emission is assigned to the tangent point annulus, while for the latter the gas is assigned to the local annulus (i.e., the one that spans R_☉ = 8.5 kpc). In addition, if the gas is placed above a certain height above the Galactic plane, it is assumed to be local. The height differs between the gas distributions and was chosen to be 1 kpc for H i and 0.2 kpc for CO. These values were chosen to be significantly larger than the scale heights of the gas distributions (e.g., Nakanishi & Sofue 2003, 2006).

The kinematic resolution of the method vanishes for directions near the Galactic center and Galactic anti-center. Therefore, we linearly interpolate each annulus independently across the ranges |l| < 10° and |180 − l| < 10° to get an estimate of the radial profile of the gas. To estimate N(H i) or W(CO) at the edge of the region, we calculate the average over a longitude range Δl = 5° on each side of the boundary. The interpolated values are then scaled to match the total N(H i) or W(CO) along each LOS in the regions that were interpolated.

Note that the innermost annulus is entirely enclosed within the interpolated region, necessitating an alternate method to estimate its column density. For H i this is accomplished by assuming the innermost annulus contains 60% more gas than its neighboring annulus. This is a conservative number considering that observations have shown that there is gas depletion in the radial range ∼1.5–3 kpc (see Ferrière et al. 2007 for a review). For CO, we assign all high-velocity emission in the innermost annulus. Here, high velocity means

$$\begin{equation} v_{\rm LSR} < (-50 + 3l)\;{\rm km\,s^{-1}}, \end{equation} \tag{ B2 }$$

and

$$\begin{equation} v_{\rm LSR} > \left\lbrace \begin{array}{@{}l@{\quad }l@{}}25\; {\rm km\,s^{-1}} &\qquad l < 0\\ (10 + 3l)\; {\rm km\,s^{-1}}&\qquad l >= 0. \end{array}\right. \end{equation} \tag{ B3 }$$

These values were found after visual inspection of the CO data. The specific distribution in the innermost 1.5 kpc does not alter the results of this paper in a significant way.

The CO data are from the 115 GHz composite survey of Dame et al. (2001) covering the latitude range |b| < 30°. The coverage is not complete for that range but it is believed that no significant emission is missing. To increase the signal to noise in the data the CO data have been filtered with the moment masking technique (Dame et al. 2001) applied to each component of the survey independently to accurately account for varying noise levels. The sampling grid spacing of the component surveys varies from 0125 to 025, but we rebin to a resolution of 025 for the annuli. This degradation of angular resolution does not affect the DGE analysis significantly for two main reasons. First, the angular resolution of Fermi-LAT below 5 GeV where the majority of photons are detected is larger than 025. Second, the N(H i) annuli are limited anyway to 05 sampling (see below), limiting any gains from better CO sampling to the inner Galactic ridge.

The H i data are from the 21 cm composite LAB survey of Kalberla et al. (2005) covering the entire sky with a 05 sampling. Limited correction has been made for absorption against bright background radio sources and pixels with large negative brightness temperature are replaced with a linear interpolation in longitude between neighboring pixels. Emission from the Small Magellanic Cloud, Large Magellanic Cloud, and Andromeda M31 is excluded from the annuli. The observed brightness temperature, T_B, is converted to column density under the assumption of a uniform spin temperature, T_S, using the equation

$$\begin{equation} N_{{\rm H\,\mathsc{i}}}(v,T_S) = -\log \left(1-\frac{T_B}{T_S-T_{{\rm bg}}}\right) T_S Ci \Delta v, \end{equation} \tag{ B4 }$$

where T_bg ≈ 2.66 K is the brightness temperature of the microwave background at 21 cm and C = 1.83 × 10¹⁸ cm⁻². In cases where T_B > T_S − 5 K, we truncate T_B to T_S − 5 K.

The Galactic ISRF is the result of emission by stars, and the scattering, absorption, and re-emission of absorbed starlight by dust in the ISM. The first calculation considering the broadband (optical to far-IR) spectral energy distribution (SED) as a function of Galactocentric distance was made by Mathis et al. (1983). Subsequently, Chi & Wolfendale (1991) extended the Mathis et al. (1983) calculation, and this work formed the basis for the ISRF model used in the EGRET-team DGE models (e.g., Bertsch et al. 1993). A new calculation of the ISRF was made by Strong et al. (2000), using emissivities based on stellar populations from COBE/DIRBE fits by Freudenreich (1998) and the SKY model of Wainscoat et al. (1992) together with COBE/DIRBE derived infrared emissivities (Sodroski et al. 1997; Dwek et al. 1997). However, no radiative transport was done for the stellar light in the Strong et al. (2000) model, and hence there was no direct coupling between the stellar emission and the output in the IR. A full radiation transport calculation using ray tracing was done by Porter & Strong (2005), which treated the scattering and absorption of the stellar light using a dust model consistent with COBE/DIRBE data, for the first time directly relating the spatial variation of the ISRF SED throughout the Galaxy. Subsequent work (Porter et al. 2008) extended the code to calculate the full angular distribution of the intensity of the ISRF from ultraviolet to far-IR wavelengths, which is essential for the calculation of the anisotropic IC emission (Moskalenko & Strong 2000). Here, we describe further the extension of the code, which has been rewritten to use a parallel Monte Carlo radiative transfer method.

In our model, we represent the stellar distribution by four spatial components: the thin and thick disk, the bulge, and the spheroidal halo. We follow Garwood & Jones (1987) and Wainscoat et al. (1992) and use a table of stellar spectral types comprising main-sequence stars, giants, and exotics to represent the luminosity function (LF) for each of the spatial components. The spectral templates for each stellar type are taken from the semi-empirical library of Pickles (1998). The normalizations per stellar type are obtained by adjusting the space densities to reproduce the observed LFs in the V and K band for the thin disk. The LFs for the other spatial components are obtained by adjusting weights per component for each of the stellar types relative to the normalizations obtained for the thin disk LF. The spatial densities of the thin and thick disk are modeled as exponential disks. For the thin disk available estimates for the radial scale length range from ∼2 kpc to ∼4 kpc while for the thick disk estimates give the range ∼3–4 kpc (e.g., Porcel et al. 1998; Drimmel & Spergel 2001; Jurić et al. 2008; de Jong et al. 2010; McMillan 2011). We use radial scale lengths of 2.5 kpc and 3.5 kpc for the thin and thick disks, respectively, in the present work. The thin disk has a hole interior of a Galactocentric radius of ∼1.7 kpc, following Freudenreich (1998). The scale height of the stellar classes in the thin disk follows Wainscoat et al. (1992), while the thick disk is characterized by a single scale height of 0.75 kpc, which is approximately the middle of values from the literature (e.g., de Jong et al. 2010). The local thick disk to thin disk normalization, ρ_thick(R_☉)/ρ_thin(R_☉), is assumed to be 5%. We take the stellar halo as described by an oblate symmetrical spheroid with axial ratio c/a = 0.7 and power-law density profile ρ_halo∝r^−2.8, intermediate between values found from Sloan data (Jurić et al. 2008; de Jong et al. 2010). The local halo to thin disk normalization, ρ_halo(R_☉)/ρ_thin(R_☉), is 0.5%. The bulge is assumed to be "boxy" following López-Corredoira et al. (2005) with geometrical parameters taken from their paper. For our nominal ISRF, the bulge input luminosity is normalized to the K-band luminosity of Freudenreich (1998) for an assumed R_☉ = 8.5 kpc.

We use a dust model including graphite, polycyclic aromatic hydrocarbons (PAHs), and silicate. Dust grains in the model are spherical and the absorption and scattering efficiencies for graphite, PAHs, and silicate grains are taken from Li & Draine (2001). The dust grain abundance and size distribution are taken from Weingartner & Draine (2001, their best-fit Galactic model). We assume a purely neutral ISM. The absorption and re-emission by the dust is treated as described below in the Monte Carlo method.

Dust follows the Galactic gas distribution and we assume uniform mixing between the two in the ISM (Bohlin et al. 1978). The dust-to-gas ratio scales with the Galactic metallicity gradient. Estimates for the Galactic [O/H] gradient vary in the range 0.04–0.07 dex kpc⁻¹ (Strong et al. 2004c, and references therein). The variation of the metallicity gradient influences the scattering and redistribution of the mainly UV and blue component of the ISRF into the infrared: increased metallicity implies more dust, and therefore increased scattering and absorption of the star light.

The ISRF is calculated for a cylindrical geometry with a maximum radial extent of 50 kpc and maximum height above the Galactic plane of 50 kpc. We use a Monte Carlo method for the photon propagation through the ISM similar to those described by Gordon et al. (2001), Jonsson (2006), and Bianchi (2008). The system volume is segmented into cells with the number of "photons" injected per cell determined according to the ratio of the cell stellar luminosity to the system stellar luminosity for a given number of total injected particles N_total. The parallelization is done using OpenMP⁷⁵ with one thread per cell, injecting all particles for the cell. Each photon emitted within a cell is released with an isotropic angular distribution uniformly over the cell with frequency sampled from the stellar luminosity spectrum at that location. Following emission, the first interaction is forced to increase sampling efficiency (Gordon et al. 2001). The interaction length is sampled according to the dust optical depth in the emitted direction, and the photon is propagated that distance. The interaction is either a scattering or absorption (or the photon is lost from the system if the interaction length is outside the boundary). The probability for a scattering or absorption depends on the frequency dependence of the scattering and absorption cross section for the assumed grain mixture. If the photon is scattered, the new direction of the photon is calculated according to the dust grain type (determined by the relative sizes of the scattering cross sections for the assumed grain mixture) using the Henyey–Greenstein angular distribution function (Henyey & Greenstein 1941). If the photon is absorbed, the grain type that absorbed the photon is determined from the relative sizes of the absorption cross sections of the grain types in the assumed mixture. If the absorbing grain type is "large," that is, it reemits in thermal equilibrium, we use the "temperature correction" method of Bjorkman & Wood (2001) where the absorbing dust at the new location is heated by the photon, raising its temperature. To enforce radiative equilibrium the dust immediately reemits a photon, where the reemitted frequency is chosen from a distribution that corrects the temperature of the dust prior to its new temperature following the absorption of the photon. If the absorbing grain type undergoes stochastic heating (e.g., the nanograin components of the dust mixture) we cannot treat it using this method, and so the amount of luminosity absorbed in the cell is recorded and the photon is removed from the system (to be dealt with in a subsequent step, as described below). The scattered or absorbed/reemitted photons propagate in the system until they either escape, or are fully absorbed. This process is then repeated for each injected stellar photon.

To treat the photons absorbed on the grains undergoing stochastic heating, the frequency-dependent absorbed luminosity by these grains is used to compute the stochastic heating emissivity for each cell using the "thermal continuous" method described by Draine & Li (2001). The procedure used above for injection of stellar photons is then followed, except the stochastic heating emissivity distribution is used in place of the stellar luminosity distribution. With this method, we achieve particle conservation better than ∼10⁻⁴ after a single pass.

We record the intensity distribution of the ISRF at selected points in the Galactic volume used for the Monte Carlo calculation. The intensities as a function of frequency at each location are recorded as HEALPix images (Górski et al. 2005). The low-energy photon number density at each location is the integrated intensity over 4π sr. The ISRF intensity and/or number density at any location in the Galaxy is obtained by linearly interpolating among the sampling positions used in the Monte Carlo simulation.

The full data set for the ISRF used in the current paper is available from the GALPROP Web site (see http://galprop.stanford.edu/code.php for access instructions).

The χ² values resulting from the CR fit are shown in Figure 35. The nuclei part of the fit results in $\chi ^2/{\rm {\rm dof}}\approx 300/131$ for both pulsar source distributions, increasing to $\chi ^2/{\rm {\rm dof}}\approx 340/131$ for the OB star distribution. The largest contribution to the χ² value comes from low-energy protons and high-energy nuclei. Note that the χ² value is not strongly dependent on the halo size and gas model, and all models provide a good representation of the nuclei data that were fitted. This can be seen in Figures 30 and 31 that compare CR observations to a few selected models. Figures for all models are in the online supplementary material.

Zoom In Zoom Out Reset image size

The models incorporated the effect of solar modulation using the force-field approximation (see Section 3.2), with the value of the modulation potential determined in the fit for each of the instrument given in Figure 36. Our results are compatible with other estimates of the modulation potential for the observational epochs considered (Wiedenbeck et al. 2005).

Zoom In Zoom Out Reset image size

The difference in χ² for the different CR source distributions is largely due to differences at low energies, below ∼10 GeV. Note that there is a correlation between the modulation potential of the ACE experiment given in Figure 36 and the nuclei χ² value, and also between the nuclei χ² and γ_{p, 1} shown in Figure 32. While the difference is statistically significant, the use of force-field approximation for the modulation and the constraints of the propagation model do not warrant model selection based on the CR fit. The accuracy of the numerical solution of the propagation equation can also make significant changes to the χ² value. For numerical fitting of the CR spectrum we had to compromise between accuracy and speed. Note that the changes in the model prediction giving rise to the differences in χ² are very small, as demonstrated in Figure 30.

The electron fit is considerably worse, giving $\chi ^2/{\rm {\rm dof}}\approx 650/95$ for z_h = 4 kpc up to $\chi ^2/{\rm {\rm dof}}\approx 850/95$ for z_h = 10 kpc. There is a strong dependence on both halo size and X_CO value,⁷⁶ which is expected because secondary electrons and positrons comprise a significant fraction (≈50%) of the total at low energies. A larger halo also increases the IC cooling of the electrons (Strong et al. 2010). The high χ² values are mostly due to the large fraction of secondaries at low energies, making the AMS and Fermi-LAT spectra incompatible, but also due to the convex shape of the observed Fermi-LAT electron spectrum. This is expected because Ackermann et al. (2010) found that the form of injection spectrum used in the current analysis does not reproduce all the features of the Fermi-LAT total CR electron data. The low-energy part could be improved by increasing the modulation potential for the Fermi-LAT observations, but this is not well motivated given the low level of solar activity during the LAT sky survey observations analyzed here. Because the IC and bremsstrahlung components are subdominant in the energy range where the statistics for studying the DGE are the greatest, this does not affect our results significantly.

Figures 32, 33, and 37 show the resulting parameters from the fit to CR data. Most of these are in agreement with parameters found in earlier studies of CR propagation (Strong et al. 2004b; Trotta et al. 2011). The main difference is the electron injection spectrum, which has now been measured accurately up to 1 TeV by the Fermi-LAT (Abdo et al. 2009c). The spectrum found in our analysis is now more akin to the one used in the optimized model by Strong et al. (2004b) although the break energy at ∼3 GV found in this analysis is much lower than the 20 GV break assumed in the optimized model. The proton spectrum is also slightly different, with a break rigidity at ∼11.5 GV instead of 9 GV in Strong et al. (2004b) and spectral indices of ∼1.9 and ∼2.40 below and above the break, respectively, instead of 1.98 and 2.42 used in the conventional model in Strong et al. (2004b). Our proton spectrum is much closer to the one determined in the more recent analysis by Trotta et al. (2011). Most of these differences can be attributed to a different CR source distribution and the inclusion of X_CO(R) in the propagation calculation. The values of the propagation parameters D₀ and v_A agree reasonably well with the older analysis if one chooses similar models for the comparison, since the value of D₀ and v_A are correlated with z_h and the gas distribution. Note that only the X_CO values were changed for the propagation; the H i distribution is constant. This causes an overestimate of the gas densities in the optically thin assumption, because the analytical model is not corrected for T_S or dust variations, and X_CO compensates to a certain degree for the change in H i column density in the γ-ray fit. This overestimate enhances the variations in D₀ and v_A but does not significantly affect the γ-ray analysis. Note also that the injection spectrum softens with increasing z_h, both for nuclei and electrons, although the variations are barely larger than the error bars. The error bars on the data are statistical only and no attempt has been made to assess the systematic error. The largest systematic uncertainty is likely associated with the absolute energy scale of the data (Ackermann et al. 2010). This can potentially change the location of the break energies along with the intensity of the modeled spectrum. A 10% change in absolute energy would shift the normalization by ∼30%, which directly translates to a change in the intensity of the DGE model.

Zoom In Zoom Out Reset image size