Diffuse Galactic gamma–ray flux at very high energy
Abstract
The observation of the diffuse Galactic gamma ray flux is the most powerful tool to study cosmic rays in different regions of the Galaxy, because the energy and angular distributions of the photons encode information about the density and spectral shape of relativistic particles in the entire Milky Way. An open problem of fundamental importance is whether cosmic rays in distant regions of the Milky Way have the same spectral shape observed at the Earth or not. If the spectral shape of protons and nuclei is equal in all the Galaxy, the dominant, hadronic component of the diffuse gamma ray flux must have an angular distribution that, after correcting for absorption effects, is energy independent. To study experimentally the validity of this factorization of the energy and angular dependence of the diffuse flux it is necessary to compare observations in a very broad energy range. The extension of the observations to energies –10 PeV is of great interest, because it allows the study of the cosmic ray spectra around the feature known as the “knee”. The absorption probability for photons in this energy range is not negligible, and distorts the energy and angular distributions of the diffuse flux, therefore a precise calculation of the absorption effects is necessary for the interpretation of the data. In this work we present predictions of the diffuse gamma ray flux at very high energy, constructed under different hypothesis for the space dependence of the cosmic ray energy spectra, and discuss the potential of the observations for present and future detectors.
pacs:
98.35Gi,95.85Pw,95.85RyI Introduction
Direct observations of cosmic rays (CR) near the Earth only measure the fluxes present in a small region of the Milky Way in the vicinity of the solar system, however it is clearly of fundamental importance to study the spectra also in other regions of the Galaxy. The most powerful method to obtain information on the spectra of relativistic particles in distant regions of the Galaxy is the study of the diffuse fluxes of ’s and ’s. Cosmic ray propagating in interstellar space can interact with gas or radiation fields, and these inelastic interactions generate gamma rays and neutrinos. These secondary particles travel along straight lines and are observable at the Earth, forming fluxes that encode the space and energy distributions of CR in the entire volume of the Galaxy.
An open question of great importance in high energy astrophysics is whether the CR spectra have the same shape in all points of the Galaxy, or there is a non trivial space dependence. To study this problem experimentally it is clearly very desirable to measure the diffuse gamma ray and neutrino fluxes in a very broad energy range.
The total flux of rays can be naturally decomposed into the sum of several components: (i) an ensemble of point–like or quasi–point like sources; (ii) an isotropic flux of extragalactic origin; (iii) a diffuse Galactic flux generated by the interactions of cosmic rays in interstellar space.
In the energy range 0.1–1000 GeV the diffuse flux is the largest of the three components, and has been detected and studied by several gamma ray telescopes on satellites. Measurements of the diffuse flux have been obtained by OSO–3 oso3 (), SAS–2 sas2 (), COS–B cosb () and EGRET egret (). More recently the Fermi telescope has obtained accurate measurements of the flux over the entire sky Ackermann:2012pya (); Acero:2016qlg ().
The diffuse flux is concentrated in a narrow region around the Galactic equator, with one half of the total in the latitude range . The flux is also larger toward the Galactic Center, and the contribution from the longitude region is approximately twice as large as the flux from the opposite hemisphere.
Ground based gamma ray detectors Atkins:2005wu (); Abdo:2008if (); Bartoli:2015era (); Abramowski:2014vox () have also obtained measurements of the diffuse ray flux in the TeV energy range for some regions near the Galactic equator. At higher energy ( TeV) there are at present only upper limits for the diffuse flux Chantell:1997gs (); Borione:1997fy (); Apel:2017ocm (), however new detectors (such as IceTop and LHAASO) have the sensitivity to extend the observations to the PeV energy range. These measurements of the diffuse Galactic gamma ray flux at very high energy can be very important in the study of a possible space dependence of the CR spectra.
Recently the IceCube detector Aartsen:2013jdh (); Aartsen:2014gkd (); Aartsen:2015rwa (); Aartsen:2017mau (). has obtained evidence for the existence of an astrophysical signal of very high energy neutrinos ( TeV). The simplest interpretation of the data is that most of the signal is of extragalactic origin, however it is also possible that there is a subdominant contribution due to Galactic emission. Some authors have also speculated that most, or even all of the IceCube signal is of Galactic origin. This however requires a non–standard emission mechanism, because the angular distribution of the neutrino events that form the signal, in contrast to gamma ray observations at lower energy, is approximately isotropic.
If the neutrinos are generated by the standard mechanism of production, the emission is accompanied by a emission of similar intensity and spectral shape. In the energy range TeV, photons traveling over extragalactic distances are completely absorbed by interactions with low energy radiation fields. On the other hand the flux of Galactic gamma rays is suppressed and distorted by absorption effects but remains observable. The conclusion is that the Galactic component of the IceCube signal should have an observable gamma ray counterpart. The comparison of the gamma ray and neutrino fluxes in the same energy range is therefore very important to clarify the origin of the IceCube signal. In this study it is important to take into account the effects of absorption of the flux of high energy photons, that depend on the space distribution of the emission and have non trivial energy and angular dependences.
The goal of this paper is to discuss the potential of existing and future observations of the Galactic diffuse flux of gamma rays at very high energy. In our discussion we will construct and compare two predictions of the diffuse flux at very high energy based on two alternative frameworks to extrapolate the measurements performed at lower energy. In both cases the dominant source of Galactic diffuse emission is the so called hadronic mechanism, where the gamma rays are generated in the inelastic collisions of protons and other nuclei. In one model we will assume that CR protons and nuclei have the same spectral shape, equal to what is observed at the Earth, in all of the Galaxy. This implies that the energy distribution of the gamma ray emission has a shape that is independent from position and therefore, if photon absorption during propagation is negligible, that the angular distribution of the flux at the Earth has an energy independent shape. In an alternative model, following some previous studies Gaggero:2014xla (); Acero:2016qlg (); Yang:2016jda (), we will assume that the CR in the inner part of the Galaxy have a harder spectrum than what is observed at the Earth. Accordingly the space and energy dependences of the emission are not factorized and the angular shape of the diffuse flux is energy dependent, with the fraction of the diffuse flux that arrives from the region around the Galactic Center increasing with energy.
Some “non–standard” models for the diffuse Galactic gamma ray flux at very high energy, recently proposed on the basis of the IceCube results, and that predict a very different angular distribution, will be also very briefly discussed.
The work is organized as follows. In the next section we review some general properties of the diffuse Galactic flux of gamma rays. In section III we compute the “local” rate of emission of gamma rays in the vicinity of the solar system. This calculation requires only a knowledge of the CR spectra that are directly observable at the Earth. In section IV, we construct a simple (cylindrically and up–down symmetric) model for the interstellar gas density that is adequate to describe the main large scale features of the distribution. Section V contains a brief discussion of the existing observations, in particular those performed by the Fermi telescope. The following two sections present two models for the extrapolations to high energy of the Fermi measurements. Section VIII discusses some “non–standard” model for the Galactic diffuse flux motivated by the IceCube results. The following section very briefly discusses the problem of the separation of the diffuse flux from the flux generated by the ensemble of all discrete Galactic sources. In section X we review the existing measurements of the diffuse Galactic flux at high energy and compare the results to our calculations. The section also discusses the perspectives and necessary conditions to extend the measurements of the diffuse flux up to the PeV energy range. A final section gives a summary and some conclusions.
Ii The diffuse Galactic ray flux
The diffuse flux of gamma rays with energy from the direction can be calculated integrating the emission along the line of sight and including a correction for absorption effects:
(1) 
In this expression is the emission rate, that is number of gamma rays of energy emitted per unit volume, unit time and unit energy from the point . The integration is over all points along the line of sight, with the position of the solar system, and the versor in the direction . The factor 1/() follows from the assumption that the emission is isotropic. The exponential factor in Eq. (1) gives the survival probability of gamma rays during propagation, and is the optical depth.
The dominant mechanism for gamma ray emission at high energy is the so called “hadronic mechanism”, that is the creation and decay of unstable mesons (mostly with smaller contribution of other particles such as and ) in the inelastic interactions of protons and other CR nuclei with interstellar gas. The largest contribution to the hadronic emission is due to interactions between relativistic protons and the hydrogen component of the interstellar gas. This contribution can be calculated as:
(2) 
where is the number density of hydrogen gas at the point , is the flux of CR protons with energy at the same point, is the inelastic cross section and is the inclusive spectrum of gamma rays generated in a interaction after the decay of all unstable particles created in the collision. The integration is over all proton energies that can generate photons with energy . Interactions where the projectile and/or the target is a nucleus (such as –helium, helium–, helium–helium, and so on) also contribute to the hadronic emission and can be calculated with expressions that have the same structure as Eq. (2) with obvious substitutions.
Smaller contributions to the gamma ray emission are generated by leptonic processes where the gamma rays are radiated by CR electrons and positrons, via bremsstrahlung and Compton scattering. For bremsstrahlung (interactions such as ) the target, as in the hadronic emission case, is interstellar gas. For Compton scattering () the target is the ensemble of the soft photons that form the radiation fields in space. In this case it is necessary to model not only the density, but also the energy spectrum and angular distribution of the target particles.
ii.1 Gamma ray absorption
The most important process that can absorb photons during propagation in interstellar space is pair production interactions () where high energy gamma rays interact with the soft photons that form the Galaxy radiation fields. The interaction probability per unit length ) for a photon of energy and direction at the space point can be calculated, integrating over the energy and angular distributions of the target photons:
(3) 
In this expression is the 3–momentum and is the energy of the target photon, is the cosine of the angle between the interacting particles, and is the pair production cross section, that can be expressed as a function of the square of the center of mass energy .
The optical depth for photons observed at the Earth with energy , direction that have traveled a distance , can be calculated integrating the interaction probability along the photon trajectory:
(4) 
The calculation of the absorption probability and of the optical depth requires a sufficiently accurate knowledge of the energy and angular distribution of the target photons. An extended discussion of this problem is contained in Vernetto:2016alq () (see also moska2006 () and references therein). Some of the main properties of high energy gamma ray absorption are illustrated in Fig. 1. The top panel of the figure shows the (angle integrated) energy distributions of the target photons in the vicinity of the solar system according to Vernetto:2016alq (). The distribution is the superposition of three main components: the cosmic microwave background radiation (CMBR), infrared emission, and stellar light. A minor contribution is given by the Extragalactic Background Light (EBL).
The CMBR fills homogeneously all space with an isotropic blackbody spectrum of temperature Kelvin; this corresponds to a total number density cm of photons with average energy eV.
Infrared photons are radiated by interstellar dust heated by stellar light. This emission can be reasonably well described as a diluted, and distorted black body spectrum [] with a temperature of approximately 20 Kelvin and a distortion parameter of order 1.5–1.7. At high energy ( eV) the spectral shape deviates from this form because of the contribution of an ensemble of emission lines radiated by the smallest dust grains that are not in thermal equilibrium. The infrared radiation has an average energy of order 0.008 eV, and a number density cm.
Stellar light can be described as the superposition of diluted black body spectra with temperatures between 3000 and 8000 Kelvin, plus a small contribution in the ultraviolet range from young hot stars. In the vicinity of the solar system the stellar light radiation field has a total number density of order 0.5 cm of photons with average energy eV. The infrared and stellar light components of the radiation field have non trivial space and angular distributions that reflect the disk structure of the Galactic sources.
The bottom panel of Fig. 1 shows the angle averaged absorption probability in the solar neighborhood. The energy dependence of this absorption probability reflects the spectral shape of the target photon distribution. The maximum at PeV, and the two shoulders at 150 and 1.6 TeV correspond to interactions with the photons of the three main components (CMBR, dust and star emission) of the target radiation field. The probability of interactions with the photons of a single component has a maximum for a gamma ray energy of order , that corresponds to the c.m. energy of the photon–photon collisions just above the kinematical threshold, where the pair production cross section has its maximum value (where cm is the Thomson cross section). At the maximum, the absorption probability takes the value where is the total number density of target photons that form the component.
Numerically this corresponds to a minimum interaction lengths (in the solar neighborhood) of order kpc at energy PeV for absorption by the CMBR, kpc at TeV for absorption by the infrared dust emission, and Mpc at TeV for absorption by starlight.
The calculation of the optical depth requires a knowledge of the target radiation field in the entire volume of the Galaxy, however, for a qualitative understanding, one can note that the spectra of the target photons have a similar shape in all points of the Galaxy. The absorption generated by interactions with the CMBR, with an absorption length of order 10 kpc, that is of same order of the linear size of the Galaxy is very important for the propagation of photons in the PeV energy range. The effects of absorption by dust emitted infrared photons, with a (space and direction dependent) absorption length ten times longer (of order 100 kpc), are smaller but not entirely negligible. The effects of absorption on stellar light remain always small.
Iii The local diffuse ray emission
As a first step, in this section we will calculate the local diffuse gamma ray emission, that is the emission in the vicinity of the solar system. This calculation requires three elements: (a) a knowledge of the CR fluxes that are directly observable at the Earth, (b) a description of the relevant targets (gas and radiation) for CR interactions, and (c) a model for the interaction cross sections. The crucial point is that the calculation does not need to model the space dependence of the CR spectra.
In the following we will discuss separately the two main (hadronic and leptonic) emission mechanisms.
iii.1 Hadronic emission
The calculation of the hadronic emission requires a description of the nuclear components of the CR flux. Fig. 2 shows our fit to the observed spectra of protons and nuclei, together with some of the available data. In the figure the spectra are shown in the form of the nucleon flux versus the energy per nucleon . The spectra exhibit two evident spectral features. The first one is a hardening at rigidity GV that has been observed by CREAM, PAMELA and AMS02 Ahn:2010gv (); pamelahardening (); ams02protons (); ams02helium (). The second feature is the well known “knee” that is prominent in the all–particle spectrum observed by air shower experiments for a particle energy PeV. Our description of the CR fluxes assumes that also the “knee” is a spectral structure present for all nuclear species at a constant rigidity, with the softening at 4 PeV corresponding to the break in the helium component that is dominant at that energy.
For GeV our model of the CR spectra is based on a generalization of the fit given in Lipari:2017jou (), that gives a good description of the data of AMS02 ams02protons (); ams02helium () and CREAM Yoon:2017qjx (). The fluxes of nuclei with are extrapolated from the measurement of the HEAO3 detector Engelmann:1990zz (), introducing a rigidity dependent hardening.
The CR fluxes observed at the Earth, for GeV are distorted by solar modulation effects. Our model of spectra in interstellar space are demodulated using the Force Field Approximation with a parameter GeV. In the present work we discuss only high energy gamma rays, and the description of the solar modulation effects is of negligible importance.
At high energy, our model of the CR spectra smoothly joins the parametrization of Gaisser, Stanev and Tilav (GST) Gaisser:2013bla (), developed to fit the measurements of Extensive Air Shower detectors, and the two models become identical for GeV.
The target of the hadronic interactions is the interstellar gas, and is entirely characterized by its density and chemical composition. We have assumed the average solar system composition estimated by Ferrière Ferriere:2001rg (), with hydrogen, helium and all other nuclei accounting for fractions 0.90, 0.0875 and 0.0125 of the atoms. This corresponds to a total mass of the interstellar gas that is a factor 1.42 larger than the mass in hydrogen. The numerical results shown below are calculated for a nominal value of the hydrogen density cm, to allow a simple rescaling. In the following (see section IV) we will estimate that the average hydrogen density on the Galactic plane at radius kpc is cm (note that this average quantity is not identical to the gas density in the vicinity of the solar system that is determined by local effects).
The calculation of the hadronic emission requires a model for particle production in inelastic hadron–hadron collisions. Because of the present limitations in the understanding of the strong interaction processes, a non negligible source of uncertainty cannot be avoided. The interaction probabilities of relativistic hadrons are well understood, because the total, elastic and inelastic cross section have been measured with good precision in the entire energy range of interest. The cross sections for –nucleus and nucleus–nucleus collisions have been estimated from the data on collisions using a standard Glauber calculation Glauber1970 (). Collisions with helium and other nuclei (nuclei with have been modeled as entirely formed by oxygen) account for a weakly energy dependent fraction of order 21–23% of all inelastic interactions.
For the description of the inclusive spectra of final state particles, in the lower energy range ( GeV) we have used the algorithms described in Lipari:2016vqk (). At higher energy we have modeled the interactions using the Pythia Montecarlo code Sjostrand:2006za ().
The local gamma ray emission from hadronic interactions is shown in Fig. 4. The spectrum exhibits a strong suppression for GeV, that has a simple and well known kinematical origin, associated to the fact that most of the gamma rays are generated in the decay of neutral pions. It is straightforward to show that if the photons are entirely generated by the decay of neutral pions, the spectrum has the symmetry: (with the mass). This property implies that the spectrum has a maximum at , and that the emission of photons with energy is strongly suppressed. The observation of this spectral feature can in principle be used to identify in a model independent way the hadronic component in the gamma ray flux.
At higher energy () the hadronic gamma ray emission falls roughly as a power law, but the spectral index is not constant and changes gradually with energy. Inspecting Fig. 4 one can note a gradual hardening of the spectrum centered at GeV, and then a more marked, but also gradual softening centered at GeV. These broad structures reflect the features present in the primary CR spectra that we have discussed above.
For a qualitative understanding one can note that for GeV, the inclusive spectra of the final state particles, created in inelastic hadronic collisions, have (in reasonably good approximation) the scaling property:
(5) 
This equation can be derived assuming the validity of Feynman scaling in the projectile fragmentation region, and has the well known consequence that if the spectrum of primary nucleons is a simple power law with exponent , the emission is then a power law with the same spectral index. The violations of Feynman scaling, and the growth of the inelastic cross sections with c.m. energy introduce logarithmic corrections.
The all nucleon flux cannot however be described as a simple, unbroken power law because of the existence of the spectral features that we have discussed above: the “Pamela hardening” and the “knee”. The existence of these structures in the primary CR spectra are reflected in more gradual features in the emission, that are centered at (a factor of order 5–10) lower energy. This can be easily understood noting that photons with energy are generated by the interactions of primary nucleons with in a broad range of energy, with an order of magnitude extension and a median value (the precise value depends on the spectral index).
iii.2 Leptonic emission
The calculation of the leptonic emission requires a description of the flux of electrons plus positrons. Our fit to the ( is shown in Fig. 3 together with some of the measurements. The flux is accurately measured for GeV by the observations of detectors on satellites like PAMELA pamelaelectrons (); pamelapositrons (); Adriani:2014pza () Fermi Abdollahi:2017nat () and AMS02 amsallelectrons (). The observations of HESS Aharonian:2008aa (); Aharonian:2009ah (); hessicrc2017 (), and later by MAGIC BorlaTridon:2011dk (), VERITAS Staszak:2015kza (), and more recently by DAMPE Ambrosi:2017wek () have shown that the spectrum has a break at GeV where the spectrum steepens from a spectral index of order 3.1 to an index of order 3.8.
Two mechanisms contribute to the leptonic emission. In bremsstrahlung the target (interstellar gas) is identical to the one discussed for the hadronic emission. For Compton scattering the target are the photons of the interstellar radiation fields discussed in sec. II.1. For our calculations we have used the model of Vernetto:2016alq ().
The leptonic mechanisms for gamma ray production are purely electromagnetic and therefore have exactly calculable cross sections (see for example Blumenthal:1970gc ()).
The leptonic emission, separated into the contributions of bremsstrahlung and Compton scattering, is shown in Fig. 4. The results can be easily understood qualitatively. In the case of bremsstrahlung, the radiates photons with an energy independent cross section and the final state photon has an energy distribution that depends only on the ratio . If the primary have a power law spectrum, the emission is then also a power law with the same spectral index (of order for GeV). The bremsstrahlung spectrum softens at higher energy because of the break in the ( spectrum at TeV.
The Compton scattering component of the emission has initially a hard spectrum (a spectral index of order 2). This reflects the well known fact that when the interactions are in the Thomson regime (that is when the product of the energies of the interacting particles is sufficiently small: ) the spectral indices of the Compton emission and the primary electron flux are related by: . This behaviour however stops for GeV when most of the interactions are in the Klein–Nishina regime. The Compton emission suffers more suppression at higher energy also because of the softening of the spectra above 1 TeV. The result is that the local Compton emission of gamma rays gives a maximum contribution of order 5% with respect to the hadronic one.
It should be noted that the estimate of the contribution of the leptonic mechanisms to the observed gamma ray flux (that is the result of the emission from the whole Galaxy) is a more difficult task, because it requires to compute the emissions in different regions of the Galaxy, where the densities of the primary particles (electrons, protons and nuclei) and of the relevant targets (gas and radiation) can be different. In particular it is possible, and indeed likely that the Compton mechanism can be a significant component of the flux for directions that point away from the Galactic equator. This is because the interstellar gas density (the target for hadronic emission) is exponentially suppressed for large , while the density of the radiation fields (the target for Compton scattering) falls more gradually with (and the CMBR component is in fact constant). However one expects that in the region of small latitude (), where the diffuse flux is largest, the leptonic mechanisms remains subdominant, with a maximum contribution of order to the diffuse flux.
iii.3 Summary
The main results obtained in this section can be summarized as follows:

The hadronic mechanism is the dominant source of gamma rays for few GeV in the Galactic disk region () with the leptonic mechanisms accounting for only a few percent of the emitted photons.

In a broad energy interval between 10 and GeV the emission spectrum is in good approximation a power law () with a spectral index that varies slowly with values between 2.8 and 2.6.

For between and GeV the energy distribution of the emission softens, and the spectral index grows gradually up a value of order 3.1. Between 1 and 10 PeV the spectral index remains approximately constant.
Iv Galactic interstellar gas density
To model the diffuse gamma ray emission in the whole Galaxy, a knowledge of the spatial distribution of the density and composition of the interstellar medium is required. The interstellar matter has been the object of many studies, and several large scale surveys of its main components have been performed in recent years. These studies have revealed that the space distributions of interstellar gas and dust have a complex form, with structures present at many scales.
Our goal in the present work is to construct a model of the interstellar gas density that captures reasonably well its large scale properties, without attempting to describes accurately finer details such as individual clouds and filaments. For this purpose we have assumed an axially and updown symmetric distribution of the gas that neglects the spiral arms, whose geometry remain controversial, and northsouth asymmetries such as the disk warp.
The model is constructed using previous studies on the distribution of neutral atomic (H), molecular (H) and ionized hydrogen, that are the most important components of the interstellar medium. The contribution of other elements is added to the hydrogen assuming that the composition of the interstellar gas is equal in the entire Galaxy, and equal to the composition of the solar system estimated by Ferrière Ferriere:2001rg () (see section III.1).
The interstellar gas is confined to a narrow region around the Galactic plane with a thickness that grows with (the so called “disk flaring”). Assuming for the dependence a Gaussian form, the hydrogen density can be written as:
(6) 
with and cylindrical coordinates. In this equation is the midplane density, and is related to the half width at half maximum (HWHM) of the distribution by the relation . An important quantity is the surface density
(7) 
that is related to the midplane density an the disk thickness by the relations:
(8) 
iv.1 Molecular hydrogen
Molecular hydrogen cannot be observed directly, but its spatial distribution can be inferred through the observation of the emission lines of carbon monoxide (CO).
To describe the distribution in the Galactic disk, excluding the central region, we use the results of RomanDuval et al. RomanDuval:2016 (), that report the radial distributions of the surface density and the disk thickness for from 2.5 to 16 kpc, obtained analyzing CO and CO data. Their evaluations are consistent with the measurements reported in the review by Heyer and Dame Heyer:2015 (). According to their results, the H surface density has a maximum at kpc. Smoothing the small scale granularity of the observed distribution, the radial dependence of the surface density can be approximately described by exponential functions, with a slope change at kpc. According to the data reported in the same paper, the disk thickness for kpc is approximately constant, pc, and grows exponentially at larger radii.
Concerning the central region of the Galaxy, we use the work of Ferriere et al. Ferriere:2007 (), who model the molecular gas distribution by combining different sets of data. According to their study, in the Galaxy Center the most prominent feature is a small region of radius 200 pc with a very high H density, the Central Molecular Zone (CMZ), actually centered at 50 pc from the Galactic Center, embedded in a lower density ring extending up to kpc. Both these structures have an elongated shape and are tilted with respect to the Galactic plane and to the line of sight. Since their precise geometry is not known and since we only need to describe the general shape of the gas distribution, for simplicity we model the central region assuming an axisymmetric behaviour, with the surface density exponentially decreasing with . To find the radial slope and normalization, we set the total H mass in the (CMZ + ring) region equal to the mass evaluated by Ferriere et al. (5.3 10 M) and the surface density at kpc equal to the value obtained by the function used for the disk, previously described.
Concerning the distribution along , according to Ferriere et al. the thicknesses of the CMZ and the ring are approximately costant, but the CMZ is thinner ( pc) than the ring ( pc). We set pc at the center, and assume a flaring in order to connect smoothly with the data at larger radii.
Our parametrizations of the surface density, midplane density and disk thickness for molecular hydrogen are listed in Table 1. These parametrizations have a discontinuous derivative for some values of . In our numerical calculations the discontinuities are smoothed with a length scale pc. The model corresponds to a total mass of molecular hydrogen of M.
[kpc]  [ pc]  [cm]  [pc] 

0  1.5  105  135  15 
1.5  4.25  0.598  0.229  50 
4.25  8  7.5  2.88  50 
8  1.0  0.383  50 
iv.2 Atomic hydrogen
The distribution of neutral atomic hydrogen is studied throught the 21cm radio line, emitted in the transition between the two hyperfine levels of the ground state. We model the atomic hydrogen distribution according to the measurements reported by Kalberla and Dedes Kalberla:2008uu () for the region outside the Galactic Bulge, who fit the midplane density for r 8 kpc with the exponential form:
(9) 
where kpc and cm is the density at the Sun radius.
The midplane density is approximately constant between 4 and 8 kpc, and decreases rapidly towards the center of the Galaxy.
The vertical scale of the atomic hydrogen distribution grows exponentially with . Kalberla and Dedes describe the radial dependence of as
(10) 
and fit the observations with kpc and kpc.
The density of atomic hydrogen in the central zone of our Galaxy is about one order of magnitude lower than the molecular component. Similarly to the molecular case, we derive the H distribution starting from the evaluations by Ferriere et al. We describe the surface density for kpc with an exponential function having the same slope obtained for the molecular gas. We fixed the normalization in order to have the total mass of atomic hydrogen as reported in the same paper (5.2 10 M). We connect the density at =1.5 kpc to the density at =4 kpc (provided by the Kalberla and Dedes fit) with a further exponential curve.
Concerning the gas thickness, we set pc at the Galaxy Center (the value given by Ferriere et al. for the CMZ) and assume a flaring in order to connect smootly with the curve of equation (10).
Table 2 summarizes our parametrizations for atomic hydrogen in the Galaxy. As in the case of molecular hydrogen, discontinuities in the derivative of the functions are smoothed in the numerical calculations. The corresponding total mass of atomic hydrogen is M.
[kpc]  [ pc]  [cm]  [pc] 

0  1.5  10.2  3.21  45 
1.5  4  0.058  0.011  150 
4  8 
7.07  1.05  150 
8  10.6  1.05  150 
iv.3 Ionized hydrogen
The density of ionized hydrogen has beeen modeled by Taylor:1993my (); Cordes:2002wz (). In most of the Galaxy this component of the interstellar gas is negligible, however in the vicinity of the Galactic Center it is comparable to the contribution of neutral atomic hydrogen.
Fig.5 shows the midplane density of atomic, molecular and total hydrogen as a function of , for the whole Galactic plane.
V Fermi observations of the diffuse Galactic gamma ray flux
As a starting point for a model of the gamma ray diffuse emission at TeV–PeV energies, we use the existing data in the GeV energy range.
The latest and most accurate measurements of the diffuse flux in the energy range 0.1–100 GeV have been obtained in the last few years by the Fermi telescope. The Fermi collaboration has published in 2012 a dedicated paper about the diffuse Galactic emission Ackermann:2012pya () , however a significant amount of data has been obtained after that, and the methods of analysis have also significantly improved. Some of these new results are discussed in Acero:2016qlg ().
The Fermi data are public, and several authors (for example Gaggero:2014xla (); Yang:2016jda ()) have performed independent studies of the Galactic diffuse flux. In the present work we will not perform an independent analysis of the Fermi data, to estimate the diffuse gamma ray flux. This is a very important but difficult task that is postponed to a future work.
The Fermi Collaboration has made available a template of the diffuse Galactic gamma ray flux to be used as a background model for the search of point sources fermibackground (). This background model gives tables of the angular distribution of the flux (in bins equispaced in Galactic latitude and longitude with a linear size 0.125) for a discrete set of 30 energies (equispaced in ) between 58.5 MeV and 513 GeV.
In this paper we will use the Fermi background model
as a first order approximation of the Galactic diffuse flux.
We have chosen the map at the energy GeV (more precisely 11.98 GeV)
as a template of the angular distribution of the real diffuse flux at the same energy.
This template will be used here as a “boundary condition” for
extrapolations to higher energies.
The reference energy has been chosen as a reasonable
optimum choice on the basis of the following considerations.
(i) The energy must be sufficiently high, so that the contributions of
the leptonic mechanisms to the gamma ray flux is small
(see discussion in sec. III).
(ii) The energy must be sufficiently low, so that the diffuse
flux is measured with good statistical accuracy.
A low value of is also desirable because it
allows to study the evolution of the diffuse flux in a broader energy range,
when constructing different models for the extrapolation to very high energy.
Vi Model 1: space independent CR spectra
In the most commonly accepted models for Galactic cosmic rays, the spectral shape of the nuclear components (protons and nuclei) are identical in the entire volume of the Milky Way. This result emerges in a large class of models where the following conditions are satisfied:

The populations of CR in the Galaxy are in a stationary state, with the sources that compensate the losses due to escape and other effects. The spectra are not significantly distorted by the contributions of near sources that are still active or have been active in the recent past (time variations of the CR spectra associated to the evolution of the Galaxy can exist for cosmological time scales).

The spectra generated by the CR sources in different regions of the Galaxy have, after an appropriate average in time, a space independent shape. This condition is immediately satisfied in models where a single class of astrophysical events (for example SN explosions or GRB’s) is the dominant CR source.

CR propagation is well described by diffusion with a diffusion coefficient that has the same rigidity dependence in all points of the Galaxy. This corresponds to the statement that the space and rigidity dependences of the diffusion coefficient are factorized: .

Escape from the Galaxy is the dominant mechanism for CR losses.

Energy losses during CR propagation are negligible.
It is straightforward to see that if these conditions are satisfied the CR spectral shapes are proportional to the ratio of the time averaged source spectrum and the rigidity dependence of the diffusion coefficient. For example, for ultrarelativistic protons . The absolute normalization of the CR fluxes is in general a function of position , and depends on the geometry of the Galactic confinement volume and on the space distribution of the sources. Several authors have published interpretations of the CR data based on these assumptions, estimating the spectral shape of the source and the rigidity dependence of the diffusion coefficient.
The hypothesis that the nucleon flux has a spectral shape that is independent from the position can be written in the form:
(11) 
where is the locally observed spectrum, and an adimensional function of the space coordinates that, by construction, satisfies the constraint . Eq. (11) implies that the emission of gamma rays (and neutrinos) generated by the hadronic mechanism has also a factorized form:
(12) 
Inserting this factorized form of the emission in the general expression for the diffuse gamma ray flux of Eq. (1), and assuming that the absorption effects are negligible (that is in the limit ) one finds that the energy and angular dependences of the ray flux are factorized:
(13) 
where is a direction dependent length, given by the integral:
(14) 
The inclusion of absorption effects introduces an energy dependent distortion of the angular distribution of the flux, and therefore breaks the factorization of Eq. (13). A more general expression for the gamma ray flux can be written in the form:
(15) 
where we have indicated with the gamma ray survival probability, averaged over all points along the line of sight.
The crucial point of this discussion is that if the CR spectra have a space independent spectral shape, then the diffuse gamma ray flux has dependences on the energy and angle that are factorized when absorption effects are negligible. The factorization is broken at high energy ( TeV) when the absorption probability becomes significant.
Starting from these assumptions, we will evaluate the absolute gamma ray flux , using the interstellar gas model discussed in section IV, and introducing a simple parametrization for the space dependence of the CR flux:
(16) 
In this expression the function is the hyperbolic secant. This function falls exponentially for large values of the argument, but its derivative vanishes at as it is expected for the CR density at and .
In the framework we are discussing in this section, if the density of interstellar gas is known, the calculation of the gamma ray flux is a straightforward exercise, with results that are entirely determined by the two parameters and associated to the space dependence of the CR fluxes.
One example of the angular distribution of the gamma ray flux for the reference energy calculated under the factorization hypothesis for the CR spectra, and using the model of interstellar gas discussed in section IV, is shown in Fig. 6 and compared to the Fermi template. The figure shows the longitude distribution of the flux after integration in the latitude range , for different values of . Since the gas target is confined in a narrow region in , the calculation is insensitive to the value of if the parameter is sufficienty large ( kpc).
Inspecting Fig. 6 one can see that the Fermi template for the diffuse gamma ray flux has a rich structure with significant variations for angular scales as small as one degree. These rapid variations of the flux are consistent with the hypothesis that the gamma ray emission is proportional to the density of a very clumpy interstellar gas. Our calculation cannot reproduce the detailed structure of the flux for small angular scales, however it can describe reasonably well the large scale structure of the flux.
The hypothesis that the CR density is constant in the Galactic disk (and equal to what is observed at the Earth) can be excluded because such a model, that corresponds to the limit , predicts a flux that it too small for directions toward the Galactic Center, and too large for directions in the opposite hemisphere. A finite , that corresponds to a space gradient of the CR flux with a larger density in the GC region, results in a better agreement of our model with the Fermi template. The value kpc gives approximately the correct ratio for the contributions of the two hemipheres toward the Galactic Center and Anticenter.
Other illustrations of the calculated gamma ray flux at the reference energy are given in Fig. 7 that shows the Galactic longitude distribution after integration over the latitude ranges , and (that is the entire sky). Fig. 8 shows the latitude distribution of the flux at the reference energy , after integration over all longitudes.
The comparison of the model with the Fermi template shows that the main features of the diffuse gamma ray flux can be described reasonably well. The largest discrepancies are observed at large , and the effect is likely associated to the existence to the structures known as the “Fermi bubbles” Su:2010qj (); FermiLAT:2014sfa (). It should be stressed that the prediction that we have constructed is absolute, and in fact it is remarkable that a very simple model such as the one we have constructed can reproduce the observations (in shape and absolute normalization) with an error of order 10–20%.
The reasonably good agreement between the model and the Fermi template for the energy gives support to the idea that hadronic mechanism is the main source of the diffuse gamma ray flux, and also indicates that the description of the interstellar gas density is adequate, however it is not sufficient to conclude that the factorization hypothesis for the CR spectra is correct. This is because the gamma ray flux at a single energy is sensitive to the spectrum of primary nucleons in a rather narrow range of energy (–30 ), and in our calculation the absolute normalization of the CR spectra can have a non trivial space dependence.
To test the factorization hypothesis encoded in Eq. (11) one needs to study the diffuse gamma ray in a broad energy range. The results of a calculation of the gamma ray up to very high energy performed assuming the validity of the factorization hypothesis are shown in Fig. 9. The top panel in the figure shows the gamma ray flux integrated in different angular regions of the Galactic plane () and plotted versus the energy. The absolute value of the flux depends on the angular region, being largest for regions toward the Galactic Center, and decreasing for larger longitudes , however the spectral shape of the flux remains approximately equal for TeV. At higher energy the absorption effects distort the spectra in an angle dependent way. The angle dependence of the absorption effects can be seen more clearly in the bottom panel of Fig. 9 that shows the average survival probability of the gamma rays, plotted as a function of , for the different angular regions.
The spectral distortion exibits a structure that reflects the fact that the absorption is generated by the interactions on two main components of the radiation field, dust emitted infrared radiation that is most effective for TeV, and the cosmic microwave background radiation (CMBR), that is most effective for PeV (see discussion in Vernetto:2016alq ()). The distortion pattern that is formed on the spectrum is qualitatively similar in different angular regions, however the amount of absorption is largest for directions toward the Galactic Center and minimum for directions toward the Anticenter. This can be easily understood, noting that the flux in directions toward the center has its origin in points that are on average further away from the Earth.
The absorption effects are illustrated in a complementary way in Fig. 10 that shows the longitude dependence of the flux, after integration in latitude in the range , for three values of the energy: GeV, where absorption is completely negligible, PeV where absorption is significant, and PeV where absorption is largest. In the figure the gamma ray flux is rescaled to have unit value at , for a better visualization of the absorption effects. As it is intuitively obvious, the flux in directions toward the Galactic Center is more suppressed by absorption than the flux toward the Anticenter.
It should be noted that the effects of absorption remain always smaller than a factor of order two, even in the case where they are most important, that is for of order 1–3 PeV, and directions toward the Galactic Center.
Vii Model 2: space dependent CR spectra
If one or more of the conditions listed in section VI are non satisfied, the spectra of cosmic rays can have a space dependent shape. Most models for the spectra assume that this is the case because the particles can lose a significant amount of energy propagating from the sources to distant regions of the Galaxy. For protons and nuclei, that have a much smaller , energy loss effects are are expected to be negligible, but a space dependence of the spectral shape can be generated by other mechanisms.
Some recent analysis of the Galactic diffuse flux Gaggero:2014xla (); Acero:2016qlg (); Yang:2016jda () conclude that there is some evidence for the fact that cosmic rays in the central part of the Galaxy have a harder spectrum than what is observed at the Earth, while cosmic rays in the periphery of the Galaxy are (moderately) softer. This effect can be described as a space dependence of the spectral index of the gamma ray emission. Fig. 11 shows some estimates of the dependence of the gamma ray emission spectral index on the distance from the Galactic Center. It has to be noted that a crucial problem in establishing the existence of these effects is to take into account the contribution of unresolved discrete Galactic sources. This problem will be discussed in section IX.
Aiming at the construction of a model as simple as possible we have assumed that the spectral index at the refence energy GeV has a space dependence specified by the functional form:
(17) 
The parameters in this equation have been chosen so that the spectral index at the Galactic Center has the value , while the index at large is . The dependence is described by the parameter kpc. This corresponds to a spectral index in the vicinity of the solar system . The studies of Gaggero:2014xla (); Acero:2016qlg (); Yang:2016jda () have not observed a dependence of the spectral index. This can be expressed as a lower limit kpc.
The gamma ray emission in the model is then described by the expression:
(18) 
It is straightforward to see that this model is identical to the one discussed in section VI for when the last factor in equal to unity, however for this model and the factorized model of Eq. (12) start to diverge.
This effect is illustrated in Fig. 12 that shows the gamma ray energy spectra calculated in the two models after integration in two angular regions toward the Galactic Center (top panel) and the Anticenter (bottom panel). For the angular region around the Galactic Center ( and ), the spectrum calculated in the non–factorized model is significantly harder, and the ratio between the two models grows with energy. The non–factorized model becomes a factor of ten larger for PeV. On the contrary, for the angular region around the Galactic Anticenter, the non–factorized model has a spectrum that is slightly softer. In this case the difference between the models is smaller (of order 20% for energies of order 1 PeV).
These points are also illustrated in Fig. 13, that shows the ratio of fluxes calculated in the two models for the two regions discussed above, and also a third intermediate region ( and ). In this third region the non–factorized model is moderately harder, with a ratio of order two in the PeV energy range.
The same information can of course be obtained studying the shape of the angular distribution of the diffuse flux at different energies in the two models. As discussed in the previous section, in a factorized model the angular distribution is energy independent, except for absorption effects. For a non–factorized model, such as the one we have constructed here, the enhancement of the flux from directions toward the Galactic Center becomes more and more significant with increasing energy. This is illustrated in Fig. 14, where the top panel shows the shapes of the longitude distribution of the gamma ray flux at energy of 1.8 PeV, in the two models. The ratio between the fluxes in the directions around the Galactic Center and Anticenter is one order of magnitude larger in the non–factorized model.
The survival probabilities for the two models are shown in the bottom panel of Fig. 14. The two probabilities are close to each other, but not identical reflecting the difference in the space distribution of the emission. This difference can be visualized inspecting Fig. 15 that shows the distribution of pathlength of the photons that form the diffuse Galactic emission at the Earth. The figure clearly shows how a very broad range of pathlengths contribute to the diffuse flux. In the non–factorized model, the contribution to the flux of points in the central region of the Galaxy becomes enhanced with increasing energy.
Viii The IceCube neutrino signal
As discussed in the introduction, the IceCube neutrino telescope has recently obtained evidence for the existence of a signal of high energy events of astrophysical origin above the expected foreground of atmospheric ’s Aartsen:2013jdh (); Aartsen:2014gkd (); Aartsen:2015rwa (); Aartsen:2017mau (). The signal is consistent with an isotropic flux of extragalactic neutrinos, generated by the ensemble of all (unresolved) sources in the universe. The flavor composition of the events in the signal (with the three flavors having approximately the same flux) is consistent with the expected composition of a flux generated by the standard mechanism of pion decay, after taking into account flavor oscillations (and averaging over a broad range of pathlengths).
Power law fits to the neutrino energy spectrum in the range – TeV, performed under the hypothesis that the signal is an isotropic extragalactic flux, have been recently presented by IceCube Aartsen:2017mau () for different classes of events and are shown in Fig. 16.
If the neutrinos of the IceCube signal are generated by a standard production mechanism, the emission should be accompanied by an emission of gamma rays with approximately equal spectral shape and normalization. If the neutrinos are extragalactic, one does not expect to observe an associated high energy photon flux because the gamma rays are (to a very good approximation) completely absorbed during propagation. On the other hand, if a significant fraction of the signal is of Galactic origin, the corresponding gamma rays flux is only partially absorbed and remains observable.
In Fig. 16 the IceCube fits to the neutrino spectrum are shown together with the measurements of the extragalactic and diffuse Galactic gamma ray fluxes obtained by Fermi, and also with the extrapolations of the diffuse Galactic flux (for the factorized and non–factorized models) that are discussed in this paper. Note that the figure shows angle integrated fluxes, and that the Galactic gamma ray fluxes have a strong angular dependence.
The comparison of the and fluxes indicates that the IceCube signal is significantly higher than the diffuse Galactic flux predicted on the basis of “natural” extrapolations of the observations at lower energy, even if one allows for the possibility that the emission of gamma rays and neutrinos is harder in the central part of the Galaxy. Similar results for the diffuse Galactic neutrino flux have been obtained by Pagliaroli:2016lgg ().
Stringent limits on the flux of astrophysical neutrinos from the Galactic disk have been obtained by ANTARES Albert:2017oba ().
Several authors have however suggested that a significant fraction of (or even the entire) IceCube signal is of Galactic origin. This requires the introduction of some new mechanism for production to explain the higher normalization and the approximately isotropic angular distribution of the neutrino signal.
For example, Esmaili and Serpico have suggested that the neutrino signal is generated by the decay of a dark matter particle with a mass of a few PeV Esmaili:2013gha (); Esmaili:2014rma (). In this case the space distribution of the emission is proportional to the dark matter density. Esmaili and Serpico adopt a Navarro–Frenk–White density profile:
(19) 
with kpc and GeV/cm.
Taylor, Gabici and Aharonian Taylor:2014hya () have suggested that the IceCube signal is generated by the interactions of cosmic rays filling a very extended spherical halo around the disk of the Milky Way with a size of order 100 kpc. In the following we will model the space dependence of the emission as a simple gaussian: with kpc (so that kpc).
Lunardini et al. Lunardini:2013gva () have suggested that a significant fraction of the IceCube signal is generated in the Fermi bubbles Su:2010qj (); FermiLAT:2014sfa (). In the following we will model the bubbles as spheres of radius kpc, with centers above and below the Galactic Center at kpc. The emission density in the Fermi bubbles grows with radius, and has been fitted with the form .
Assuming that the gamma ray emission in the three models listed above (dark matter decay, large halo and Fermi bubbles) has the same spectrum for all points, it is straightforward to compute the angular distribution of the emission and the absorption probability as a function of the energy and direction.
The latitude and longitude distributions of the flux observable at the Earth (calculated neglecting the effects of absorptions) are shown in Fig. 17 for the three models, together with the distribution of the factorized model discussed previously.
The gamma ray survival probabilities for the three non–standard models (averaged over all directions) are shown in the top panel of Fig. 18. The probability of absorption is significantly larger than for conventional emission models, reflecting the fact that the average pathlength of the photons is longer. The pathlength distributions for the different models are shown in the bottom panel of Fig. 17.
It should be noted that the knowledge of the pathlength distribution is not sufficient to estimate exactly the gamma ray absorption effects because the absorption probability per unit length has a non trivial space dependence. However, the most important component of the target radiation fields (the CMBR) is homogeneous and isotropic, and therefore it is possible to quickly obtain a reasonably accurate estimate of the absorption probability considering only the CMBR and the pathlength distribution.
Ix Galactic point sources
The measurement of the diffuse Galactic flux requires the subtraction of the extragalactic flux and of the contribution of all Galactic sources.
The extragalactic flux is (to a very good approximation) isotropic and can be accurately measured in polar regions (at large latitudes ) where the Galactic foreground is small and then subtracted also in angular regions where the Galactic flux is dominant. The subtraction of the contribution of the Galactic sources is more problematic, because it requires an estimate of the contribution of unresolved sources.
Many discrete (point–like and extended) gamma ray sources have been identified by Fermi in the 0.1–1000 GeV energy range Acero:2015hja (); TheFermiLAT:2017pvy (), and at higher energy by Cherenkov telescopes and air shower arrays tevcat (); hesscat (); Bartoli:2013qxm (); Abeysekara:2017hyn ().
The cumulative flux of all detected Galactic sources estimated from the published source catalogs is shown in the top panel of Fig. 19 and is compared to our extrapolations of the diffuse Galactic flux. In the case of the two Fermi catalogs (3FGL Acero:2015hja () and 3FHL TheFermiLAT:2017pvy ()) the separation of Galactic and extragalactic sources is performed statistically. The flux from extragalactic sources is estimated from a study of the sources at large latitude (). This contribution is then assumed to be isotropic and subtracted from the flux of all sources in the Galactic equatorial region to estimate the Galactic component. To estimate the flux at higher energy we have summed the flux of 116 sources with in the TeVCat online catalog tevcat ().
The bottom panel of Fig. 19 gives the ratio between the sum of all discrete resolved sources and the diffuse Galactic flux. For the diffuse flux we use the Fermi template for GeV and at higher energy the two extrapolations calculated in this work.
The important point of Fig. 19 is that the flux of the ensemble of the detected discrete sources has a harder spectrum than the diffuse flux, and the ratio between the flux of the ensemble of all resolved sources and the diffuse flux grows with energy from a value of order 0.06 at GeV to a value of order 0.3–0.5 at 1 TeV.
The estimate of the contribution of unresolved discrete sources to the Galactic flux requires models for the luminosity function and space distribution of the sources. In the present work we will not address this problem, but we can note here that it is likely that also the ratio will grow with energy.
Disentangling the contributions from gamma ray emission in interstellar space (the diffuse component) from emission inside or in the vicinity of sources is therefore an important problem, and an incorrect estimate of the flux of unresolved sources can lead to incorrect conclusions. One should also note that the contribution of unresolved sources is expected to have a non trivial angular dependence, with most of the flux arriving from directions toward the Galactic Center.
The problem of separating the “source” and “diffuse” components of the Galactic emission is likely to become more difficult with increasing . At very high energy, when cosmic rays escape rapidly from the sources and then from the Galaxy, the separation of the two components could in fact become impossible.
X Measurements of the diffuse gamma ray flux at high energy
Some measurements of the TeV diffuse Galactic gamma ray flux have been obtained in the recent past by high altitude air shower detectors located in the Northern hemisphere.
After a long collection of upper limits, the first measurement of the diffuse Galactic emission has been reported by the Milagro detector Atkins:2005wu () that measured the flux from the Galactic plane region – and at energies above 3.5 TeV. The measurement was higher than the expectations, suggesting the possible existence of a “TeV excess”, perhaps connected to the “GeV excess” previously reported by EGRET (later recognized as an instrumental effect). It has been later suggested that this measurement could include the contributions of several discrete sources and therefore overestimate the diffuse flux Bartoli:2015era ().
The Milagro telescope has later measured the diffuse flux at a median energy of 15 TeV in the region – and Abdo:2008if (). This measurement is shown in the top panel of Fig. 20 where it is compared with the predictions of our models for the same angular region. The Milagro observation is consistent with both the factorized and non–factorized models, but the second one seems to be slightly favoured.
The ARGO–YBJ detector has measured the energy spectrum of the diffuse Galactic emission in the region –100 and for energies in the interval between GeV and TeV Bartoli:2015era (). To evaluate the diffuse emission, the known gamma ray sources in the region have been masked, however a small residual contribution from non resolved sources cannot be excluded. The ARGO–YBJ measurement is shown in the middle panel of Fig. 20 and compared with our models. The data points are consistent with the predictions, but also in this case the large error bars do not allow to discriminate between the two models that in this energy and angular region give predictions that are close to each other.
At higher energy ( TeV) the available measurements of the diffuse emission have only provided upper limits. The most stringent results have been obtained by CASA–MIA, and constrain both the isotropic emission Chantell:1997gs () and the emission from the Galactic plane Borione:1997fy (). The bottom panel of Fig. 20 shows the CASA–MIA flux upper limits for the Galactic region –200 and , in the energy interval between 140 and 1300 TeV, compared to our predictions. The CASAMIA limits are a factor 2–5 higher than our models (depending on the energy). It can be noted that measurements in this angular region do not allow to discriminate between the factorized and non–factorized models because the two predictions are very close to each other.
In summary, the existing measurements of the diffuse Galactic gamma ray flux above 1 TeV are consistent with our extrapolations of the Fermi observations, but are not capable to discriminate between the two models discussed in this work. Future measurements with improved sensitivity, and a more complete coverage of the Galactic plane have however the potential to reach more firm conclusions on this problem.
x.1 Detector sensitivity
The study of the diffuse gamma ray flux at very high energy is probably best performed with air shower detectors with the capability to discriminate between electromagnetic cascades generated by photons and hadronic cascades generated by protons and nuclei.
For an order of magnitude estimate of the sensitivity of an air shower detector to a diffuse gamma ray flux, one can note that the signal of gamma ray events from a regions in the sky of angular size around a point of celestial coordinates and with energy in an interval of size centered on is:
(20) 
In this equation, is the diffuse gamma ray flux, is the detector area, is the observation time (much longer than a sidereal day), is the gamma ray detection efficiency, and is an adimensional factor that takes into account the visibility of the sky region under study at the detector geographical position. The quantity takes into account the fraction of a sidereal day that a point of celestial declination spends in the zenith angle range () where observations are possible:
(21) 
In this equation is the hour angle, is the zenith angle of a point in the sky with declination as seen by a detector of latitude , and the factor accounts for the geometrical reduction of the detection area when the source point has zenith angle (assuming a flat, horizontal detector). The quantity ), independently from the detector latitude, satisfies the condition:
(22) 
that expresses the obvious fact that the integral of the exposure over the entire visible part of the celestial sphere is independent from the detector geographical position. The top panel of Fig. 21 shows the quantity as a function of declination for some examples of the detector latitude. It is obviously of great interest to study the gamma ray flux in the Galactic equatorial region, in a broad interval of longitude. This is better achieved combining the observations of more than one detector. The bottom panel of Fig. 21 shows the exposure factor for points on the Galactic equator and different detector locations.
The gamma ray signal is detected together with a background of events generated by cosmic rays:
(23) 
In this equation is the cosmic ray flux, and is the fraction of the cosmic ray showers that survives after cuts designed to select photon showers (for example muon multiplicity and/or structure of the shower front).
Note that in general an air shower detector will measure the primary particle energy using an observable, such as the shower size, that is correlated to . The correlation between the observable and the energy will be in general different for showers generated by gamma rays or protons/nuclei, and it is important to take into account this difference. For example, if the shower size is used as an estimate of the energy, the selection of showers in a fixed interval of size selects photons and proton shower of different energy. Showers generated by photons have on average a size larger than hadronic showers of the same primary energy. This effects reduces the background, and can be included in the definition of the factor of Eq. (23).
The requirement that an observable signal must be larger than the background fluctuations () results in the minimum detectable flux:
(24) 
For , the condition for the minimum detectable flux becomes simply:
(25) 
where is the minimum number of events for a detection.
In Fig. 22 the sensitivity estimated with Eqs. (24) and (25) is compared with the expected diffuse flux estimated in the present work. One can conclude that a detector with an area or order one km with an hadron rejection factor of order has the potential to perform very interesting studies up to energies of order of several PeV’s.
In this discussion we have assumed that the diffuse gamma ray signal from the desired angular region can be estimated subtracting a background that is measured observing other regions in the sky where the signal is absent (or much smaller). To study (or set limits to) a diffuse flux that is quasi isotropic, this method cannot be applied. In this case the identification of the gamma ray signal must rely on a (Montecarlo based) absolute prediction of the rejection power for hadronic showers after use of the appropriate gamma ray selection algorithms. In this case the sensitivity of a telescope could be limited by systematic uncertainties in the description of hadronic cascades in the atmosphere.
Xi Outlook
The extension of the observations of the Galactic diffuse gamma ray flux to higher energy, in the TeV and PeV range, is a very important scientific goal that can give essential information on Galactic cosmic rays.
In this work we have focused the discussion on the dominant contribution to the flux, generated by the hadronic mechanism. The study of this dominant component allows to measure the spectra of protons and nuclei in distant regions of the Galaxy, and to determine their space dependence. At the moment it is known that in the energy range 10–100 GeV the CR density near the Galactic Center is a factor two to three times larger than what is observed at the Earth, however the question of the space dependence of the shape of the spectra remains uncertain.
The ratio between the fluxes predicted in models where the CR spectra have identical shape in the entire Milky Way and in models where the spectra in the central region of the Galaxy are harder, can become as large as one order of magnitude for gamma rays in the PeV energy range and directions close to the Galactic Center. These effects can be studied by air shower arrays with sufficiently good detection capabilities (area km and rejection for hadronic shower ). The energy range TeV is particularly important because it allows to obtain information about CR in the “knee” region also for distant parts of the Galaxy.
Measuring the subdominant leptonic contribution to the diffuse Galactic emission is a difficult but very important task, that can give fundamental information on the space dependence of the spectra. The leptonic component could become observable in the flux at large Galactic latitudes,
The effects of absorption of high energy gamma rays propagating over Galactic distances are important, and are largest for of order 13 PeV. The mean free path has its minumum value (of order 7 Kpc) for PeV. Photons of this energy emitted near the Galactic Center have a survival probability of order 0.29. This implies that a large part of the Galactic volume is effectively unobservable with PeV gamma rays, however a large fraction of the diffuse flux has its origin at shorter distances, and therefore observations in the PeV range can give very important information on the CR space distributions. The Galactic Center itself can in fact be studied even if the flux is suppressed.
The measurement of the diffuse gamma ray flux in a large angular region is very important for an understanding of the space dependence of the CR spectra, of particular importance is to observe the entire Galactic disk. Since ground based detectors can view only a fraction of the sky, it is desirable to have multiple telescopes at different latitude.
The problem of disentangling the diffuse flux from the flux of unresolved discrete sources is not completely solved even at low energy, and it is likely that this question will be more important in an energy range (0.1–10 PeV) where little is now known about the sources. It should however be noted that also a measurement of the sum of the two (diffuse and source) components can be very useful to develop an understanding of the Galactic cosmic rays.
The same information of the CR spectra that can be inferred from the observation of gamma rays is also contained in the spectra of Galactic neutrinos. The simultaneous measurement of the gamma ray and neutrino fluxes is clearly very desirable and very important to study the existence of a variety of different effects, including the possibility of non standard mechanism of production, or non–standard propagation for ’s and/or ’s.
Acknowledgments. During the preparation of this work we benefitted from discussion with several colleagues. We are especially grateful to Ralph Engel, Tom Gaisser, Elisa Resconi, Todor Stanev, Francesco Villante, Francesco Vissani and Cao Zhen.
References
 (1) W.L. Kraushaar et al., Astrophys. J. 177, 341 (1972).
 (2) D.A. Kniffen et al., Astrophys. J. 186, L105 (1973).
 (3) H.A. MayerHasselwander et al. Astron. Astrophys. 104, 164 (1982).
 (4) S. D. Hunter et al., Astrophys. J. 481, 205 (1997). doi:10.1086/304012
 (5) M. Ackermann et al. [FermiLAT Collaboration], Astrophys. J. 750, 3 (2012) doi:10.1088/0004637X/750/1/3 [arXiv:1202.4039 [astroph.HE]].
 (6) F. Acero et al. [FermiLAT Collaboration], Astrophys. J. Suppl. 223, no. 2, 26 (2016) doi:10.3847/00670049/223/2/26 [arXiv:1602.07246 [astroph.HE]].
 (7) R. Atkins et al. [Milagro Collaboration], Phys. Rev. Lett. 95, 251103 (2005) doi:10.1103/PhysRevLett.95.251103 [astroph/0502303].
 (8) A. A. Abdo et al., Astrophys. J. 688, 1078 (2008) doi:10.1086/592213 [arXiv:0805.0417 [astroph]].
 (9) B. Bartoli et al. [ARGOYBJ Collaboration], Astrophys. J. 806, 20 (2015) doi:10.1088/0004637X/806/1/20 [arXiv:1507.06758 [astroph.IM]].
 (10) A. Abramowski et al. [H.E.S.S. Collaboration], Phys. Rev. D 90, no. 12, 122007 (2014) doi:10.1103/PhysRevD.90.122007 [arXiv:1411.7568 [astroph.HE]].
 (11) M. C. Chantell et al. [CASAMIA Collaboration], Phys. Rev. Lett. 79, 1805 (1997) doi:10.1103/PhysRevLett.79.1805 [astroph/9705246].
 (12) A. Borione et al., Astrophys. J. 493, 175 (1998) doi:10.1086/305096 [astroph/9703063].
 (13) W. D. Apel et al. [KASCADE Grande Collaboration], Astrophys. J. 848, no. 1, 1 (2017) doi:10.3847/15384357/aa8bb7 [arXiv:1710.02889 [astroph.HE]].
 (14) M. G. Aartsen et al. [IceCube Collaboration], Science 342, 1242856 (2013) doi:10.1126/science.1242856 [arXiv:1311.5238 [astroph.HE]].
 (15) M. G. Aartsen et al. [IceCube Collaboration], Phys. Rev. Lett. 113, 101101 (2014) doi:10.1103/PhysRevLett.113.101101 [arXiv:1405.5303 [astroph.HE]].
 (16) M. G. Aartsen et al. [IceCube Collaboration], Phys. Rev. Lett. 115, no. 8, 081102 (2015) doi:10.1103/PhysRevLett.115.081102 [arXiv:1507.04005 [astroph.HE]].
 (17) M. G. Aartsen et al. [IceCube Collaboration], arXiv:1710.01191 [astroph.HE].
 (18) D. Gaggero, A. Urbano, M. Valli and P. Ullio, Phys. Rev. D 91, no. 8, 083012 (2015) doi:10.1103/PhysRevD.91.083012 [arXiv:1411.7623 [astroph.HE]].
 (19) R. Yang, F. Aharonian and C. Evoli, Phys. Rev. D 93, no. 12, 123007 (2016) doi:10.1103/PhysRevD.93.123007 [arXiv:1602.04710 [astroph.HE]].
 (20) S. Vernetto and P. Lipari, Phys. Rev. D 94, no. 6, 063009 (2016) doi:10.1103/PhysRevD.94.063009 [arXiv:1608.01587 [astroph.HE]].
 (21) I. V. Moskalenko, T. A. Porter and A. W. Strong, Astrophys. J. 640, L155 (2006) [astroph/0511149].
 (22) H. S. Ahn et al., Astrophys. J. 714, L89 (2010) doi:10.1088/20418205/714/1/L89 [arXiv:1004.1123 [astroph.HE]].
 (23) O. Adriani et al. [PAMELA Collaboration], Science 332, 69 (2011) doi:10.1126/science.1199172 [arXiv:1103.4055 [astroph.HE]].
 (24) M. Aguilar et al. [AMS Collaboration], Phys. Rev. Lett. 114, 171103 (2015) doi:10.1103/PhysRevLett.114.171103.
 (25) M. Aguilar et al. [AMS Collaboration], Phys. Rev. Lett. 115, no. 21, 211101 (2015) doi:10.1103/PhysRevLett.115.211101.
 (26) P. Lipari, Astropart. Phys. 97, 197 (2018) doi:10.1016/j.astropartphys.2017.11.008 [arXiv:1707.02504 [astroph.HE]].
 (27) P. Lipari, Phys. Rev. D 95, no. 6, 063009 (2017) doi:10.1103/PhysRevD.95.063009 [arXiv:1608.02018 [astroph.HE]].
 (28) Y. S. Yoon et al., Astrophys. J. 839, no. 1, 5 (2017) doi:10.3847/15384357/aa68e4 [arXiv:1704.02512 [astroph.HE]].
 (29) J. J. Engelmann, P. Ferrando, A. Soutoul, P. Goret and E. Juliusson, Astron. Astrophys. 233, 96 (1990).
 (30) T. K. Gaisser, T. Stanev and S. Tilav, Front. Phys. (Beijing) 8, 748 (2013) doi:10.1007/s1146701303197 [arXiv:1303.3565 [astroph.HE]].
 (31) K. M. Ferriere, Rev. Mod. Phys. 73, 1031 (2001) doi:10.1103/RevModPhys.73.1031 [astroph/0106359].
 (32) R. J. Glauber and G. Matthiae, Nucl. Phys. B 21, 135 (1970) doi:10.1016/05503213(70)905110.
 (33) T. Sjostrand, S. Mrenna and P. Z. Skands, JHEP 0605, 026 (2006)
 (34) O. Adriani et al. [PAMELA Collaboration], Phys. Rev. Lett. 106, 201101 (2011) doi:10.1103/PhysRevLett.106.201101 [arXiv:1103.2880 [astroph.HE]].
 (35) O. Adriani et al. [PAMELA Collaboration], Phys. Rev. Lett. 111, 081102 (2013) [arXiv:1308.0133 [astroph.HE]].
 (36) O. Adriani et al., Phys. Rept. 544, 323 (2014) doi:10.1016/j.physrep.2014.06.003.
 (37) S. Abdollahi et al. [FermiLAT Collaboration], Phys. Rev. D 95, no. 8, 082007 (2017) [arXiv:1704.07195 [astroph.HE]].
 (38) M. Aguilar et al. [AMS Collaboration], Phys. Rev. Lett. 113, 221102 (2014) doi:10.1103/PhysRevLett.113.221102.
 (39) F. Aharonian et al. [HESS Collaboration], Phys. Rev. Lett. 101, 261104 (2008) [arXiv:0811.3894 [astroph]].
 (40) F. Aharonian et al. [HESS Collaboration], Astron. Astrophys. 508, 561 (2009) [arXiv:0905.0105 [astroph.HE]].
 (41) D. Kerszberg for the HESS collaboration, in Proc. of the 35th ICRC (2017).
 (42) D. Borla Tridon et al. [MAGIC Collaboration], arXiv:1110.4008 [astroph.HE].
 (43) D. Staszak [VERITAS Collaboration], in Proc. of the 34th ICRC (2015), arXiv:1508.06597 [astroph.HE].
 (44) G. Ambrosi et al. [DAMPE Collaboration], Nature 552, 63 (2017) doi:10.1038/nature24475 [arXiv:1711.10981 [astroph.HE]].
 (45) G. R. Blumenthal and R. J. Gould, Rev. Mod. Phys. 42, 237 (1970) doi:10.1103/RevModPhys.42.237.
 (46) J. Roman–Duval, et al. ApJ 818, 144 (2016) doi:10.3847/0004637X/818/2/144.
 (47) M. Heyer and T.M. Dame, Ann. Rev. Astr. Astroph., 53, 583 (2015) doi:10.1146/annurevastro082214122324.
 (48) K. Ferriere, W. Gillard and P. Jean, Astron. Astrophys. 467, 611 (2007) doi:10.1051/00046361:20066992 [astroph/0702532].
 (49) P. M. W. Kalberla and L. Dedes, Astron. Astrophys. 487, 951 (2008) doi:10.1051/00046361:20079240 [arXiv:0804.4831 [astroph]].
 (50) J. H. Taylor and J. M. Cordes, Astrophys. J. 411, 674 (1993) doi:10.1086/172870.
 (51) J. M. Cordes and T. J. W. Lazio, astroph/0207156.

(52)
https://fermi.gsfc.nasa.gov/ssc/data/access/lat
/BackgroundModels.html
We have used the models contained in the file gll_iem_v06.fits.  (53) M. Su, T. R. Slatyer and D. P. Finkbeiner, Astrophys. J. 724, 1044 (2010) doi:10.1088/0004637X/724/2/1044 [arXiv:1005.5480 [astroph.HE]].
 (54) M. Ackermann et al. [FermiLAT Collaboration], Astrophys. J. 793, no. 1, 64 (2014) doi:10.1088/0004637X/793/1/64 [arXiv:1407.7905 [astroph.HE]].
 (55) M. Ackermann et al. [FermiLAT Collaboration], Astrophys. J. 799, 86 (2015) doi:10.1088/0004637X/799/1/86 [arXiv:1410.3696 [astroph.HE]].
 (56) G. Pagliaroli, C. Evoli and F. L. Villante, JCAP 1611, no. 11, 004 (2016) doi:10.1088/14757516/2016/11/004 [arXiv:1606.04489 [astroph.HE]].
 (57) A. Albert et al. [ANTARES Collaboration], Phys. Rev. D 96, no. 6, 062001 (2017) doi:10.1103/PhysRevD.96.062001 [arXiv:1705.00497 [astroph.HE]].
 (58) A. Esmaili and P. D. Serpico, JCAP 1311, 054 (2013) doi:10.1088/14757516/2013/11/054 [arXiv:1308.1105 [hepph]].
 (59) A. Esmaili, S. K. Kang and P. D. Serpico, JCAP 1412, no. 12, 054 (2014) doi:10.1088/14757516/2014/12/054 [arXiv:1410.5979 [hepph]].
 (60) A. M. Taylor, S. Gabici and F. Aharonian, Phys. Rev. D 89, no. 10, 103003 (2014) doi:10.1103/PhysRevD.89.103003 [arXiv:1403.3206 [astroph.HE]].
 (61) C. Lunardini, S. Razzaque, K. T. Theodoseau and L. Yang, Phys. Rev. D 90, no. 2, 023016 (2014) doi:10.1103/PhysRevD.90.023016 [arXiv:1311.7188 [astroph.HE]].
 (62) F. Acero et al. [FermiLAT Collaboration], Astrophys. J. Suppl. 218, no. 2, 23 (2015) doi:10.1088/00670049/218/2/23 [arXiv:1501.02003 [astroph.HE]].
 (63) M. Ajello et al. [FermiLAT Collaboration], Astrophys. J. Suppl. 232, no. 2, 18 (2017) doi:10.3847/15384365/aa8221 [arXiv:1702.00664 [astroph.HE]].
 (64) http://tevcat.uchicago.edu

(65)
https://heasarc.gsfc.nasa.gov/W3Browse
/all/hesscat.html  (66) B. Bartoli et al. [ARGOYBJ Collaboration], Astrophys. J. 779, 27 (2013) doi:10.1088/0004637X/779/1/27 [arXiv:1311.3376 [astroph.HE]].
 (67) A. U. Abeysekara et al., Astrophys. J. 843, no. 1, 40 (2017) doi:10.3847/15384357/aa7556 [arXiv:1702.02992 [astroph.HE]].