Galactic diffuse gamma rays meet the PeV frontier

The Tibet AS$\gamma$ and LHAASO collaborations recently reported the observation of a $\gamma$-ray diffuse emission with energy up to the PeV from the Galactic plane. We discuss the relevance of non-uniform cosmic-ray transport scenarios and the implications of these results for cosmic-ray physics. We use the {\tt DRAGON} and {\tt HERMES} codes to build high-resolution maps and spectral distributions of that emission for several representative models under the condition that they reproduce a wide set of local cosmic-ray data up to 100 PeV. We show that the energy spectra measured by Tibet AS$\gamma$, LHAASO, ARGO-YBJ and Fermi-LAT in several regions of interest in the sky can all be consistently described in terms of the emission arising by the Galactic cosmic-ray"sea". We also show that all our models are compatible with IceTop $\gamma$-ray upper limits. Our results favor transport models characterized by spatial-dependent diffusion although some degeneracy remains between the choice of the transport scenario and that of the cosmic-ray spectral shape above 10 TeV. We discuss the role of forthcoming measurements in resolving that ambiguity.


Introduction
The past few years witnessed remarkable observations in the field of multi-messenger astrophysics, especially in the cosmicray (CR) and γ-ray channels.The Tibet ASγ experiment reported the observation of a diffuse γ-ray emission from the Galactic plane (25 • < l < 100 • and 50 • < l < 200 • , both for |b| < 5 • ) with energy up to the PeV (Amenomori et al. 2021).This finding was confirmed by the LHAASO collaboration (Zhao et al. 2021) for the innermost region.If not originated in currentlyactive unresolved sources (see e.g.Casanova & Dingus 2008 and the more recent Vecchiotti et al. 2021 for an alternative interpretation), the presence of truly diffuse γ-rays at ∼ PeV is likely due to ∼ O(10) PeV CRs injected by PeV accelerators that were active in the past, the so-called PeVatrons.The ability to explore the knee region (E CR ∼ few PeV s) of the CR spectrum is of outstanding importance for our understanding of CR physics.Indeed, if we assume the conventional scenario of Supernova Remnants (SNRs) as the bulk of Galactic CRs (see Gabici et al. 2019 for a recent review), it is a theoretical challenge to even achieve particle acceleration at the level of ∼ O(100) TeV, and only for a short time in the evolution of the remnant (Lagage & Cesarsky 1983a,b).To overcome this problem, stellar clusters have recently come back as a viable explanation for such high-energy particle acceleration (Cesarsky & Montmerle 1983;Morlino et al. 2021), although it is not clear up to what extent in the locally observed CRs.Therefore, whether the knee in the CR spectrum is due to a change in the CR acceleration mechanism or to a transport effect is still matter of debate.
For what just described, the Tibet ASγ and LHAASO measurements offer a new valuable handle to study the origin and the propagation of Galactic CRs, which complement direct measurements, that above ∼ O(100) TeV are performed with Extended Air Shower (EAS) experiments.As witnessed by the large discrepancies in the energy spectra observed by different collaborations (see Fig. 1), these measurements suffer from large systematic errors, mostly associated with modelling hadronic interaction within the Earth atmosphere.Another valuable advantage of measuring the E γ ≥ TeV γ-ray diffuse emission of the Galaxy is that it may allow to test whether features in the CR spectrasuch as the hardening at ∼ 300 GeV/n found by CREAM (Ahn et al. 2010), PAMELA (Adriani et al. 2011), AMS (Aguilar et al. 2015c,a) and the softening at ∼ 10 TeV/n measured by DAMPE (An et al. 2019a) and CREAM (Yoon et al. 2017b) -are due to local sources or if they are representative of the whole Galaxy, namely due to transport effects (Blasi et al. 2012).Moreover, the detection of the diffuse emission at low Galactic longitudes gives crucial information on how CR propagation behaves in the inner regions of the Galaxy, where the conditions for CR transport are expected to be different than in the average disc, noticeably leading to a radial (Galactocentric) dependence for the CR spectra.It should be noted that, while in the last years a huge experimen-tal effort allowed to strongly reduce the uncertainties on the CR transport parameters within a few kpc's horizon around the Solar System, very little is still known about CR propagation beyond that distance.With this regard, the study of weakly interacting messengers such as γ-rays and neutrinos originated in the scattering of CRs with the InterStellar Medium (ISM) and the IS Radiation Field (ISRF) offers a valuable probe of CR transport in those distant regions.
In Gaggero et al. (2015b); Cerri et al. (2017) the authors argue that the hardening of the γ-ray diffuse emission found by Fermi-LAT above 10 GeV (Ackermann et al. 2012b;Acero et al. 2016) -see also Yang et al. (2016) -can be explained in terms of a radially dependent CR spectrum.Based on that interpretation, in Gaggero et al. (2015a) the authors predict a diffuse flux that is much larger than what is expected by conventional transport models at low Galactic longitudes and very high energies E 1 TeV.They show that this also allows to consistently reproduce Fermi-LAT and Milagro (Abdo et al. 2008) measurements in the inner Galaxy, as well as (Gaggero et al. 2017) the H.E.S.S. observed diffuse-emission in the Central Molecular Zone (Abramowski et al. 2016) (although the CR population in that small region may be affected by one or more PeVatrons close to the Galactic Center (GC)).
Recently, a phenomenological model inspired by that scenario -where the energy and spatial dependence of the CR spectrum are non-factorized -has been studied (Lipari & Vernetto 2018) (LV non-factorized model hereafter) and used by Tibet ASγ to interpret their results (Amenomori et al. 2021).For reference, the authors also considered a conventional (LV factorized) model with a radially independent CR spectrum.LV models also account for γ-ray absorption -which is significant above 100 TeV -due to pair production by scattering onto the CMB and the dust emitted infra-red photon background (Vernetto & Lipari 2016).Both LV models are built to reproduce AMS-02 and CREAM proton and Helium data below 100 TeV.As shown in Amenomori et al. (2021) the LV non-factorized (spatial dependent) model seems to be in better agreement with Tibet ASγ results than the factorized (conventional) one.Preliminary LHAASO data (Zhao et al. 2021), however, are slightly lower than Tibet in some energy bins.Most importantly, they have significantly smaller statistical errors and extend down to 10 TeV where they are less scattered and the uncertainties on the CR spectral elemental shapes are smaller.As we will show this may have relevant implications for their interpretation.
In this work we model the γ-ray diffuse emission of the Galaxy from 10 GeV up to the PeV trying to assess the effect of the uncertainties due to the choice of the CR spectra, hence on the choice of the transport scenario.Differently from what has been done in Lipari & Vernetto (2018), we make extensive use of Fermi-LAT results -showing their relevance to constrain CR transport properties -and use physically motivated models.In fact, CR transport is modelled with the DRAGON2 numerical code (Evoli et al. 2017(Evoli et al. , 2018) ) in two transport scenarios: either assuming a spatially independent or dependent diffusion coefficient.Additionally, we use the newly released HERMES code (Dundovic et al. 2021) to integrate along the line of sight the CR spatial and energy distributions obtained with DRAGON2 using detailed IS gas emission maps, updated ISRF models and γray production cross-sections, to get high resolution sky maps of the diffuse emission at all relevant energies, which allow to provide more accurate and robust results than with the analytical gas models used in (Lipari & Vernetto 2018).Moreover we account for a wider set of CR data above 10 TeV (see, however, Vernetto & Lipari 2021) and, in order to cope with the large discrepancies among different data sets (see also Koldobskiy et al. 2021), we consider two set-ups for the CR injection spectra.We will then compare our predictions with Tibet ASγ and LHAASO results as well as with ARGO-YBJ (Bartoli et al. 2015) at lower energy, and with IceCube (Aartsen et al. 2019b) and CASA-MIA (Borione et al. 1998) upper limits between 100 TeV and 2 PeV in regions closer to the GC, besides with Fermi-LAT data in several Galactic Plane (GP) quadrants.
We will show that Tibet ASγ and, if confirmed, LHAASO data in combination with Fermi-LAT favor a spatially dependent transport scenario.
As shown in Fig. 1, measurements are in good agreement with each other for both protons and Helium nuclei, and show a power-law structure up to ∼ O(100 − 500) GeV, where a spectral hardening is reported by all the experiments, typically interpreted as a consequence of a new phenomenon occurring in particle transport rather than in the acceleration (see Gabici et al. 2019 and references therein).Another structure that is clearly visible in the proton and helium spectra is a pronounced softening reported by DAMPE at a rigidity of ∼ O(10) TV (An et al. 2019b;Alemanno et al. 2021), confirming previous claims by CREAM-III (Yoon et al. 2017a) and NUCLEON (Atkin et al. 2018).
It is important to notice that, since the highest-precision measurements come from experiments that are not ground based -for which it is challenging to achieve energies above ∼ 10 TeV with significant statistics -, the highest energies currently reached do not allow to discriminate whether the softening is a (i) local feature or (ii) a global structure.In the first scenario (i), the particle spectra in the Galaxy would follow its trend after the hardening at a few hundreds of GeV and proceed up to the knee at ∼ 5 PeV, where it undergoes a cut-off in the form of a sharp softening -the "bump" is then typically interpreted as a signature of a nearby source contributing to the primary CRs (Fornieri et al. 2021).Throughout the paper, we refer to this setting as the Max setup.As for the second case (ii), we would expect a softening above E CR ∼ O(10) TeV and a cutoff as a pronounced break at the knee (at ∼ a few PeV).We refer to this other setting as the Min setup.
For what concerns the highest-energy band (1 PeV E γ 1 EeV), such energies can be achieved at the moment only by ground-based experiments measuring extensive air showers (EAS), which imply much larger uncertainties and strong dependence on models of the hadronization occurring in the Earth's atmosphere, specifically on their Monte Carlo-based simulations (see Pierog et al. 2006;Parsons et al. 2011 for a discussion on the model-dependent uncertainties of the hadronic showers).In particular, we are considering data points from KASCADE (Antoni et al. 2005), KASCADE-Grande (Apel et al. 2013) or, in alternative, from IceTop (Aartsen et al. 2019a), being the two datasets incompatible with each other.
In what follows, we compute the γ-ray emission in the Galaxy for both the Min and the Max spectral setups.As shown in Fig. 1, our procedure is addressed to bracket the uncertainties due to the poor knowledge of the CR spectrum in the highestenergy range.
As an additional comment, we do not account for CR species heavier than Helium since, under reasonable conditions, their contribution to the gamma-ray source term is subdominant (less than 10%) with respect to the expected flux due to proton + Helium production rate (see e.g.Breuhaus et al. 2022 for a recent analysis).

