Submitted D Ruse Galactic Continuum Gamma Rays. A model compatible with EGRET data and cosmic-ray measurements Andrew W. Strong Max-Planck-Institut fur extraterrestm’sche Physik, Postfach 1603, 0-85740 Garching, Germany aws@mpe. m pg .d e ----- /---- --__I~ \ Igor V. Moskalenkol I 1 \--- NASA/Goddard Space Flight Center, Code 661, Greenbelt, MD 20771 [email protected] and Olaf Reimer Ruhr- Universitat Bochum, 0-447 80 Bochum, Germany ABSTRACT We present a study of the compatibility of some current models of the diffuse Galactic continuum y-rays with EGRET data. A set of regions sampling the whole sky is chosen to provide a comprehensive range of tests. The range of EGRET data used is extended to 100 GeV. The models are computed with our GALPROP cosmic-ray propagation and gamma-ray production code. We confirm that the “conventional model” based on the locally observed electron and nucleon spectra is inadequate, for all sky regions. A conventional model plus hard sources in the inner Galaxy is also inadequate, since this cannot explain the ‘Joint Center for Astrophysics, University of Maryland, Baltimore County, Baltimore, MD 21250 - 2 - GeV excess away from the Galactic plane. Models with a hard electron injection spectrum are inconsistent with the local spectrum even considering the expected fluctuations; they are also inconsistent with the EGRET data above 10 GeV. We present a new model which fits the spectrum in all sky regions adequately. Secondary antiproton data were used to fix the Galactic average proton spectrum, while the electron spectrum is adjusted using the spectrum of diffuse emission it- self. The derived electron and proton spectra are compatible with those measured locally considering fluctuations due to energy losses, propagation, or possibly de- tails of Galactic structure. This model requires a much less dramatic variation in the electron spectrum than models with a hard electron injection spectrum, and moreover it fits the y-ray spectrum better and to the highest EGRET energies. It gives a good representation of the latitude distribution of the y-ray emission from the plane to the poles, and of the longitude distribution. We show that sec- ondary positrons and electrons make an essential contribution to Galactic diffuse y-ray emission. Subject headings: diffusion - cosmic rays - ISM: general - Galaxy: general - gamma rays: observations gamma rays: theory - 1. Introduction Diffuse continuum y-rays from the interstellar medium are potentially able to reveal much about the sources and propagation of cosmic rays (CR), but in practice the exploitation of this well-known connection is problematic. While the basic processes governing the CR propagation and production of diffuse y-ray emission seem to be well-established, some puzzles remain. In particular, the spectrum of y-rays calculated under the assumption that the proton and electron spectra in the Galaxy resemble those measured locally reveals an excess above 1 GeV in the EGRET spectrum (Hunter et al. 1997). The Galactic diffuse continuum y-rays are produced in energetic interactions of nucleons with gas via neutral pion production, and by electrons via inverse Compton (IC) scattering and bremsstrahlung. These processes are dominant in different parts of the spectrum, and therefore if deciphered the y-ray spectrum can provide information about the large-scale spectra of nucleonic and leptonic components of CR. In turn, having an improved under- standing of the Galactic diffuse y-ray emission and the role of CR is essential for unveiling the spectra of other components of the diffuse emission and is thus of critical importance for the study of many topics in y-ray astronomy, both Galactic and extragalactic. . , -3- The puzzle of the “GeV excess” has lead to an attempt to re-evaluate the reaction of .Iro-production in pp-interactions. However, a calculation made using modern Monte Carlo event generators to simulate high-energy ppcollisions has shown (Mori 1997) that the y- ray flux agrees rather well with previous calculations. A flatter Galactic nucleon spectrum has been suggested as a possible solution to the “GeV excess” problem (Gralewicz et al. 1997; Mori 1997). Explaining the excess requires an ambient proton spectrum power-law index of about -2.4-2.5, compared to -2.75 measured locally (for a summary of recent data see Moskalenko et ai. 2002). Such large variations in the proton spectrum are, however, improbable. A sensitive test of the large-scale-average proton spectrum has been proposed by Moskalenko, Strong, & Reimer (1998) based on the fact that secondary antiprotons and y-rays are produced in the same ppinteractions. The secondary antiprotons sample the proton spectrum in a large region of the Galaxy, and a flatter nucleon spectrum in distant regions would lead to overproduction of secondary antiprotons and positrons. The “hard ~ucleons pectrum” hypothesis has effectively been excluded by recent measurements ‘of p/p ratio at high energies (e.g., Beach et al. 2001). In addition, new accurate measurements of the local proton and Helium spectrum allow less freedom for deviations in the .Iro-decay component. A “hard electron spectrum” hypothesis has been investigated by Porter & Protheroe (1997), Pohl & Esposito (1998), and Aharonian & Atoyan (2000). An essential idea of this approach is that the locally-measured CR spectrum of electrons is not a good constraint because of the spatial fluctuations due to energy losses and the stochastic nature of the sources in space and time; the average interstellar electron spectrum responsible for y-rays via IC emission (and bremsstrahlung) can therefore be quite different from that measured locally. An extensive study of this hypothesis has been made by Strong, Moskalenko, & Reimer (2000); in this model a less dramatic but essential modification of the proton and Helium spectrum (for the TO-decay component) was also required. The latter was still consistent with the locally-observed proton spectrum, as it should be since the proton fluctuations are expected to be small (Strong & Moskalenko 2001a) as the result of their negligible energy losses. The “hard electron spectrum” hypothesis suffers however from the following problems: It is hardly compatible with the local electron spectrum even considering the fluctu- 0 ations due to stochastic sources and energy losses, as shown by a 3D time-dependent study (Strong & Moskalenko 2001b); The fit to the shape of the y-ray spectrum is still poor above 1 GeV (Strong et al. 0 2000); It cannot reproduce the spectrum in the inner and outer Galaxy and intermediate/high 0 latitudes simultaneously (Strong, Moskalenko, k Reimer 2003a); -4- These problems were already evident before the present study, but now we show in addition that it predicts significantly higher intensities than the EGRET data above 10 GeV. 0 Another suggestion which has been made (Berezhko & Volk 2000) is that the y-ray spectrum contains a contribution from accelerated particles confined in SNR. The SNR proton and electron spectra, being much harder than the interstellar CR spectra, produce ro-decay and IC y-rays adding to the apparently diffuse y-rays, while the SNRs themselves are too distant to be resolved into individual sources. A shortcoming of previous analyses was that the comparison with EGRET data was limited to particular regions, and the rich EGRET data have remained not fully exploited. Hunter et al. (1997) made an extensive comparison of the spectra near the Galactic plane Ibl < 10”. Other analyses have concentrated on particular molecular clouds: Ophiuchus (Hunter et al. 1994), Orion (Digel et al. 1999), Cepheus and Perseus (Digel et al. 1996), Monoceros (Digel et al. 2001), high latitudes (Sreekumar et al. 1998). The previous anal- ysis by Strong et al. (2000) was limited to the inner Galaxy at low latitudes, and profiles integrated over large regions of longitude or latitude. In that study, we compared a range of models, based on our CR propagation code GALPROP, with data from the Compton Gamma Ray Observatory. Relative to the work of Hunter et al. (1997) we emphasized the connection with CR propagation theory and the importance of IC emission, and less to fit- ting to structural details of the Galactic plane. The study confirmed that it is rather easy to get agreement within a factor -2 from a few MeV to 10 GeV with a “conventional” set of parameters, however, the data quality warrant considerably better fits. In the present paper we attempt to exploit the fact that the models predict quite spe- cific behaviour for different sky regions and this provides a critical test: the “correct” model should be consistent with the data in all directions. We show that a new model, with less dramatic changes of electron and nucleon spectra relative to the “conventional” model, can well reproduce the y-ray data. The changes consist in renormalization of the intensities of the electron and proton spectra, and a relatively small modification of the proton spectrum at low energies. The model is compatible with locally observed particle spectra considering the expected level of spatial fluctuations in the Galaxy. We extend the y-ray data compar- isons over the entire sky and to 100 GeV in energy. We also exploit the recent improved measurements of the local proton, Helium, as well as antiproton, and positron spectra which are used as constraints on the proton spectrum in distant regions. Our approach differs from that of Hunter et al. (1997) in that it is based on a model of CR propagation while Hunter et al. use CR-gas coupling and a relatively small IC component. . -5- It is also different from Strong & Mattox (1996) in that it is based on a physical model, while that work was based on model-fitting to gas surveys to determine the y-ray emissivity spectrum as a function of Galactocentric radius. The selection of a good model for the diffuse Galactic emission is critical to another topic, the extragalactic ?-ray background (EGRB). We have argued in Strong et al. (2000) and Moskalenko & Strong (2000) that IC from a large halo can make up a substantial fraction of the high-latitude emission and hence reduce the residual EGRB (and modify its spectrum relative to Sreekumar et al. 1998). In a companion paper (Strong, Moskalenko, & Reimer 2004) we present a comprehensive discussion of the EGRB, with a new estimate which is used in the present paper. 2. Models The principles of the GALPROP code for CR propagation and y-ray emission have been described in Strong & Moskalenko (1998) and Strong et al. (2000). Since then the codel has been entirely re-written in C++ (Moskalenko et al. 2002, and references therein) using the experience gained from the original (fortran) version, with improvements in particular in the generation of y-ray skymaps. Both 2D (radially symmetric) and full 3D options are available, the latter allowing also explicit time-dependence with stochastic SNR source events (Strong & Moskalenko 2001b). For this paper the 2D option is sufficient since we need only kpc-scale averaged CR spectra (even if these differ from local CR measurements). An important point to note is that even in the 2D case, the symmetry applies only to the CR distribution; for the gas-related components (.lro-decay and bremsstrahlung) of the -+ y-ray skymaps we use 21-cm line survey data for H I and CO (J = 1 0) survey data for Hz, in the form of column densities for Galactocentric "rings," using velocity information and a rotation curve (see Appendix A for details). In this way details of Galactic structure are included at least for the gas, at a sufficient level for the present limited state of knowledge on e.g. the relation of cosmic rays to spiral structure. The longitude range 350" < I < 10" is not included in the H I and CO survey data due to lack of kinematic information; for the analysis interpolated values are used, and this is found to be fully consistent will the y-ray data. The interstellar radiation field (ISRF) for computing IC emission and electron energy losses is the same as that described and used in Strong et al. (2000); pending a new calculation (an ambitious project) this is the best we havc available. Although the uncertainty in the ISRF is a shortcoming, note that since we fit to the y-ray data by adjusting the electron spectrum, As usual the code and documentation is available at http://www.mpe.mpg.de/-aws/aws.html - 6 - inaccuracies in the ISRF spectrum will tend to be compensated. The radial distribution of CR sources used is the same as in Strong et al. (2000), since we find this empirically-derived form still gives a good reproduction of the y-ray longitude distribution. Although flatter than the SNR distribution (e.g., Case & Bhattacharya 1998), this may be compensated by the gradient in the CO-to-H, conversion factor whose metallicity and temperature dependences have the net effect of causing the factor to increase with R (Papadopoulos, Thi, & Viti 2002; Israel 1997). We use a uniform value of XCO= 1.9~10~’ molecules cm-2/(K km s-l) as in Strong et al. (2000) and Strong & Mattox (1996); this is consistent with the value (1.8 f O.3)x1O2O molecules cm-2/(K km s-l) from a recent (non-y-ray) CO survey analysis by Dame et al. (2001). The parameters of the models are summarized in Table 1. The models differ only in the injection spectra of protons (and He) and electrons, while the injection spectra of heavier nuclei are assumed to have the same power-law in rigidity, for all models. For propagation, we use essentially the same diffusion reacceleration model, model DR, as described in Mos- kalenko et al. (2002). The propagation parameters have been tuned to fit the B/C ratio (Fig. 1) using improved cross-sections (Moskalenko & Mashnik 2003). The spatial diffu- sion coefficient is taken as pDo(p/p0)’, where Do = 5.8 x cm s-’ at po = 4 GV, and 6 = 1/3 (Kolmogorov spectrum). The Alfvh speed is WA = 30 km s-’. The halo height is taken as zh = 4 kpc as in Strong et al. (2000), in accordance with our analysis of CR sec-, ondary/primary ratios (Moskalenko, Mashnik, & Strong 2001; Moskalenko et al. 2002, and references therein). However, values zh differing by 50% (the estimated error) with corre- sponding adjustment of Do would provide essentially similar results since the IC contribution scales mainly with the electron spectrum which is here treated as a free parameter. The spectra are compared in the regions summarized in Table 2. Region A corresponds to the “inner radian”, region B is the Galactic plane excluding the inner radian, region C is the “outer Galaxy”, regions D and E cover higher latitudes at all longitudes, region F is “Galactic poles”. Region H is the same as in Hunter et al. (1997) and is used for comparison with results of Hunter et al. In addition to spectra, profiles in longitude and latitude are an essential diagnostic; our latitude profiles are plotted logarithmically because of the large dynamic range from the Galactic plane to the poles. The EGRB used here is based on the new determination by Strong et al. (2004). Since this was derived for the EGRET energy bands, it is interpolated in order to produce a contin- uous spectrum for combining with the model Galactic components. The present analysis is however not sensitive to the details of the EGRB. Since our COMPTEL data do not contain the EGRB (see Section 3.2), we do not extrapolate the EGRB beyond the EGRET energy range when comparing with data. - 7 - The output of the GALPROP runs is in the form of FITS files; the visualization2 in the form of spectra and profiles, and comparison of the results with data involves integrations over sky regions and energy as well as convolution. The predicted model skymaps are convolved with the EGRET point-spread function as described in Strong et al. (2000). For the profiles the convolved model is directly compared with the observed intensities. For the spectra, the procedure is slightly different: the predicted (unconvolved) intensities are compared with intensities corrected for the effect of convolution as given by the model under study. This procedure has the advantage that the spectra are spatially deconvolved, allowing for more direct interpretation and also the combination of data with other experiments, such as COMPTEL, with different instrument response functions. The effect of this procedure on the spectra is only significant below 500 MeV. 3. ?-ray and cosmic ray measurements 3.1. EGRET data We use the co-added and point-source removed EGRET counts and exposure maps in Galactic coordinates with 0.5" x 0.5" binsize at energies between 30 MeV and 10 GeV, as described in Strong et al. (2000). Apart from the most intense sources, the removal of sources has little influence on the comparison with models. For the spectra, the statistical errors on the EGRET data points are very small since the regions chosen have large solid angle; the systematic error dominates and we have conservatively adopted a range 515% in plotting the observed spectra (Sreekumar et al. 1998; Esposito et al. 1999). For longitude and latitude profiles only the statistical errors are plotted. In addition we use EGRET data in the energy ranges 10-20, 20-50 and 50-120 GeV. Because the instrumental response of EGRET determined at energies above 10 GeV is less certain compared to energies belaw 10 GeV, it is necessary to account for additional uncertainties. In particular the EGRET effective area can only be deduced by extrapolation from the calibrated effective area at lower energies (Thompson et al. 1993a). We accordingly adopt values of 0.9, 0.8, and 0.7 times the 4-10 GeV effective area, respectively. On top of the statistical and systematic uncertainties as described above we account for the uncertainties due to the uncalibrated effective area of the EGRET telescope above 10 GeV with an additional systematic error of f5%. However, the actual number of photons >10 GeV is small: 1091, 362 and 53 events respectively, and concentrated mainly in the inner Galaxy; hence the comparison with models above 10 GeV *An additional program GALPLOT has been developed for this purpose, with flexible plotting options and convolution; this will be made available with future versions of GALPROP. b -8- can only be made in this region. At low energies the EGRET effective area includes the so-called “Kniffen factor” (Thomp- son et ai. 1993b) derived by fitting the Crab spectrum; this additional uncertainty (factor = 2 3.4 for 30-50 MeV and 1.2 - 1.6 for 50-70 MeV) should be borne in mind when - comparing models with EGRET data. 3.2. COMPTEL data The intensities are based on COMPTEL maximum entropy all-sky maps in the energy ranges 1-3, 3-10 and 10-30 MeV, as published in Strong et al. (1999). The intensities are averaged over the appropriate sky regions, with high latitudes being used to define the zero level. COMPTEL data is only used for the inner Galaxy spectra, since the skymaps do not show significant diffuse emission elsewhere. For this reason, the COMPTEL data shown in the figures does not include the EGRB. 3.3. Cosmic rays In our conventional model we use the locally-observed proton, He, and electron spec- tra (solid lines in Figs. 2, 3). The nucleon data are now more precise than those which were available for Strong et al. (2000). The proton (and Helium) injection spectra and the propagation parameters are chosen to reproduce the most recent measurements of primary and secondary nuclei, as described in detail in Moskalenko et al. (2002). The error on the - dominant proton spectrum in the critical (for TO-decay) 10 - 100 GeV range is now only 5% for BESS (Sanuki et al. 2000). Relative to protons, the contribution of He in CR to the y-ray flux is about 17%, and the CNO nuclei in CR contribute about 3%. The He nuclei in the ISM contribute about 25% relative to hydrogen for the given ratio He/H = 0.11 by number. The total contribution of nuclei with 2 > 1 is about 50% relative to protons. In our optimized model we use the proton and He high-energy spectral shape derived from the local data (dotted lines in Figs. 2, 3). We allow however for some deviations in the normalization. The antiproton (Orito et al. 2000; Beach et al. 2001) and positron data provide an important constraint (Moskalenko et al. 1998; Strong et al. 2000) on the proton spectrum on a large scale. Since the low-energy protons are undetectable in the ISM, we allow more freedom in the proton and He spectrum below 10 GeV. We introduce a break at 10 GeV which enables us to fit the y-ray spectrum while still remaining within the constraints provided by the locally-observed antiproton and positron spectra. The deviations from the . - 9 - local measurements at low energies can be caused by the effect of energy losses and spatial fluctuations in the Galaxy. The modification of the low-energy proton spectrum may also be partly a compensation for errors in the models of neutral pion production at low energies (e.g., Stecker 1970), which rely on the data of 1960’s (see Dermer 1986, and references therein) and do not provide the required accuracy now. Besides, the low-energy protons are strongly affected by solar modulation; while the effect of solar modulation is not fully understood, it is essential below 10 GeV. We refer a reader to Section 8 where various aspects of the uncertainties are discussed in more detail. For electrons, the injection index near -1.8 at -1 GeV is consistent (see Strong et al. 2000) with observations of the synchrotron index /I = 2.40 - 2.55 for 22-408 MHz (Roger et al. 1999) and /3 = 2.57 f 0.03 for 10-100 MHz (Webber, Simpson, & Cane 1980). Secondary and tertiary antiprotons are calculated as described in Moskalenko et al. (2002). Secondary positron and electron production is computed using the formalism de- scribed in Moskalenko & Strong (1998), which includes a reevaluation of the secondary &- and K*-meson decay calculations. Antiprotons, positrons, and electrons including secondary electrons are propagated in the same model as other CR species. 4. Conventional model We start by repeating the test of the “conventional” model; the y-ray spectra in the 7 test regions are shown in Fig. 4. As required by the ‘[conventional” tag, the proton and electron spectra are consistent with the locally observed spectra (Figs. 2,3). This is the same “conventional” model as in Strong et al. (2000), with updated nucleon spectra, but because we compare with a more complete set of EGRET data than in Strong et al. (2000), the discrepancies become more explicit, and we can check whether they arise only in particular sky regions. Note that IC plays only a minor role in this type of model. As found in previous work, the GeV energy range shows an excess relative to that predicted; what is now evident is that this excess appears in all Zatitudes/Zongitude ranges. This is consistent with the results of Hunter et al. (1997) and Digel et al. (2001). It already shows that the GeV excess is not a feature restricted to the Galactic ridge or the gas-related emission. Further it is clear that a simple upward rescaling of the no-decay component will not improve the fit in any region, since the observed peak is at higher energies than the TO-decay peak. In other words, since the spectrum is very different from 7-ro-decay even at intermediate latitudes, a substantial IC component is required. Note that this version of the [Lconventiona17s’p ectrum is nevertheless in rather better agreement with EGRET data than in Strong et al. (2000), due to inclusion of secondary - 10 - positrons/electrons, general improvements in the model (e.g., iTTO-decay, improved gas data) and the EGRET data treatment (direct use of the count and exposure data instead of the model-fitting analysis of Strong & Mattox 1996). The improvement is especially evident in the 30-100 MeV range, where secondary positrons/electrons make a substantial contribution (see Section 7). A test against antiproton and positron data also shows “excesses”. The conventional model with reacceleration is known (Moskalenko et al. 2002) to produce a factor of -1.5 (~2.50)le ss antiprotons at 2 GeV than measured by BESS (Orito et al. 2000). The antipro- ton spectrum for the conventional model is shown in Fig. 5. Positron data, though scattered, also show some “excess” at high energies (Fig. 6). It is thus clear that the “excesses” in GeV y-rays in all directions, in GeV antiprotons, and in positrons above several GeV found in the conventional model indicate that the average high-energy proton flux in the Galaxy should be more intense or our reacceleration model is invalid or there is a contribution from uncon- ventional sources (e.g. dark matter). For more discussion of antiproton and positron tests see Section 7. In the “SNR source” scenario of Berezhko & Volk (2000) the y-ray spectrum in the inner Galaxy is attributed to an additional population of unresolved SNR, but this component cannot explain the excess observed at high latitudes, and hardly in the outer Galaxy3. The presence of the GeV excess in all sky regions is also a problem for the suggestion by Aharonian & Atoyan (2000) of a hard proton spectrum in the inner Galaxy. These explanations are therefore by themselves insufficient, although they could give a contribution. 5. Hard electron injection spectrum model This model is essentially as in Strong et al. (2000), recomputed with the current GAL- PROP code. The main feature is the electron injection index of 1.9. Comparison of the spectra in the 7 sky areas (Fig. 7) show that this model reproduces the GeV excess except in the inner Galaxy (region A) where it is still too low. However the spectral shape is not well reproduced. More significant, comparing with the new EGRET data above 10 GeV in the inner Galaxy, the spectrum is much too hard. Fig. 3 compares the locally observed electron spectrum with that from the model; the deviation at high energies is much larger than expected from the 3D study by Strong & Moskalenko (2001b). As discussed in the Introduction, there are therefore a number of reasons to lead us to consider this model as after all untenable. 3Berezhko & Volk (2000) did not address the question of regions away from the inner Galaxy.