The Steady-State Multi-TeV Diffuse Gamma-Ray Emission Predicted with GALPROP and Prospects for the Cherenkov Telescope Array

For an assumed propagation phenomenology, the critical inputs for a GALPROP run are the CR source density distribution, together with the distributions of the other ISM components that result in the energy losses and corresponding secondary particle and emissions production. The CR source distribution is subject to considerable uncertainty, and we construct a representative set of 3D models varying the respective content of sources distributed in the Galactic disc and spiral arms. Meanwhile, GALPROP has well developed 3D models for the interstellar gas and radiation fields, and a collection of GMF models taken from the literature. We fix the interstellar gas model and make predictions for combinations of the interstellar radiation and magnetic field distributions over the CR source density distributions. We describe in detail below these elements for our calculations.

The propagation model parameters for each source distribution (Table 1) are used to calculate steady-state CR intensity solutions across the Galaxy using GALPROP. The spatial grid used for the CR propagation calculations is also employed for determining $\gamma$-ray emissivities. The $\gamma$-ray intensity maps at the solar system location are obtained by line-of-sight (LOS) integration of the $\gamma$-ray emissivities for the standard processes ($\pi^{0}$-decay, IC scattering), where the emissivities are determined for a logarithmic grid from 1 GeV to 100 TeV using five bins per decade spacing. We use a HEALPix (Górski et al., 2005) order 9 isopixelisation for the skymap generation²²2This gives a pixel size of $6.9^{\prime}$, which is similar to the point-spread function of H.E.S.S. ($4.8^{\prime}$).. All calculations of the IC component use the anisotropic scattering cross section (Moskalenko & Strong, 2000), which accounts for the full directional intensity distribution for the R12/F98 ISRF models. The $\gamma\gamma\rightarrow e^{\pm}$ attenuation that affects the $\gamma$-ray intensities mainly $\gtrsim$10 TeV, is included in the LOS integration, where the optical depth calculation includes the directionality of the ISRF (Moskalenko et al., 2006; Porter et al., 2018). For the R12/F98 ISRF models, the relevant precomputed optical depth maps included with the v57 GALPROP release are employed.

Fig. 1 shows the predicted fluxes from 1 GeV to 100 TeV averaged over different Galactic quadrants about the plane and within $30^{\circ}$ of the north/south polar regions. For clarity, we show only the modelled fluxes for the SA0 and SA100 distributions for the R12/PBSS ISRF/GMF combination (other combinations display qualitatively similar trends). Along the plane, for the central longitudes ($l\leq|90^{\circ}|$) and central latitudes ($b\leq|10^{\circ}|$) the IC emission starts to account for a large fraction of the total $\gamma$-ray emission around 100 GeV, and dominates the TeV Galactic emission after 10 TeV, for all source distributions. For the outer-Galactic region, where there are fewer CR sources ($l\geq|90^{\circ}|$), the IC provides a strong relative contribution until $\sim$10 TeV energies, when there is less variation overall between the SA distributions. Meanwhile, towards the polar regions ($b\geq|60^{\circ}|$) there is a much lower column density for the ISM gas, which results in a reduction of the $\pi^{0}$-decay emissions intensity. Because the scale height of the ISRF is larger than that of the gas, the IC emission is reduced to a lesser degree over these parts of the sky. This results in it being the strongest component of the $\gamma$-ray sky brightness for all SA distributions across the entire spectrum, from 1 GeV to 100 TeV.

Fig. 2 shows the total flux for the various SA distributions averaged over the HGPS analysis region ($l\leq|40^{\circ}|$ and $-1.5^{\circ}<b<1.0^{\circ}$). As the cooling time for VHE CR electrons is short, the impact of the source distribution is most easily seen for the IC emission contribution. Source density distributions with a higher weighting for CR sources in the arms generally boost the IC signal. This is due to the greater concentration of CR injection about the arms together with the higher ISRF/GMF densities that are also similar in distribution to the sources. This ‘density squared’ effect (Porter et al., 2017) enhances the emissions over the more smoothly distributed model combinations. The impact is such that for the SA100 distribution the IC flux at 100 TeV is a factor $\sim$2 brighter than that from the SA0 distribution.