Fermi-LAT data
The LAT is a pair-conversion telescope that provides gammaray data from August 2008 in a range of energies ranging from ∼ 20 MeV to a few hundreds GeV (Atwood et al. 2009).In this work, we make use of ∼ 149 months of data (from 2008-8-2008 to 2020-12-31), selecting CLEAN events from the PASS8 data.We select events from good quality time intervals and remove the intervals when the LAT was subtended at rocking angles θ r < 52 deg ((DATA_QUAL>0) && (LAT_CONFIG==1) && ABS(ROCK_ANGLE<52)), and in addition apply the zenithangle cut θ z < 100 deg, in order to avoid contamination with the Earth albedo.We consider front-and back-converted events in order to minimize statistical errors.We employ the P8R3_CLEAN_V3 version of the instrument response functions.The extraction of Fermi data and calculation of exposure maps is performed using the most up-to-date version of the Sci-enceTools1 (2.0.8; released on 01/20/2021).In this work, our γray model is optimized to reproduce Fermi-LAT data in different windows of galactic longitude in the energy range between ∼ 10 to ∼ 300 GeV, under the hypothesis of a spatially-dependent spectral index of the diffusion coefficient (see the details in section 3.2), considering the recent Fermi 4FGL-DR2 catalog for accounting for the source's γ-ray emission.

