{'fraction_non_alphanumeric': 0.029093795969832116, 'fraction_numerical': 0.0155976370046933, 'mean_word_length': 4.590392899099334, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the progressive hardening of the cosmic-ray proton spectrum in the inner Galaxy

Spatial variations of the average properties that characterize the hadronic component of the diffuse Galactic cosmic-ray sea, in particular the spectral slope and normalization, may unveil critical information about their confinement mechanism in the Galaxy. In the first part of this paper we perform an analysis of the Fermi-LAT gamma-ray data with the SkyFACT package, which combines image reconstruction techniques with standard template fitting, isolate the hadronic emission and decompose it into Galactocentric rings. We find a significant hardening of the hadronic spectral index towards the molecular ring. We study this hardening in different energy ranges, and assess its resilience with respect to different prescriptions in the analysis setup. In the second part we quantify the contribution to the diffuse gamma-ray flux coming from unresolved point sources with a dedicated Monte Carlo simulation, and consider whether the trend characterized in the first part can be mimicked by a progressively more relevant flux associated to this component in the inner Galaxy. We find that the observed hardening of the hadronic spectral index cannot be due to unresolved sources in the sub-TeV energy range, especially outside the molecular ring, given reasonable assumptions about the unresolved source population.


Introduction
The study of the gamma-ray emission from the Galactic plane provides a unique opportunity to understand the transport properties of cosmic rays (CRs) in different regions of the Galaxy. The basic picture for the production and propagation of Galactic CRs, which dates back to pioneering studies in the early 1960s [1], is based on the presence of a diffuse sea of high-energy particles confined in a large diffusive halo. The motion of those particles is described by a random walk, relatively easily modeled in terms of a homogeneous, isotropic diffusion equation, with the addition of advective and loss terms. Although the actual physical processes characterizing the interaction between the CR sea and the magnetized, turbulent interstellar plasma, which are supposed to be responsible for their random walk, and ultimately for their confinement, are far from being understood, the phenomenological picture stated above was considered adequate to reproduce the available data for a long time.
Currently, given the high statistics of the most recent CR measurements, and the sharper view of the gamma-ray sky in the multi-GeV domain provided by the Fermi-LAT instrument, it is possible to challenge this standard picture. It is therefore important to identify clear trends in the data, for example associated to spectral changes across the Galaxy, since these features are likely to provide useful insights on the physics of CR transport across the Galaxy.
A remarkable example of a large-scale trend with possible relevant implications for the microphysics of CR transport is the progressive hardening towards the inner Galaxy of the proton slope inferred by gamma-ray data. This feature was first outlined in [2], and later confirmed by two different model-independent analyses [3,4].
The first question that naturally arises from those results is whether the trend is actually related to the diffuse component of the gamma-ray sky, and can therefore be the signature of some anomalous or non-standard CR transport process at work, or instead can be ascribed to a progressively larger number of unresolved sources in the inner regions of the Galaxy.
Once this aspect is clarified, it is important to understand possible physical scenarios able to reproduce such a hardening. In this context, two different explanations were proposed in the literature.
The scenario presented in [5] is based on two pieces of information: the presence of a poloidal component of the regular Galactic magnetic field which is directed in the perpendicular direction with respect to the Galactic plane in the inner Galaxy; and the harder scaling of the parallel (w.r.t. the regular field) diffusion coefficient with rigidity that can be inferred from numerical simulations [6,7]. Given these inputs, the authors build a model based on anisotropic CR transport featuring a progressive transition from perpendicular to parallel vertical escape from the Galaxy, and demonstrate that the picture is compatible with the data.
In [8], instead, the role of CR self-confinement via the streaming instability is taken into account. GeV CRs are expected to trigger the growth of the Alfvén waves that are in turn responsible for their own confinement, and this effect is expected to be stronger where stronger CR gradients exist, i.e. in the inner Galaxy. As a consequence, the diffusion coefficient in the inner regions of the Galaxy should be much lower (reflecting the more effective confinement), so advection should take over up to large energies and the propagated spectrum should be therefore closer to the injected one, explaining the hardening.
The two interpretations predict different behavior at high energies: In fact, the former is expected to hold at all energies, while the latter is based on a physical effect that becomes progressively less relevant with increasing energy. The high-energy CR flux steeply declines with rigidity. Therefore, at energies 100 GeV, their ability to trigger the streaming instability should be negligible [9,10] 1 .
Given this important distinction, a natural way to disentangle the two interpretations is to analyze the high-energy gamma-ray data (above 30 GeV), and identify whether the same trend pointing towards a hardening in the inner Galaxy is present. If so, an interpretation in terms of transport properties, such as anisotropic diffusion, that hold at all energies would be favored. This paper addresses the two important questions outlined above. In the first part we adopt, for the first time in this context, the SkyFACT package [11], a flexible statistical tool designed to model and decompose the gamma-ray sky into its various components. It allows to account for a very large number (∼ 10 5 ) of nuisance parameters in order to marginalize over the imperfections of current models built with CR transport codes. We use this analysis tool to robustly characterize the dependence of the proton spectral index as a function of the Galactocentric radius in different energy ranges, with the ultimate goal to shed light on the origin of this anomaly. In particular, we perform for the first time an energy-dependent analysis in order to assess the presence of such a trend in Fermi-LAT data above 30 GeV, and discuss different systematic effects that can affect the result, including the degeneracy with Inverse Compton emission.
In the second part of the paper, we quantify the relevance of unresolved point sources (UPS) via a population study based on currently available resolved point source catalogs and our current knowledge on the luminosity functions of the most relevant classes of sources, and discuss whether a significant contribution from UPS can mimic the progressive hardening.
2 Hadronic gamma-ray slope analysis: Method The gamma-ray emission from the Galaxy is due to a variety of mechanisms. CR nuclei and leptons, during their random walk through the interstellar medium (ISM), interact with interstellar gas, magnetic fields and the diffuse low-energy interstellar radiation field (ISRF) emitted by stars and further reprocessed by dust grains. As a result, hadrons produce neutral pions which in turn decay into energetic gamma rays; leptons, on the other hand, are able to up-scatter the ISRF photons, and emit bremsstrahlung radiation as well due to the interaction with the electric fields associated to the ionized component of the ISM.
In order to gain information on the hadronic CR component and its spectral properties, it is therefore crucial to disentangle these different types of emission, associated to different targets and characterized by different morphologies. To this aim, we perform a novel analysis of the gamma-ray emission from the Galactic plane by exploiting the recently developed SkyFACT tool [11].

SkyFACT
SkyFACT uses a combination of template fitting and image reconstruction methods to decompose the gamma-ray sky into its various components. It can account for the (still numerous) imperfections of diffuse emission models by giving the individual model components some flexibility to better match the data. We give a brief overview of the procedure of fitting the gamma-ray sky using SkyFACT but refer the reader to [11] for a detailed description.
Traditional analysis of diffuse gamma-ray emission relies on fitting gamma-ray data with a linear combination of templates that capture the different emission processes mentioned above. These templates, however, come with significant systematic uncertainties that are often difficult to quantify, which create relevant biases when interpreting the gamma-ray sky. SkyFACT moves beyond the "static" template fitting paradigm and allows the templates to vary in both morphology and spectrum. Even if the deviation from the original templates remain small, this requires a large number of extra parameters (as many as the number of pixels plus the number of energy bins), which are treated as regularized nuisance parameters. Regularization terms in the model likelihood are then implemented such that deviations from the original templates get penalized in a way that reflects the expected systematic uncertainties. SkyFACT also allows for the possibility to leave complete spectral or spatial freedom for some of the templates, in which case it effectively behaves like an image reconstruction tool for components that cannot be modeled a priori.
In more detail, the diffuse gamma-ray photon flux in pixel p and energy bin b is modeled as: where the sum is over k model emission components and T and S are, respectively, the spatial and spectral templates. The set of parameters τ (k) and σ (k) are spatial and spectral nuisance parameters, respectively, and ν (k) is the overall normalization parameter, for model component k. These parameters essentially represent the normalizations per spatial pixel, per energy bin, and overall, for each component. The model flux in Equation 2.1 is then convolved with the point spread function and exposure to give the expected number of photons in each bin.
As mentioned above, in order to constrain the variability of the nuisance parameters τ , σ, and ν, additional regularization terms are added to the likelihood. These terms are controlled by three hyper-parameters (λ, λ and λ ) that constrain the spatial, spectral and overall modulation. We adopt the maximum entropy method (MEM) as a regularization algorithm, a technique commonly used in the context of astronomical image reconstruction [12]. This method is designed in order not to introduce any spurious structures or patterns in the data. Hence, the reconstructed images we get appear as featureless as possible, while still being consistent with the input data.
SkyFACT adopts the L-BFGS-B [13] algorithm to perform maximum likelihood estimation. The final likelihood takes the following form: where L P is a standard Poisson likelihood and L R represents the set of regularization terms (MEM) that constrain the nuisance parameters τ (k) , σ (k) , and ν (k) , with the beforehand set regularization hyper-parameters. We also include some additional smoothing conditions on the spatial and spectral modulation parameters that control the pixel-to-pixel and bin-by-bin variation; this introduces an additional pair of hyper-parameters per component (see [11] for further details). Point sources have no spatial modulation or smoothing, so there are 3 more hyper-parameters that control the point source spectral modulation, spectral smoothing, and overall normalization.
In this way, the strength of the allowed modulation can be tuned for each model component, effectively allowing us to include the uncertainties intrinsic to the model. A table of the hyper-parameter values for each model component is shown in Table 1.
The statistical uncertainties on the modulation parameters are estimated from the Fisher information matrix, the inverse of which is a lower limit on the covariance [14]. It is computationally difficult to directly perform this inversion, so in practice we use Cholesky decomposition to efficiently sample from a multivariate Gaussian such that the covariance of the sampled model parameters is equivalent to the inverse of the Fisher matrix (see [11] for more details). This procedure yields uncertainties that should be interpreted as conservative estimates on the true statistical uncertainties.

Fermi-LAT data
For this study we use 9.3 years, from 4 August 2008 to 6 November 2017, of Fermi-LAT Pass 8 data. We select only ULTRACLEAN events (evclass=512) and FRONT and BACK (evtype=3) events. A zenith angle cut at 90 • was used to reduce photons from the Earth's limb and the data quality cuts recommended by the Fermi-LAT collaboration (DATA_QUAL>0)&&(LAT_ CONFIG==1) were applied. We use standard LAT data analysis tools as available from the NASA website 2 to create counts and exposure maps of 720 × 81 pixels with a pixel size of 0.5 degrees centered at the Galactic center. The photon energy was binned in 25 logspaced energy bins between 0.34-228.65 GeV and our region of interest spans the whole disk (|l| < 180 • ) with latitudes |b| < 20.25 • . The model components and regularization parameters used in this study are very similar to the settings of the model labeled RUN4 in [11] (see Table 1 therein). We will summarize the most important components here and note important differences between the analysis in this study and the one in [11].

Model components
We use the gas column densities (HI+CO) as available from the GALPROP public release to build the π 0 template and bin them in 9 Galactocentric annuli. This binning facilitates comparison with previous work [3]. The spin temperature for the HI maps is fixed to T S = 125 K and the molecular hydrogen proportionality constant X CO = 1.9 × 10 20 cm −2 /(K km/s).
The initial spectrum for the π 0 templates is taken from Ref. [15]. This spectrum is in agreement with local cosmic-ray observables and approximately follows a power-law scaling with an index of 2.7 at energies > 2 GeV. The spectral regularization parameter (λ ) of the gas rings is chosen so that the spectra are constrained within 25% of the initial spectrum.
The spatial regularization is set so that variations are within ∼ 32% of the initial template. The best-fit gas templates are shown in Figure 1.
Note that we do not implement any correction to take into account dark gas, i.e. gas that is not resolved by the 21cm or CO line emission, because it was shown that this can be recovered through the use of the modulation parameters [11].
The ICS model we used was created with DRAGON and GammaSky [16][17][18]. We use the ISRF as documented in [19] with the Ferrière source model [20] and the 'KRA4' cosmic-ray propagation model from [21]. The initial ICS spectrum is again taken from [15] and allowed to vary by 25%. The ICS spatial template is constrained to within a factor two.
A notable difference between this work and [11] is that we do not add a template for the GeV excess in the Galactic center and use a more conservative template for the Fermi bubbles, following the region as defined in [22], and allow no spatial variation for this bubble template. The spectrum of the bubbles is taken from [23] and is constrained to within 1%. We explore the effects of these choices on the final results in section 4. Furthermore, we add the necessary point sources and extended sources within our ROI from the 3FGL catalog [24] and treat them as described in [11].
9× [ 10 16 0 25 0 · ] Table 1. An overview of the model components and the regularization hyper-parameters as used in the gamma-ray sky fit with SkyFACT . The second column contains a description of the modeling uncertainties that we allow in the fit. The matrix contains the corresponding regularizing hyperparameters, where λ, λ and λ are spatial, spectral and overall modulation parameters respectively and η and η the spatial and spectral smoothing.
3 Hadronic gamma-ray slope analysis: Results

Hardening of the proton spectrum in the full energy range
In Figure 2, the best-fit spectra for each Galactocentric annulus from SkyFACT are shown with black points. In order to identify any trends in the spectra across different rings, we then fit the black points with a single power law in the 2 − 228.65 GeV energy range. We are ultimately interested, however, in trends present in the CR proton spectrum. Above a gamma-ray energy E γ ∼ 2 GeV (the location of the π 0 bump), the spectrum of the gammaray emission from π 0 decay scales with the underlying CR proton spectrum, with an overall 10 0

− 50 kpc
Best fit index: 2.73 2.7+UPS Figure 2. Hadronic gamma-ray spectra in different gas rings. The black points represent the best-fit spectra to each gas ring template recovered from SkyFACT . The red line corresponds to the best-fit power-law from 2-228.65 GeV, and the pink band is the 68% credible interval. The gray dotted line shows a power-law with an index of 2.7 for comparison. The blue dashed line shows the contribution of unresolved point source for our reference scenario, added to the reference power-law (see section 5). The first radial bin is grayed out because the fit in that bin is strongly affected by systematic uncertainties (see subsection 3.1).
factor of E γ ≈ E CR /10. Therefore, the photon spectral index of the π 0 -induced gamma-ray emission for E γ > 2 GeV is directly comparable to the CR proton spectral index.
To perform a power-law fit to the black points above 2 GeV in Figure 2, we assume the photon flux in each energy bin is normally distributed with a standard deviation as sampled from the inverse Fisher matrix and set up a Markov Chain Monte Carlo (MCMC) using the emcee package [25] to sample the posterior distribution of the power-law index. The mean and 68% credible intervals of the best-fit power-law is shown in Figure 2 for every Galactocentric annulus.
A trend of hardening spectra towards the inner Galaxy is already apparent in Figure 2. The black points and the best-fit power-law index (the red line) are clearly harder than the reference power-law with index of 2.7, which is the observed slope of the local CR emission  We show the spectral index fitted as explained in the text from 2-220 GeV compared to the trend found in [3,4]. Horizontal error bars indicate bin width in R and vertical error bars are 68% credible intervals.
The best-fit photon index from Figure 2 is plotted versus radial distance from the Galactic center in Figure 4 and compared with previous analyses [3,4]. We also show the hadronic gamma-ray emissivity integrated over energies above 1 GeV (for straightforward comparison with [4]; [3] shows the emissivity at 2 GeV), which is a proxy of the hadronic CR flux, in Figure 3. Overall, we find a reasonable agreement with both studies. Figure 3, the gamma-ray emissivity decreases as a function of the Galactocentric radius in the region beyond the molecular ring, which is located in the second ring corresponding to the 1.7 − 4.5 kpc range. The decline in emissivity is not as steep as naively predicted in the context of isotropic and homogeneous diffusive CR transport, as already pointed out several times even before Fermi-LAT data (the so-called gradient problem: see e.g. [28] and the more recent discussion in [8,29] Spectral index π 0 2-30 GeV π 0 30-220 GeV π 0 2-220 GeV Figure 5. Energy-dependent analysis of the hadronic gamma-ray slope. We show the spectral index of gamma rays associated with π 0 decay at different radii from the Galactic center. The power-law fit is first performed in the full range, and then restricted to both the low-energy (from 2-30 GeV), and the high-energy (30-220 GeV) domain; the low-energy and high-energy points are artificially shifted to the right by 0.25 kpc and 0.5 kpc respectively, so that the error bars do not overlap. Figure 4 clearly shows the progressive hardening towards the Galatic center in a wide range of radii outside the Galactic bulge. This trend is qualitatively compatible with previous studies. The uncertainties on the power-law index in this work are smaller than previous analyses due to a combination of different effects, in particular more statistics with respect to previous papers, and the different technique implemented in SkyFACT . We discuss these issues in more detail in section 4.

As shown in
The first radial bin from 0-1.7 kpc, which contains the Galactic bulge, is quite problematic, and different studies find different results. This is not too surprising, given the small number of pixels associated to this first ring, which results in a low normalization of the hadronic component associated to this bin as a result of our fitting procedure due to lack of constraining power. More importantly, there is large degeneracy between the different components in the Galactic center region, which is by far the most challenging issue both from the observational and modeling point of view. However, a lower-than-average CR flux in that region may be compatible with a scenario characterized by a very fast diffusive escape in the vertical direction, parallel with respect to the direction of the regular magnetic field, if the diffusion tensor is highly anisotropic (see Fig. 4 in [5]).

Energy dependence of the hardening of the proton spectrum
We want to assess whether the trend in the spectral hardening is still present at higher energies. This would help to discriminate between the two theoretical interpretations discussed above. In particular, the streaming instability, mostly taking place in the hot ionized phase of the interstellar medium, is not expected to be the dominant confinement mechanism for CR hadrons more energetic than ∼ 100 GeV, because the CR flux, and therefore the associated pressure, is not large enough to effectively trigger the growth of Alfvénic perturbations. As mentioned in the Introduction, this was already noticed in [9], building on previous works (e.g. [30,31]), and recently discussed in [10,32] in the context of the local spectral feature in CR nuclei detected by PAMELA and AMS. As a consequence, higher energy CRs should be confined by the interaction with pre-existing turbulence, and the argument presented in [8] cannot be invoked. The model based on anisotropic transport, on the other hand, should hold at all energies [5]. The assessment of the trend for gamma-ray energies above few tens of GeV (which trace CR hadrons more energetic than ∼ 100 GeV) is therefore crucial in this context.
In order to search for the presence of a hardening trend at high energies, we model the gamma-ray sky as described above, but due to the increasingly lower photon counts at higher energies, we use a smaller number of gas rings: We merge rings 3 and 4 (4.5-6.5 kpc), and rings 5 and 6 (6.5-8 kpc) and therefore adopt a total of 7 rings instead of 9. We model the gamma-ray data with SkyFACT following the method outlined above; however, we now fit the π 0 emission with a broken power-law featuring a break at E γ = 30 GeV.
The spectral index from 2-30 GeV and from 30-220 GeV are shown in Figure 5 along with our result over all energies from Figure 4. The plot shows that the trend is clearly present in the high-energy domain, in particular outside of the Galactic bulge, and is even more pronounced than the one inferred by low-energy gamma rays. The reason why the highenergy spectral slopes are systematically harder at all radii, with a more pronounced trend in the molecular ring region, may be due to a combination of different effects. The spectral feature in the CR primary spectrum observed by PAMELA and AMS could be present all through the Galaxy (in particular if interpreted as the transition between the scattering on self-generated and pre-existing turbulence), and is expected to be imprinted in the gammaray slope as well. Additionally, unresolved sources are expected to contribute more at larger energies, altering the trend. This effect will be quantified in section 5.

Systematic uncertainties
Let us now consider the systematic uncertainty associated to our analysis. To this aim, we present several SkyFACT analyses performed with different choices in components and modulation parameters. We refer to Figure 6 for a visualization of these results.
First, we consider a fit to the gamma-ray data where the ICS component is split over different radii from the Galactic center: 0-3 kpc, 3-8.3 kpc, 8.3-50 kpc. This splitting is similar to the one applied to the hadronic component, but restricted to three rings only in order to keep the number of fitted parameters under control. We also restrict the ICS spectrum to be more constrained towards the input spectrum. The idea here is that the harder ICS spectrum can pick up some of the residual high energy photons in the inner Galaxy that were previously being absorbed by the gas templates. We see however (blue circles in Figure 6) that there is still a clear trend of a hardening towards the inner Galaxy.    Table 2.
We also check whether adding a template for the Galactic center excess and allowing more freedom on the Fermi Bubbles affects the trend in the spectral index. In Figure 6 we show (red crosses) the effect of including a Galactic center excess template that follows Red Cluster Giants and the Nuclear Bulge (RCG NB), which according to Ref. [33] gives the best fit to the data. The addition of these components slightly affects only the first radial bin, but still within the error bars of the baseline model.
Then, we perform the analysis with opposite prescriptions regarding the freedom on the different templates: we do not allow any freedom on the spatial modulation and, conversely, Model Comment Baseline model As described in section 2. Split ICS ICS template split over three rings: 0-3, 3-8.3 and 8.3-50 kpc. Free Bubbles + RCG NB Spatial modulation of Fermi Bubbles free to vary with template for the Galactic Center excess as in the best fit from [33].

No spatial modulation on Gas+ICS
Similar to the traditional analysis as in [3,4]. The spectral modulation is more free: λ = 1. Free spectrum on Gas Spectral modulation on each gas ring is more free: λ = 4. Free spatial + free spectrum on Gas Spectral modulation as above, with more free spatial regularization λ = 4. 4 years of data Same analysis as baseline model but now only with data from the first 4 years of Fermi-LAT. Table 2. Overview of the several systematics of the gamma-ray fit that were tested in this paper. more freedom in the spectrum. This setup essentially reproduces the usual rigid template fitting technique adopted in [3,4]. With this prescription, we obtain larger residuals than in the baseline model but the trend in the slope is not altered very much (green triangles), except for the increase of the estimated uncertainties. We also show that allowing more freedom in the spatial and spectral regularization still qualitatively reproduces the same trend (small purple crosses) -there seems to be some degeneracy between radial bins when allowing more spatial freedom -with slightly larger error bars, as expected. The size of the error bars is mostly affected by the spectral regularization so that increasing the spectral variation from 25% to 50% leads to a twofold increase of the error bar on the spectral index (orange pentagons). Finally, we consider again the baseline model, but only using 4 years of data, in order to check the impact on the error bars. As expected, the uncertainty increases by the square root of the difference in observation time (cyan diamonds).

Unresolved sources
In order to draw more robust conclusions about the nature of the hardening trend outlined in the previous sections, it is important to quantitatively assess the role of unresolved point sources along the Galactic plane. Indeed, a significant population of unresolved hard gammaray sources that lie below the detection limit of Fermi-LAT could in principle mimic the progressive hardening towards the inner Galaxy. This would naturally explain why the radial dependence of the proton slope is very similar to the source distribution profile (which peaks at 5 kpc), and why the inferred CR spectrum is even harder at higher energies. Therefore, in this section we model the expected flux from a population of unresolved hard gamma-ray sources and quantify its relevance with respect to the diffuse π 0 emission via a dedicated Monte Carlo simulation.

Gamma-ray source catalogs
Population studies of gamma-ray sources in the Galaxy -mainly pulsars, pulsars wind nebulae (PWNe), and shell-type supernova remnants (SNRs) -are available from the 3FGL [24] and 1FHL [34] catalogs. In particular, the latest reference model from the 3FGL catalog, based on the spatial distribution of pulsars and on their spectra in the 1-100 GeV range, features 2.9 × 10 4 Galactic sources drawn from a power-law luminosity function characterized by a slope of −1.8 and a luminosity range 2 × 10 34 -2 × 10 39 ph s −1 . The conclusion of that study is that the unresolved sources in the inner Galaxy below 1 × 10 −9 ph cm −2 s −1 amount to 3% of the total diffuse emission at 1 GeV.
In this section we aim to characterize the relevance of unresolved gamma-ray sources with particular focus on the contribution at higher energies. Given that the average spectrum of point sources is harder than the diffuse emission, the contribution of point sources is expected to increase with energy and can become much larger than the value of 3% quoted above.

Estimating the contribution from unresolved point sources
It is useful to make a back-of-the-envelope estimate of the amount of flux from unresolved sources that one needs to induce a harder observed spectrum. Let us assume that the true π 0 spectrum has a declining slope with index γ, while the measured π 0 spectrum has a slope with index α and is contaminated by unresolved sources: where Φ UPS (E) is the flux from unresolved point sources, and the '∼' indicates that the same measured spectrum can be fitted with different model spectra. Furthermore, Φ 0 and Φ 1 are free normalization factors, which can in principle be obtained from fits to the data. The fraction f of the flux from unresolved point sources with respect to the measured π 0 spectrum can then be written as follows: This estimate holds assuming the true and the measured π 0 spectrum both give a good fit to the gamma-ray data at E 0 = 1 GeV so that Φ 0 Φ 1 (which seems to be the case, see [15]). Setting α = 2.5 and γ = 2.7 (see Figure 4) we find that, at 10 (100) GeV, 40 (60)% of the measured flux would need to be from unresolved sources to induce a hardening from 2.7 to 2.5.
We can also use the equations above to estimate the contribution of unresolved point sources to the additional hardening at E γ > 30 GeV on top of the already-present hardening at E γ < 30 GeV -that is, we compare the purple points to the green points in Figure 5. In this case, we set E 0 = E break = 30 GeV, and assume that that Φ 0 Φ 1 , (where Φ 0 and Φ 1 are now normalizations at 30 GeV). Then, at E = 100 GeV, unresolved sources would need to be 20% of the measured flux to induce a change in the power-law index from 2.5 to 2.3.
Following the same train of thought we can estimate the fraction of the integrated flux above the break energy E b one needs to induce a break in the spectrum at E b with index ∆γ. Let us assume that the diffuse flux is described by a single power law. The integrated flux above E b from a single power law simply reads: We assume that the addition of the UPS component induces a break: The integrated flux for a broken power law can be approximated by: Hence, the integrated flux from UPS above E b can be estimated as Ψ UPS ≈ Ψ 1 − Ψ 0 , and the ratio of this flux with respect to the initial single power law can be written as: Again, under the assumption that both the single and broken power law give a good fit to the data at ≤ E b , as seen in Figure 5, we get an expression for the fraction F as a function of the spectral index change ∆γ: .
This implies F 13% for a spectral hardening from a slope of 2.5 to 2.3 above 30 GeV. In the following section, we describe a more quantitative analysis of the estimated fractional contribution of unresolved point sources with a dedicated simulation.

Our simulation of Galactic sources
Our approach to simulate Galactic sources is based on [24,34] (which in turn refer to [35]). The code associated with this simulation can be found here 3 .
The simulation procedure is as follows: We produce ensembles of gamma-ray point sources by drawing the positions from a given spatial distribution in R, φ, and z. For each source, an integrated luminosity (in the 1-100 GeV range) is drawn from a luminosity function. Then, given distances and luminosities, we calculate a flux for each source. We then assign a random spectrum (in the 1-200 GeV domain) drawn from fitted spectra in the 3FGL catalog from the following classes of sources: associated or identified pulsars (PSRs), supernova remnants (SNRs), pulsar wind nebulae (PWNe), or unidentified sources with latitudes |b| < 5 • . The spectra of all considered sources is shown in Figure 7. Finally, we tune our result to the actual number of resolved point sources in the catalog.
Let us now describe the different astrophysical inputs mentioned above in more detail. The luminosity function (which refers to the 1-100 GeV range) is a power-law distribution of the form:  Table 3 are compared with the flux from resolved sources from the 3FGL catalog provided by the Fermi-LAT collaboration. The colored bands indicate the estimator of the standard deviation computed on an ensemble of 100 simulations.
which features a low-luminosity cutoff L min . The integral from L min to L max gives the total number of sources N . This simple power-law form for the luminosity function is widely used [34][35][36][37]; in particular, an index a of 1.5 was derived for pulsars in Ref. [35]. Let us remark that the true form of the luminosity function is uncertain, and the estimated flux from unresolved sources is highly dependent on the exact shape of this ingredient. However, since our only goal is to ascertain whether unresolved sources could potentially mimic a harder spectrum of the diffuse emission and we do not aim at finding the most accurate luminosity function, we adopt the reference model from [24] and define two additional scenarios that provide a reasonable bracketing of the uncertainties. A lower limit model can be obtained by choosing a larger value for L min , and a smaller (i.e., harder) value of the slope a, and vice versa for an upper limit model. The parameter L max has a minor impact on our results. We adopt the value 2 × 10 39 ph s −1 .
The exact spatial distribution of sources has a percent-level impact on our estimate up to energies of several GeV. Therefore, we only report here the results obtained from the distribution of pulsars taken from [38]. However, we tested other options [20,39] as well. The R dependence of the Lorimer type C source distribution can be written as:  Table 3. Different scenarios considered in this work for the parameters associated to the gamma-ray source luminosity function.
The threshold for unresolved sources is set at 1×10 −9 ph cm −2 s −1 . Given this threshold, along with a slope a and L min for the luminosity function, the total number of sources can be determined by requiring that the number of sources above this threshold match what is observed in the 3FGL catalog.
We consider as a reference a simulation based on the baseline model described in the 3FGL paper [24], and label it model 1 in Table 3. For this model, the slope and minimum luminosity are set to the following values: a = 1.8 and L min = 2 × 10 34 ph s −1 . Tuning the number of sources so that the number of the ones above threshold correspond to the catalog yields 29000 sources for this model.
In addition, we test two different scenarios that provide estimates of the upper and lower limits of our predictions. Model 2 incorporates a harder luminosity function with index a = 1.5 and a higher minimum luminosity L min = 4 × 10 34 ph s −1 . We find 3000 sources for this model. Model 3 has a softer luminosity function a = 2.1 and a lower value for L min = 1 × 10 34 ph s −1 . This yields 210000 sources. Note that with a supernova rate in our Galaxy of about 2 per century [40], and given that these sources are active in gamma rays for about 1 Myr, we expect only 2 × 10 4 of these sources in the Galaxy. Model 3 with 2.1 × 10 5 sources can thus be considered as an extreme upper limit, under the assumption that most of these sources result as byproducts of supernova explosions.
In order to estimate statistical fluctuations we produce 100 simulations for each model and report the mean and standard deviation in the following section.

Results: Contribution of unresolved sources at different energies
In Figure 8 we show the number of sources per flux interval for the three reference models described above, compared with all the identified Galactic sources from the 3FGL catalog (orange circles), and with all the identified sources (including the extra-Galactic ones) with |b| < 5 (black boxes): All the models considered are compatible with those data.
In Figure 9 and Figure 10 we show the longitude and latitude profile of the flux from unresolved sources compared to the diffuse π 0 emission outlined by our analysis (section 3): From these plots we can see that our model 1 reproduces the few percent contribution of unresolved sources to the diffuse emission found in [24]. We remark that the fraction from unresolved sources can be an order of magnitude larger in the inner Galaxy compared to the outer plane. Figure 11 and Figure 12 show the fraction of the flux associated to unresolved sources to the diffuse pion emission as a function of energy, averaged over the full ROI (|l| < 180 and |b| < 20) and the inner Galaxy (|l| < 40 and |b| < 5), respectively. The second region in Figure 12 corresponds to the region where the hardening is particularly pronounced (namely the R = 4.5-5.5 kpc range). In both plots, the dashed line represents the fraction that is   needed to mimic the hardening by inducing a change of the power-law associated to the diffuse emission from 2.7 to 2.5, with normalization set at 1 GeV, as estimated in subsection 5.2 (Equation 5.2). The dotted line is similar, but is normalized at 30 GeV and should be interpreted as the fraction needed to induce the additional hardening from low to high energies present in Figure 5. Over the entire Galactic plane, all the models considered predict fractions well below the "warning" lines in the whole energy range considered (1-100 GeV). In the inner Galaxy, instead, model 3 clearly overshoots the dotted line. However, we emphasize that such a scenario largely overproduces the total number of Galactic gamma-ray point sources sources expected in the Galaxy given our current knowledge on the SN rate, and can to be considered as a rather extreme upper limit on our estimate.
Given these results, we can conclude that the overall trend outlined with the SkyFACT -based analysis cannot be mimicked by a population of unresolved point sources with increasing density towards the inner Galaxy, especially beyond the molecular ring (R 4.5 kpc). The result is substantially unchanged if the tuning is performed on the 3FHL catalog.
A useful visualization is provided in Figure 2, where we show the UPS flux associated to model 1 added to the reference power law (compatible with the local slope). For R > 4.5 kpc, according to our reference model, the impact of the unresolved sources is clearly not enough to account for the harder gamma-ray emission. However, the UPS component can be responsible for the hardening in the second radial bin.
Finally, we remark that, since the dotted "warning" line associated to the 30 GeV normalization is actually close to our fiducial model, we can not exclude the possibility that the trend in the high-energy domain ( Figure 5), which is more pronounced with respect to the overall one, could partly be due to a contribution from the UPS population in the inner Galaxy.

Summary
In this paper we have decomposed the gamma-ray emission measured by Fermi-LAT along the Galactic plane (|b| < 20 • ) in the 0.3-230 GeV energy range into different components, including hadronic emission associated to 9 Galactocentric rings, Inverse Compton emission, isotropic diffuse background, point sources and extended sources, and Fermi bubbles. The study was based on a penalized maximum likelihood analysis performed with the SkyFACT package, which is a hybrid of template fitting and image reconstruction methods and allows the templates to vary in both morphology and spectrum.
We provided evidence for a progressive hardening of the gamma-ray hadronic emission towards the inner Galaxy, which corresponds to a hardening in the CR proton spectrum. We showed that this feature is robust with respect to different variations of the analysis setup, and demonstrated that it can be recovered at high energies even if we allow for a break in the power-law fit at 30 GeV. The latter point seems to disfavor an interpretation based on a preeminent role of advection with respect to diffusion in the inner Galaxy due to a significant contribution of streaming instability that lowers the diffusive escape time.
In the second part we estimated the contribution of unresolved point sources via a dedicated Monte Carlo simulation based on our current knowledge of the luminosity function and spatial distribution of Galactic sources (mainly supernova remnants, pulsar wind nebulae and pulsars) and tuned on the 3FGL catalog. We considered a fiducial model and extreme scenarios designed to bracket the large uncertainties involved (in particular on the slope and lower limit of the luminosity function). The ratio between the flux associated to the simulated population of unresolved sources and the diffuse hadronic emission along the Galactic plane stays below 20% for all those models up to 100 GeV. We can safely conclude that the progressive hardening under scrutiny is unlikely to be the result of increasing contamination of unresolved point sources in the inner Galaxy for R > 4.5 kpc, while further study is needed in the very inner regions of the Galaxy.