Enhancements due to the ISRF distribution mainly affect the IC emissions for energies $\lesssim$10 TeV. In Fig. 3 we show predicted fluxes over the HGPS region with variation for the ISRF model. Below $\sim$100 GeV the stronger optical ISRF about the spiral arms in the R12 model enhances the IC emissions when compared to the F98 model. At higher energies the Klein-Nishina (KN) suppression reduces this effect, and the predictions for both become closer due to how the IR component of the respective ISRF models is distributed relative to the CR electrons. For energies $\sim$1–2 TeV the IC emission for the SA0 distribution varies $\sim$10% between the R12 and F98 ISRFs, while there is little variation switching between the ISRF models for the SA100 distribution. This variation for the SA0 distribution is due to the peaking of the source density about the GC. There is a correspondingly higher electron intensity across the inner Galaxy that can produce IC emissions. While for the SA100 distribution there is little variation as the inner cut-off for the spiral arms effectively produces a central ‘hole’ in the electron intensity distribution. Therefore, there is much less IC emission coming from the inner Galaxy for this source density distribution. For energies $\gtrsim$10 TeV even the IR component becomes KN supressed as well, and only scattering of the CMB is contributing to the emission. The result is a negligible difference between the ISRF models at higher energies, where the different predictions is due solely to the varying source density distributions.

Fig. 4 shows the variation of the IC flux when the GMF model is varied for fixed ISRF distribution. Varying the GMF distribution has a significant effect on the IC emission as energy increases into the VHE range due to the synchrotron losses dominating the cooling time scale for the electrons. At 1 TeV the difference in the IC emission between the PBSS and GASE distributions is approximately 30%, increasing to 100% by 100 TeV. For synchrotron emission, the energy radiated by TeV electrons over time ($\dot{E}/E$) is proportional to the square of the magnetic field strength ($B$) of the medium that they diffuse through and the energy of the particle, i.e. $\dot{E}/E\propto E^{2}B^{2}$. As the radiative losses of electrons depends strongly on the magnetic field strength, regions with higher GMF energy density will rapidly be depleted of high energy electrons. The GASE distribution has the most intense region of the field at the GC, leading to the IC emission being greatly reduced in the GC when compared to the PBSS distribution.

Fig. 5 shows the longitudinal profile of the IC and $\pi^{0}$-decay emission for the various source distributions for energy levels between 0.1 TeV and 100 TeV. The $\pi^{0}$-decay dominates the $\gamma$-ray emission for all longitudes for the energies $E<$1 TeV. Because the energy losses for CR nuclei are slow and the diffusion is relatively fast, the propagated distributions for them are fairly smoothly distributed over all source distributions. There is, correspondingly, minimal impact on the $\pi^{0}$-decay component in the longitudinal profiles. Meanwhile, the IC emission becomes dominant for $\gtrsim$10 TeV energies. Because the energy losses are much faster for the electrons, its contribution to the profiles is much more sensitive to changes in the source distribution over most of the Galactic plane. As the majority of the spiral arms are located between Earth and the GC, with one of the arms located in close proximity to the Earth, we see that for a larger fraction of CRs being injected into the spiral arms (e.g. SA100) the IC component becomes more intense for the central region (i.e. $l\leq|50^{\circ}|$). The IC emission is also boosted along the spiral arm tangents located at $l\approx|20^{\circ}|$. For the negative longitudes between the anti-centre and perpendicular to the GC, i.e. for $-90^{\circ}<l<-180^{\circ}$, the IC emission decreases with an increasing fraction of CR injection in the spiral arms, and the difference between the source density distributions increases with increasing energy.

Across the entire plane for all source distributions the IC emission dominates the $\gamma$-ray emission at 100 TeV. It is not a surprise that the spiral structure in the source density distributions including arms is evident in the profiles. However, at these energies there also appears to be structure for the SA0 density distribution IC profile similar to those models including arms that is not apparent at lower energies. Because of the KN suppression for the optical and IR components of the ISRF, the only photon field is the (uniform) CMB for IC processes. The energy density of the GMF varies, but is much higher than that for the CMB about the spiral arms, and hence determines the energy loss time scale affecting the electron distribution. The observed structure in the 100 TeV IC profile for the SA0 distribution is therefore entirely due to that encoded in the PBSS GMF distribution, with its influence imprinted on the electron energy density.