Very High Energy gamma-ray data
A dedicated section is necessary for the very-high energy band, namely for the photons detected with energy E γ ≥ 100 TeV, that can be probed via the deconvolution of the EAS's at groundbased detectors.In particular, we want to comment on the observations of PeV γ-rays by the LHAASO (Cao et al. 2021) and Tibet ASγ (Amenomori et al. 2021) Collaborations.LHAASO reports the detection of ∼ 530 photons above E γ = 100 TeV at σ ≥ 7 significance from 12 regions with overlapping known sources.Among these, the source J2032+4102 -associated with the stellar cluster Cygnus-OB2, considered a probable PeVacceleration site (Aharonian et al. 2019) -deserves special at-tention, as it is where the highest-energy photon (E γ = 1.42 ± 0.13 PeV) is originated.Alongside, Tibet observes 10 events in the energy bin 398 ≤ E γ (TeV) ≤ 1000, 4 of which originate around the OB2 cluster.However, in order to remove the cluster contribution to the total detected photons, in Tibet they use a Gaussian profile with width σ OB = 0.5 • -this is consistent with the angular extension exploited by LHAASO (Cao et al. 2021) -, motivated by the profiles σ HAWC ≤ 0.3 • used by the HAWC Collaboration for sources above E γ = 56 TeV (Abeysekara et al. 2020).On the other hand, for the specific case of the Cygnus cocoon, HAWC reports an extension of σ OB, HAWC = 2.1 • (Abeysekara et al. 2021), which implies that most of the diffuse emission from the Cygnus region could actually be coming from the single source.This information has been used in Liu & Wang (2021) to interpret the neutrino counterpart expected from the region, that would otherwise be in tension with the IceCube observations on the GP (Aartsen et al. 2017).
As a result, it is possible that the highest energy γ-ray point in the Tibet observations actually should show a flux lower than what is currently estimated.Future analyses will certainly help to clarify the picture.