For the spiral arms in the source distributions, the closest arm in the sky region $-90^{\circ}<l<-180^{\circ}$ is the Perseus arm located $\sim$3 kpc away. For the source distributions with a high fraction of CR injection in the spiral arms (e.g. SA75 and SA100), there is a large void of electron sources between the Solar System and the Perseus arm, resulting in reduced TeV IC emission in that direction. However, as the Galactic disc component of the source distributions is smoothly varying, the source distributions with a low fraction of CR injection in the spiral arms (e.g. SA0 and SA25) have boosted TeV IC emission in these regions, as CRs are injected in the voids between the Solar System and the Perseus arm. This is not seen in the positive longitudes far from the GC ( $90^{\circ}<l<180^{\circ}$) as the Perseus arm is close enough to the Solar System such that there is no void in the electron densities in any of the source distributions. The proximity of the Perseus arm leads to a lower degree of variation in the IC emission in this region of the Galactic plane across the SA distributions. Figs. 4 and 5 therefore show that the source distributions with a larger fraction of sources in the spiral arms show increased IC emission along the arm tangents, with the IC emission dominating over the pion-decay at 1 TeV along the arms.

Fig. 6 shows the envelope of longitudinal profiles at 1 TeV of all the GALPROP runs across our model grid, i.e. any combination of the source distribution, ISRF, and GMF. Due to using a mix of distributions that include the spiral arms and those that do not (e.g. combining the SA0 with the PBSS distribution), no single GALPROP run describes the lower or upper bound of the envelope. However, from the upper/lower bounds we infer that there is a $\pm 50$% variation over the entire grid of GALPROP predictions for models that are consistent with the local CR data. Changing the source distribution alone tends to alter the longitudinal profile most significantly toward the spiral arm tangents, where the variation can be up to $\sim$40% when comparing SA0 to SA100. Toward the inner Galaxy there are smaller differences, due to both the lower number of sources about the GC for distributions with spiral arms, and to the smoothing from the LOS integration over the long path lengths for these directions. For similar reasons, changing the ISRF model impacts the results only modestly (up to $\sim$4%) around the central plane $|l|\leq 30^{\circ}$. Meanwhile, changing the GMF model has the strongest effect around the central $|l|\leq 100^{\circ}$ where the impact can be up to $\sim$20% variation. We show overlaid in Fig. 6 two individual results – one with no spiral arms in the CR injection distribution, ISRF, or GMF (i.e. the SA0, F98, and GASE configuration), and the other having spiral arms in all three (i.e. the SA100, R12, and PBSS configuration). The variation between the predictions is generally driven by the increasingly complex 3D structures for the respective components of the models that we have considered.

Our work is the first systematic study using the GALPROP framework to predict the diffuse emissions into the TeV energy range for a grid of physical modelling configurations. For the steady-state models that we have considered, generally the $\pi^{0}$-decay component produced by the CR nuclei interactions with the interstellar gas has only a small variation over the HGPS survey area (Fig. 2). Because the energy losses of these CRs are slow, they are distributed fairly smoothly across the MW with only mild imprints of the source density distribution. However, this is not the case for the IC emissions. Averaged over the HGPS survey area there is a much stronger variation due to the choice of source distribution. The energy losses for the CR electrons are much more rapid than those for the nuclei. The concentration of the ISRF and GMF energy densities about the spiral arms combines to affect the $\gamma$-ray intensities toward those regions, because they approximately coincide toward where the source densities with arms tend to also be enhanced. Into the TeV energy range, the IC emissions have much higher sensitivity to the ISRF and GMF energy density distributions.

There are three composite source density distributions in our source density model grid (both Galactic disc and spiral arm components: SA25, SA50, and SA75). For a fixed target density distribution (i.e. holding the ISRF/GMF constant), comparing the 1 TeV $\gamma$-ray longitudinal profiles for these three composite source density distributions shows a variation of $\sim$10%. When comparing all five source density distributions, we find a much larger variation up to $\sim$50%. Altering the target densities as well as the source density distributions, as is done in Fig. 6, shows a variation up to $\sim$100%. The effect of progressively more complex 3D structures for the source density and ISRF/GMF components leads to more variation for the TeV emissions from the CR electrons for the propagation model phenomenology that we have used in this paper. However, the envelope of model predictions resulting from variations in the target densities is larger than that from solely changing the source density distribution. While our current understanding of the spatial distributions of the ISM components that are the targets for producing the VHE emissions is incomplete, we have used models for them that are reasonably agreeing with data. The GALPROP predictions from our study are therefore representative of the likely uncertainty associated with current modelling for the VHE diffuse emissions.