The models
The diffuse γ-ray sky is the superposition of a number of different components.They include the isotropic background: This component captures the contribution from extra-galactic diffuse emission, unresolved extra-galactic sources, and residual (mis-classified) cosmicray emission.We adopt here a reference isotropic spectral template provided by the Fermi-LAT collaboration (iso_P8R3_CLEAN_V3_v1)2 .the unresolved point-like, or extended, source emission: The contribution of unresolved sources is derived from population synthesis, using the model for the Galactic gamma-ray population following a four-arm spiral distribution described in Steppa & Egberts (2020).This source model for fluxes at 1 TeV is translated to the energy range of the Fermi-LAT data using a spectral index of -2.4 for the TeV energy range as the mean spectral index of the H.E.S.S. Galactic source catalog (Abdalla et al. 2018) and for lower energies a spectral index of -2.27 as mean of the Fermi-LAT 4FGL (Abdollahi et al. 2020) with a transition regime between 100 GeV and 1 TeV.Source populations are simulated and their sources classified as detectable or unresolved according to the Fermi-LAT sensitivity of the 4FGL.The sample of detectable sources is used as verification: both the cumulative spectrum and the log N − log S distribution of detectable sources are well matched in simulations and data in the energy range considered in the Fermi analysis.As model for the unresolved emission we use the mean of 200 realisations, which results in a flux in unresolved sources of around 1% of the detected sources.Individual realisations exhibit large fluctuations, manifested in a standard deviation in the simulated fluxes of > 100%.
the secondary emission of Galactic CRs: at the energies relevant for this paper (E 10 GeV) this can be further divided in the hadronic component, due the decay of neutral pion produced by CR scattering onto the IS gas (mostly hydrogen and Helium), and the Inverse Compton (IC) emission of CR electrons and positrons onto the ISRF.The relative contributions of these components depend on the Galactic coordinates and on the energy.On the GP and at energies larger than 10 TeV, the hadronic emission by CRs is expected to be dominant although a significant -see Linden & Buckman (2018) -contribution due to IC emissions cannot be excluded.
Here we focus mainly on modeling the secondary diffuse emission due to interaction of Galactic CRs during their propagation.We do that with the HERMES (Dundovic et al. 2021) code which, at each given energy bin and for each relevant CR species, performs a numerical integration along the line-of-sight of the product of the CR differential energy flux, of the IS gas density (or the ISRF for the IC emission) and of the γ-ray production cross section.More details on the cross-sections and the gas (Hydrogen and Helium) distributions used in this work will be given in Secs.3.2 and 3.1 respectively.
In the following subsection we rather discuss how the CR energy and spatial distributions are computed.

The interstellar gas
Our model consists of a set of column density maps in (l, b) Galactic coordinates for atomic and molecular gas, associated to Galactocentric rings.The atomic gas model is based on the 21cm line emission data observed by the recent HI4PI survey that covers the whole sky with a 1/12 degree binning (HI4PI Collaboration et al. 2016).As far as molecular gas is concerned, the decomposition is based on the observations of the CO rotational line at 115 GHz from the CfA survey (Dame et al. 2001;Dame & Thaddeus 2004).The profile decomposition is discussed in Remy et al. (2021); Fermi-LAT (2021).In our framework, every Galactocentric ring can be associated to a value of the CO-H 2 conversion factor (X CO ).In our model, we adopt the values of [1.8, 3.5, 4.0, 4.5, 7.5, 8.0] in units of 10 20 cm −2 K −1 / (km s −1 ) in the following Galactocentric radial intervals: [0 − 3 kpc; 3 − 5 kpc; 5 − 6 kpc; 6 − 7 kpc; 7 − 15 kpc; 15 − 30 kpc].We assume here that the ISM gas is a mixture of Hydrogen and Helium nuclei with uniform density ratio f He = 0.1.