We have used a grid of source density models where the normalisation is obtained by requiring consistency of the propagated CR intensity distributions with the data obtained at the solar system. This is a well-motivated normalisation condition, given the considerable uncertainties associated with correctly defining CR source population characteristics. However, there is an effective CR ‘horizon’ beyond which additional components may be added without any noticeable effect when normalising to the local CR data. For example, CRs injected about the GC likely have a negligible impact on the intensities at the solar system location. Correspondingly, a source density localised about the inner Galaxy, e.g. distributed similar to the stellar bulge/bar (Porter et al., 2017; Macias et al., 2019), will have its injected CR power as a free parameter. Such additions localised about the inner Galaxy require tuning at other frequencies, e.g. the GeV energy data, to correctly optimise parameters for use as a ‘background’ in the VHE regime. While the uncertainty band for the diffuse emissions that we estimated above is below the 5$\sigma$ significance level of the HGPS (see Fig. 8), such additional components have the possibility to also contribute toward the inner Galaxy. At the moment it is likely that a unique solution is not possible for precise estimation of the diffuse background for the HGPS data.

Our modelling has also shown that the IC emissions across the plane starts to become dominant over the $\pi^{0}$-decay above 10 TeV (Fig. 5). The profiles for $\gtrsim$10 TeV $\gamma$-ray energies display structure in the leptonic component that is due to the combination of the spatial distributions employed for our models. However, at the higher energies ($\sim$100 TeV) the situation becomes simpler, because the complications of the uncertainty in the ISRF distribution effectively disappear by the KN suppression of the optical and IR components. Only the GMF and source density distributions are determining the structure visible in the longitudinal profiles for these energies. Currently, the GMF models for the MW are poorly constrained, and are most commonly constructed by using rotation measures of extragalactic sources. Given that the GMF structure is observable at 100 TeV, the future VHE/UHE $\gamma$-ray observations provide an interesting possibility to constrain the GMF configurations. However, $\gamma$-ray emission is created by the convolution of the source density distribution and the target density distributions. The alignment of the GMF and source density distribution is not complete, which is not surprising given that the different tracers for the Galactic spiral structure show offsets (e.g. Vallée, 2020). Therefore using $\gamma$-ray observations to constrain the GMF will likely require $\gamma$-ray facilities many times more sensitive than what is currently available.

We note that for this paper we have only considered steady-state models. Into the TeV energy range the time-dependent solutions, particularly for the CR electrons, are likely necessary (Porter et al., 2019). For example, the Tibet AS${\gamma}$ data indicate that many of the $\gamma$ rays detected within their $25^{\circ}<l<100^{\circ}$, $|b|<5^{\circ}$ window (Amenomori et al., 2021) are not close on the sky to known VHE $\gamma$-ray sources. The observational data from Tibet AS${\gamma}$ appears inconsistent with the simple models of Lipari & Vernetto (2018) and Cataldo et al. (2019)⁴⁴4These efforts do not use propagation codes, instead employing simple parametrizations for the CR spectra assuming steady-state conditions, and do not consider the 3D structure of the CR sources and ISM., and a discrete origin has been suggested (e.g. Dzhatdoev, 2021; Fang & Murase, 2021). Recent results from Vecchiotti et al. (2022) found that the discrepancy between models and the Tibet AS${\gamma}$ flux could be due to a population of unresolved PWNe, which would likely increase the leptonic component of the $\gtrsim 100$ TeV $\gamma$ rays. PWN halos have also been proposed as a large contributor to the unresolved sources for both H.E.S.S. and HAWC (Martin et al., 2022), further suggesting that electrons contribute to the diffuse emission. The increasing contribution by the leptonic component with energy that we have shown in this paper, and its likely significant fluctuations, may help explain the lack of correlation between the PeV $\gamma$ rays observed by Tibet AS${\gamma}$ and neutrinos seen along the Galactic plane. We also have results from HAWC (Abeysekara et al., 2017) and LHAASO (Zhao et al., 2021) that show for energies $\gtrsim$10 TeV the localisation of emissions about individual sources is stronger. Incorporating time-evolution into the modelling seems to be necessary next step to accurately connect the GeV to 100s TeV energies emissions that are presumably coming from a common origin.

CTA will perform a Galactic plane survey (GPS) over its first ten years of operation, which will cover Galactic latitudes in the range $b\leq|5^{\circ}|$ and all Galactic longitudes, averaging 32 hours of observation per square degree (CTA Consortium et al., 2019). The CTA GPS will have an energy range from 100 GeV to 100 TeV, with a resolution of around $3^{\prime}$ for energies above 1 TeV. The $5\sigma$ point source sensitivity of CTA’s proposed 10-year survey will reach 1.8 mCrab for the inner $l\leq|60^{\circ}|$ and 3.8 mCrab for the outer Galactic longitudes. For comparison, the HGPS has a point source sensitivity of 4–20 mCrab between $40^{\circ}\leq l\leq 80^{\circ}$ with an average of 146 hours of observation per square degree. The sensitivity of the proposed CTA GPS after adjusting for a containment radius of $R_{c}=0.2^{\circ}$ is also shown in Fig. 8.