CR transport: the conventional and γ-optimized scenarios
We determine the energy and spatial distribution of each relevant CR species solving numerically the transport equation with the DRAGON2 code (Evoli et al. 2017(Evoli et al. , 2018)).We assume that the observed CR spectrum can be approximated as a steady-state solution for a smooth distribution of continuous sources, which we fix on the basis of SNR catalogues (here we use the SNR distribution reported in Ferriere ( 2001)).For a given source spectrum -generally a broken power-law tuned against locally measured CR spectra -as an output the code provides the propagated spectra of each primary and secondary species in each point of the Galaxy.Besides several astrophysical quantities, as an input the code needs to receive the CR diffusion coefficient D(ρ, x) as a function of the particle rigidity ρ and of the spatial coordinates.
In the conventional scenario this is assumed to be a single power law function of the particle rigidity with a spatially dependent slope, parameterized as follows: , where D 0 is its normalization at a reference rigidity ρ 0 = 4 GV3 , and β is the velocity of the particles in units of the speed of light.The index δ, a priori being poorly known, is inferred from the comparison with the measured secondary to primary CR flux ratios, the boron-to-carbon (B/C) ratio being the most common.Works based on multi-channel analysis (Génolini et  found that at the Solar System δ(R ) 0.5.A different scenario arises if δ = δ(x) which turns into a non-factorized dependence of the propagated CR spectra on energy and position.For the models studied here, the Alfvèn velocity is taken to be V A = 13 km s −1 , the normalization of the diffusion coefficient is D 0 = 6.1×10 28 cm 2 s −1 and the halo size is H = 6.7 kpc, in agreement with recent analyses of 10 Be ratios (De La Torre Luque et al. 2021).We checked that passing to the γ-optimized scenario has no effect on the local B/C (see e.g.Gaggero et al. (2015b)), as well as on other secondary-to-primary CR ratios, which indeed are correctly reproduced with this setup.We notice that adopting a different value of H and rescaling D 0 so to keep their ratiohence the CR grammage -unchanged would have no significant effect on the CR propagated spectra and the γ-ray diffuse emission along the GP.Concerning the vertical -with respect to the GPz dependence of the diffusion coefficient, we assume here a simple step function D(ρ, R, z) = D(ρ, R)Θ(H − |z|).Since the target gas is concentrated in a thin disk with z gas H, varying the value of H within errors, or using a smoother vertical profile, would not have significant effects on the line-of-sight integrated γ-ray emission.
Probing only the spectra of CR reaching the Solar System, the detection of charged CRs does not allow to discriminate among the factorized (conventional) and non-factorized scenarios.Indeed, the first evidence supporting the latter scenario was found in Gaggero et al. (2015b) on the basis of the Fermi-LAT results showing an excess of the diffuse γ-ray emission of the Galaxy above 10 GeV in the inner GP respect to the predictions of the conventional one.In that work it was shown that the Fermi-LAT results are reproduced if δ has a linear dependence on the Galactocentric radius R which turns into a harder spectrum of CR protons at low R, hence of secondary γ-rays at low Galactic longitudes (see also Acero et al. 2016).For this reason, in the following we will call this γ-optimized scenario.Afterwords, it was shown (Cerri et al. 2017) that this scenario is theoretically motivated and arises as a consequence of the growing poloidal component of the Galactic magnetic field at small Galactocentric radii (see e.g.Jansson & Farrar 2012).
In the following we will consider two transport setups: the Base one, which is representative of the conventional scenario, and the γ-optimized one.The main parameters of those models are reported in Table 1.For the γ-optimized setup we find δ(R) = 0.04(kpc −1 ) • R(kpc) + 0.17, for R < R = 8.5 kpc and δ(R) = δ(R ) for R ≥ R .As discussed at the end of Sec.2.1, for each model we will use two different choices -Min and Max setups -of the CR proton and Helium source spectra, which enter to determine the spectra propagated to the Earth, so to bracket the experimental uncertainties above 10 TeV (see Fig. 1).This arrangement effectively accounts also for the scenario in which the features in the propagated spectra are originated by transport rather than by the acceleration mechanism close to the sources.CR electron source spectra are also tuned to match local experimental data up to the few TeV's -see Both models are implemented with the DRAGON2 code (Evoli et al. 2017(Evoli et al. , 2018) ) -which was designed to enforce spatially dependent diffusion -using the same cross-sections and IS gas and ISRF as discussed below.Hadronic emission maps due to CR interactions on IS Hydrogen and Helium, as well the IC emission by electrons, are then computed with the HERMES code (Dundovic et al. 2021).We adopt the γ-ray production cross-section described in Kelner & Aharonian (2008) with the updated parameterization of the proton-proton total inelastic cross-section reported in Kafexhiu et al. (2014).For Helium, we assume that it contributes with a geometrical weight estimated considering its atomic radius as a function of the atomic mass A, R(A) A 1/3 R 0 , being R 0 the reference proton radius.Since its contribution appears due to the pp interaction cross-section, then its geometrical weight becomes R(A He ) R(A p ) A 1/3 He 2 = 4 2/3 , with respect to protons.At high energies the absorption by pairproduction on the spatially-independent CMB photon field is computed as shown in Dundovic et al. 2021.The effect of the infrared Galactic Background is negligible compared to the uncertainties and we neglect it in this work.In Fig. 2 we compare our models with the spectrum of the diffuse emission below 300 GeV obtained from Fermi-LAT data as described in Sec.2.2.Not to overcrowd the figure, for the Base transport model only the Max setting is shown.A more comprehensive comparison is presented in Fig. A.1 in Appendix A. As mentioned in the above, the γ-optimized setups provide a better agreement with Fermi-LAT results a low Galactic longitudes.
We notice that simplified versions of the Base and γoptimized models considered in this work -adopting an exponential cutoff at 5 PeV/n and neglecting absorption -where used in Acharyya et al. (2021).

LHAASO IC_86
3.24993e-21 2.86481e-16 J [cm 2 s 1 GeV 1 sr 1 ] Injection parameters  1: Spectral indexes at injection for the Max and Min models.These spectral indexes are tuned to CR local data as described above and correspond to spectral breaks at the following energies: 335 and 6 • 10 6 GeV for the Max models and 335, 2 • 10 4 and 4 • 10 6 GeV for the Min models.
We compute the full-sky maps of the diffuse gamma-ray emission associated to π 0 emission, Inverse Compton scattering and Bremsstrahlung with the HERMES code (Dundovic et al. 2021).We choose an angular resolution characterized by the Healpix resolution pararameter nside = 512, corresponding to a mean spacing between pixel of 0.11 • (Górski et al. 2005), nicely matching the angular resolution of the gas models adopted to compute the hadronic emission.For illustrative purpose, we show the Mollweide projection of the total emission associated to the γ-optimized Min model in Fig. 3, in a lower resolution.
In order to directly compare our models to the different experimental results described above, we consider several regions of interest, directly associated to the spectral data provided by the experiments focused on the very-high-energy domain.In particular, we show in the same Figure the contours of the regions observed by LHAASO (coincident with Tibet ASγ and ARGO) and IceCube-86.
We obtain the integrated flux in these regions, which we compare to the experimental data without any further ad-hoc tuning and post-processing.We emphasize once again that all the details of the setup (in particular, the ring-by-ring normalization of the molecular gas density, and the CR transport setup) are set by the comparison with both local data on charged CRs and Fermi-LAT data in the GeV-TeV domain, as commented in more details in the Appendix.The results are presented in Fig. s (4) and ( 6).The absorption due to γ − γ scattering is accounted as described at the end of Sec.3.2.Its effect is shown in Fig. 7 for the γ-optimized scenario.
Fig. 4, in particular, clearly represents the main result of this paper.This plot demonstrates that the diffuse emission models presented in this work -obtained under the assumption that the emission is fully originated by the diffuse Galactic CR "sea"are able to capture the main features of the observed data in a remarkably large range of energies, from 10 GeV all the way up to the PeV domain.This is already a major result.
However, since we are willing to go beyond this first level of interpretation and use our results to learn something about Galactic CR properties we face two main problems: there is a significant degeneracy between the choice of the CR transport setup and that of the source spectra (which, as we shown, depends also on the CR data systematics); there is a significant scatter of the Tibet and LHAASO data above 50 TeV.
While this situation is likely to improve with the next data releases we may already get some valuable hints limiting ourselves to consider only the lowest energy bin of both experiments which should be affected by lower systematics.Interestingly we notice that the four lowest energy LHAASO points -below 50 TeVare well aligned among themselves and the Tibet ones.We notice that those data favour the γ-optimized Max model.Even if we were to disregard Tibet data, or assume them to be contaminated by the emission of the Cygnus cocoon (see Sec. 2.3), the γ-optimized scenario would remain the preferred one though in its Min realization (see also Fig. 7).Although the Base -Max model is also in reasonable agreement with LHAASO data it is disfavored by Fermi-LAT and ARGO results.This shows the importance of using data over the widest possible energy range.(Amenomori et al. 2021) and LHAASO (Zhao et al. 2021) (preliminary) data in the window |b| < 5 • , 25 • < l < 100 • .The Galactic diffusion emission spectrum measured by Fermi-LAT and extracted as discussed in Sec.2.2, as well as ARGO-YBJ data (Bartoli et al. 2015) in the same region, are also reported.The models account for the effect of γ-ray absorption onto the CMB photons (see Sec. 3.2).
We also consider the Tibet ASγ data in the window |b| < 5 • , 50 • < l < 200 • (Fig. 5).We notice that in this more external region the predictions of the γ-optimized and Base scenarios are quite similar so that those data may help to remove the degeneracy between the choice of the transport scenario and the shape of the source spectrum.Remarkably, even accounting for a possible contamination due to Cygnus-OB2, Tibet results seems to neatly favour the Max setup for the latter unknown.It will be very interesting, therefore, to see if LHAASO will possibly confirm Tibet results in that region.This will be also relevant to scrutinize an alternative interpretation of Tibet results given in terms of the emission of unresolved pulsar wind nebulae (Vecchiotti et al. 2021).
We also performed a comparison of our models with Ice-Top and CASA-MIA upper limits which refer to regions different from those probed by Tibet and LHAASO (see Fig. 3).As evident from Fig. 6, where we also report ARGO-YBJ data, although those limits do not constrain any of our models yet, the IceTop sensitivity is close to the level required to test the γ-optimized Max model.Fermi-LAT systematic uncertainties dominate above ∼ 200 GeV, while the systematic error associated to TIBET data in this region is estimated to be around 30% (Amenomori et al. 2009).CASA-MIA (Borione et al. 1998) upper limits in the same region are also reported.

Conclusions
In this paper we presented a set of gamma-ray diffuse emission models that are able to consistently reproduce the available measurements from the GeV domain all the way up to PeV energies in the Galactic Plane region.
In particular, we discussed a reference model based on the assumption of homogeneous transport properties in the whole Galaxy, and an optimized model aimed at capturing the progressive hardening of the proton slope inferred from Fermi-LAT data in the GeV domain.Both scenarios are tuned on local CR data, and are presented in two different versions, that correspond to a different fitting strategy of local CR data in the very-high-energy part of the spectrum, which results in different choices of the injection spectra.
We found a relevant degeneracy between the choice of the CR propagation setup and that of the injection spectrum.In spite of that, we argued that the comparison between our models and the combination of different γ-ray data sets is able to provide valuable hints and may help to shed light on CR transport properties in different regions of the Galaxy.
We highlighted in particular that the Galactic diffuse emission measured by most experiments from 10 GeV to the PeV can almost entirely be explained as truly diffuse emission stemming from the Galactic CR "sea", and showed that the Tibet ASγ and LHAASO data, in combination with Fermi-LAT results, seem to favour a transport scenario characterized by spatially dependent diffusion.
Our results give a stronger support to a similar claim of the Tibet collaboration -based on the comparison of their data with the phenomenological models in Lipari & Vernetto (2018) which however did not take into account the degeneracy mentioned in the above.
The comparison of our models with forthcoming experimental results at lower energies by HAWC (Nayerhoda et al. 2020), H.E.S.S. (Abramowski et al. 2014), SWGO (Albert et al. 2019), CTA (Cherenkov Telescope Array Consortium et al. 2019;Acharya et al. 2018) and ALPACA (Takita et al. 2017;Takashi et al. 2021) will be crucial to further validate this scenario and to probe the CR shape throughout the Galaxy, especially by those experiments, which will have a better view of its most central regions.
In order to facilitate the comparison with these forthcoming data, we provide the scientific community with high resolution all-sky-maps of the diffuse γ-ray emission of the Galaxy from few GeVs to few PeVs computed for our benchmark models4 .They can be valuable tools for experimental collaborations and can be used as "background models" in different contexts, from the generation of Galactic and extra-Galactic source catalogues to indirect dark matter searches.
As a final discussion point, let us return to the potential role of unresolved sources.In general, the relative weight of truly diffuse emission and unresolved source contribution depends on a wide range of parameters, that characterize: the nature of the sources, the capability of the experiment to identify and resolve individual sources, the transport/escape of the high-energy particles that constitute the diffuse CR sea, and of course on the total amount of target gas and photon background that is directly responsible for the truly diffuse signal.In this work, we choose to consider a specific model for the unresolved contribution and take it into account consistently when tuning our prediction for the diffuse emission on the Fermi-LAT data.Due to the uncertainty associated in the modeling of this component, and the impossibility to directly extrapolate the model to the PeV domain (a task that would require an accurate characterization of each PeV-oriented experiment to resolve individual sources), we cannot exclude a degeneracy between our scenario and an alternative model featuring a more prominent role of the unresolved component.However, in order to safely disentangle different interpretations, to have more data will be crucial.In particular, we emphasize the role of neutrino measurements in this context (see also Ahlers et al. 2016;Gaggero et al. 2015a;?;Breuhaus et al. 2022).In fact, the truly diffuse component is likely to be entirely hadronic above few TeV, while the sources that could potentially contribute to the unresolved flux are mostly leptonic (namely, pulsar wind nebulae).Hence, a future firm discovery of a neutrino diffuse component associated to the Galactic plane would be a clear signature of the scenario we propose here.

Fig. 1 :
Fig.1: The proton (left panel) and Helium (right panel) local spectra computed for the γ-optimized scenario are plotted against a representative set of data.For each species the spectra as predicted using the Max and Min source spectrum set-ups are shown.We do not show here the corresponding lines computed for the Base scenario since they are almost coincident with those reported here above 10 GeV/n (at the Solar System position).
Fig. B.1 in the Appendix B.

Fig. 2 :
Fig. 2: The average spectrum of the γ-ray diffuse emission along the galactic plane (|b| < 5 • ) is compared with Fermi-LAT data in three longitude intervals.The "Sources" component comprises the sources included in the 4FGL Fermi Catalog as well as unresolved sources and the Interstellar Galactic background light ("IGB") component comprises the extra-galactic background light and Fermi's instrumental background.The errorbars represent just the statistical error of the measurements.

Fig. 3 :
Fig. 3: The Mollweide projection of the all sky map of the Galactic diffuse emission flux above 100 TeV obtained for the Min γ-optimized model is reported in this figure.The contours of the regions probed by Tibet/LHAASO and by IceCube (IC-86) are reported.J, in the legend, means the differential flux of γ-rays per unit of solid angle.

Fig. 4 :
Fig.4: The γ-ray spectra computed within the conventional (base) and γ-optimized scenarios are compared to Tibet ASγ(Amenomori et al. 2021) and LHAASO(Zhao et al. 2021) (preliminary) data in the window |b| < 5 • , 25 • < l < 100 • .The Galactic diffusion emission spectrum measured by Fermi-LAT and extracted as discussed in Sec.2.2, as well as ARGO-YBJ data(Bartoli et al. 2015) in the same region, are also reported.The models account for the effect of γ-ray absorption onto the CMB photons (see Sec. 3.2).

Fig. 5 :
Fig.5: Predicted γ-ray spectra for the different scenarios studied in this work and compared to Tibet ASγ(Amenomori et al. 2021) and Fermi-LAT data in the window |b| < 5 • , 50 • < l < 200 • .The experimental errorbars show the 1σ statistical uncertainty of the measurement.Fermi-LAT systematic uncertainties dominate above ∼ 200 GeV, while the systematic error associated to TIBET data in this region is estimated to be around 30%(Amenomori et al. 2009).CASA-MIA(Borione et al. 1998) upper limits in the same region are also reported.

Fig. 6 :
Fig. 6: The γ-ray spectra computed in the conventional (base) and γoptimized CR transport scenarios are compared to IceCube (Aartsen et al. 2019b) and CASA-MIA (Borione et al. 1998) upper limits.Since those data refer to different regions of the sky, they are rescaled as described in Aartsen et al. (2019b) (see Fig. 12 in that paper).

Fig. 7 :
Fig. 7: In this figure we show the effect of γ-ray absorption onto the CMB photons (see Sec. 3.2) for the γ-optimized scenario.