As CTA is expected to detect and resolve on the order of a thousand of new sources (Actis et al., 2011), these exclusion regions will grow to encapsulate a significant fraction of the Galactic plane, if not the entire plane. Allowing nearby sources to contaminate the background, such as would be very likely along the Galactic plane without careful consideration, would significantly reduce the sensitivity of the observations (Ambrogi et al., 2016).

CTA presents a leap forward in terms of the angular resolution it will achieve, with a PSF reaching approximately $1.8^{\prime}$ for $\gamma$-ray measurements at 1 TeV. The all-sky H i and CO survey data that is used for building the neutral ISM model used in this work have resolutions of $27.5^{\prime}$ and $13.7^{\prime}$, respectively. This is adequate for comparisons to the HGPS, but not to CTA’s expected Galactic plane survey. New high-resolution ISM surveys, such as Mopra CO (Braiding et al., 2015, 2018) in the Southern hemisphere, and Nobeyama 45m CO (Minamidani et al., 2015; Umemoto et al., 2017) and THOR H i/OH (Beuther et al., 2016; Wang et al., 2020) in the Northern hemisphere have become available. High-resolution ISM surveys that can observe the small, dense clumps at comparable or better angular resolution to CTA will become a necessary modelling ingredient for future high-resolution comparisons to its GPS. Due to the high density of these clumps it will also become necessary to adapt the diffusion coefficient to the density of the ISM gas (e.g. Gabici et al., 2007, 2009).

From Fig. 8 we can see that CTA-South should be able to detect the TeV diffuse $\gamma$-ray emission to the 5$\sigma$ level given even our conservative estimates of the diffuse $\gamma$-ray emission predicted in this paper for the central $80^{\circ}$ of longitude ($|l|\leq 40^{\circ}$). As has already been seen with the sensitive Fermi-LAT data at lower energies, it will therefore be necessary to develop accurate models for the diffuse emission to separate individual source characteristics from it. This will be essential for resolving and detecting faint, extended TeV emissions such as those such as those expected from PWN haloes, and for probing complex morphological structures of TeV sources. This will potentially allow the identification of currently unidentified sources in the HGPS and other current and future observations, such as those from CTA.

For the Galactic longitudes open to Northern-hemisphere observations, i.e. $l>|50^{\circ}|$, our estimates state that the diffuse emission flux at 1 TeV is below 5 mCrab/deg². This emission is below the extended source sensitivity for CTA-North, which is planned to perform the deepest Galactic plane survey of any IACT in the Northern hemisphere at 1 TeV (CTA Consortium et al., 2019). It is therefore unlikely that any IACT in the Northern hemisphere, such as VERITAS (Vassiliev, 1999) and MAGIC (Aleksić et al., 2012), will observe the diffuse emission. Ground-level particle detectors in the Northern hemisphere, such as LHAASO (Sciascio, 2016), HAWC (Abeysekara et al., 2013), and Tibet AS${\gamma}$ (Sako et al., 2009), are optimised for multi-TeV to PeV observations. This has allowed Tibet AS${\gamma}$ to successfully detect extended, presumably diffuse, emission at PeV energies (Amenomori et al., 2021). However, as these particle detectors are optimised for multi-TeV to PeV observations, it is unlikely they will observe this diffuse emission around 1 TeV within their first ten years of operation. The Northern hemisphere detector ARGO-YBJ was constructed to be more sensitive at TeV energies, and has measured the diffuse emission at 1 TeV after five years in the region $25^{\circ}<l<100^{\circ}$. The emission predicted by our GALPROP models is within the uncertainty band estimated for these data (Bartoli et al., 2015).

Observation of the Galactic diffuse emission below the PeV regime (e.g. around 1 TeV) will likely remain difficult with the current and the next generation of $\gamma$-ray observatories. For the Southern hemisphere, the planned ground-level particle detector SWGO (Albert et al., 2019) has the possibility of observing the Galactic diffuse emission at 1 TeV around the GC within the first ten years of operation, complementing the higher-resolution CTA-South observations.