© T. Siegert et al. 2022

1 Introduction

The diffuse emission spectrum of the Milky Way at photon energies of a few megaelectronvolts (MeV) is one of the least explored phenomena in astrophysics. Despite its rich scientific connections to fundamental nuclear, particle, and cosmic-ray (CR) physics, only one instrument measured the Galactic emission in the 1–30 MeV band: the Compton Telescope (COMP-TEL) onboard the Compton Gamma Ray Observatory (Strong et al. 1999) − more than 20 yr ago. The MeV spectrum provides invaluable and otherwise unavailable insight: the magnitude and shape of the interstellar radiation field (ISRF) is determined through inverse Compton (IC) scattering of giga-electronvolt (GeV) electrons (e.g. Moskalenko & Strong 2000), resulting in an MeV continuum. The low-energy CR spectrum (≲ 100 MeV) outside the Solar System can be measured throughout the Galaxy via nuclear excitation of interstellar medium (ISM) elements, which produce de-excitation γ-ray lines (e.g. Benhabiles-Mezhoud et al. 2013). This is otherwise only possible with the Voyager probes (Stone et al. 2013), which, however, are only sensitive to CR spectra in the local ISM. Annihilation of positrons in flight, which determines the injection energy of their sources in a steady state, shows γ rays from 0.26 MeV up to the particles’ kinetic energy (Beacom & Yüksel 2006). Dark matter candidates could also leave an imprint of their nature in the MeV band (e.g. Boehm et al. 2004; Fortin et al. 2009; Siegert et al. 2022a).

Currently, only one instrument is able to measure the extended emission along the Galactic plane: the Spectrometer aboard the International Gamma-Ray Astrophysics Laboratory, INTEGRAL/SPI (Vedrenne et al. 2003; Winkler et al. 2003). SPI measures photons in the range between 20 keV and 8 MeV through a coded aperture mask. Although INTEGRAL is currently in its 20th mission year and has performed deep exposures in the Galactic bulge and disc, the upper decade of SPI’s spectral bandpass has barely been touched in data analysis.

In this paper we determine the spectrum of diffuse emission in the Milky Way between 0.5 and 8 MeV based on 16 yr of INTEGRAL/SPI data. Our approach is based on spatial template fitting of GALPROP (Strong et al. 2011) models and relies on the success of recent developments in modelling the instrumental background of SPI (Diehl et al. 2018; Siegert et al. 2019). This paper is structured as follows: in Sect. 2 we explain the challenges of SPI data above 2 MeV, the impact of the Sun and Earth’s albedo, and how we handle the background. Our dataset and analysis is presented in Sect. 3. We assess the fit quality and estimate systematic uncertainties in Sect. 4. The resulting spectrum and residuals are found in Sect. 5. We discuss our findings in terms of the Galactic electron population that leads to the IC spectrum and summarise in Sect. 6.

2 SPI data above 2 MeV

In the ‘high-energy’ (HE) range of SPI, between 2 and 8 MeV, only a few targets have been characterised spectrally: the Crab (e.g. Jourdain & Roques 2009, 2020) and the Sun (e.g. Gros et al. 2004; Kiener et al. 2006)¹. The latter has a huge impact on the instrumental background behaviour as a function of time, which is provided by the enhanced particle flux during solar flare events.

SPI data between 2–8 MeV are recorded in 16 384 channels, corresponding to a channel resolution of 0.52 keV (Vedrenne et al. 2003). By default, in official processing (ISDC/OSA; Courvoisier et al. 2003) the HE range is binned into 1keV bins. Per detector, the count rate drops from 10⁻³ to 10⁻⁵ cnts s⁻¹ keV⁻¹ from 2 to 8 MeV. This describes a notoriously noisy spectrum during one observation pointing, lasting typically 0.5–1.0 h, and is one of the main reasons why these data are difficult to analyse.

Most of the measured counts are due to instrumental background radiation, originating from CR interactions with the satellite material. In addition, the γ-ray albedo spectrum from Earth, also induced by CRs, begins to contribute significantly at these energies because SPI’s anti-coincidence shield becomes more and more transparent. Share & Murphy (2001) identified about 20 atmospheric γ-ray lines between 0.5 and 7MeV and their impact on the spectra from the Solar Maximum Mission (SMM). Owing to the composition of Earth’s atmosphere, all these lines are related to either O or N, weak lines of C and B, and the positron annihilation line. The authors also find a significant contribution of unresolved lines as well as an electron bremsstrahlung continuum. The absolute numbers from SMM cannot directly be translated to instrumental background rates for SPI because of INTEGRAL’s eccentric and high inclination orbit as well as its different design. However, the relative rates among the lines and, in particular, between solar quiescence and flares serve as a good proxy for data selection (Sect. 2.1).

Fig. 1 Background count rate of instrumental lines. The atmospheric 12C line at 4.4 MeV shows strong variations (one to three orders of magnitude) when a solar flare occurs. The instrumental 69Ge line at 882.5 keV is barely affected by solar events. We exclude all revolutions in which the 12C rate is 3σ above the running median.

2.1 Solar impact

In Fig. 1 we show the measured count rate of the atmospheric ¹²C line at 4438 keV as a function of the INTEGRAL mission time in units of satellite revolutions around Earth (~3 d). For comparison, we also show an instrumental line of ⁶⁹Ge at 882.5 keV. The long-term behaviour is determined by the solar cycle, being inversely proportional to the sunspot number and therefore the magnetic activity of the Sun (see Diehl et al. 2018, for more examples). The short-term behaviour of the two lines is clearly different: Both atmospheric lines and instrument lines are additionally excited by solar particle events, but, since the shield is more transparent at 4.4 MeV than at 0.9 MeV, the impact of solar events is much stronger for the former. We note that for X-class solar events, such as during INTEGRAL revolutions 128 or 1861, even the 882.5 keV line (and, correspondingly, most other lines) showed a significant increase with respect to the running mean. For the ¹²C line, and similar lines such as ¹⁶O (6129 keV) or ¹⁴N (3674 keV), solar flares increase the received count rate by up to three orders of magnitude.

We wanted to avoid entire revolutions in our data selection and background modelling altogether and used the ¹²C line, which shows the largest rate ratios between flares and quiescence, as a proxy for enhanced short-term solar activity. We applied a running median of 30 revolutions to estimate the baseline rate at 4438 keV. Then, we removed any revolution from our data in which the measured ¹²C line rate is more than three standard deviations above the median rate. This was applied to all energies and processing chains (see Sect. 3.1).

2.2 Handling instrumental background

Based on the reduced HE database, with solar particle events filtered out as described in Sect. 2.1, we applied the method from Siegert et al. (2019) to construct a high spectral resolution instrumental background database. First, we integrated over the entire filtered data archive and all detectors to also identify the weakest lines. We found 610 lines with rates per detector between 10⁻⁷ and 10⁻² cnts s⁻¹. We then split the energy range between 2 and 8 MeV into multiple smaller bands to determine the spectral parameters of the background lines on top of a multiple broken power law for each detector. The integration time to extract the spectral information was set to the time interval between two annealing periods, which is typically half a year. This has the advantage of enough counts per spectrum to determine the spectral shapes reliably, but it only estimates an average degradation of the detectors, typically broadening lines by up to 15% between two annealings.

In Fig. 2 we show the spectrum of detector 00 measured between INTEGRAL revolutions 777 and 795, together with the fit to determine the flux ratios for the final background model, and the residuals. Over the full energy range, this method provides adequate fits. In the right panel of Fig. 2 we show the zoomed-in version between 2.4 and 2.9 MeV and detail the instrumental background lines.

In Siegert et al. (2019), this technique was applied to the diffuse emission of the 511 keV and 1809 keV lines and to pointlike emission from continuum sources beyond 3 MeV. Here we extend this approach to diffuse continuum emission up to the boundaries of the SPI HE response of 8 MeV. While the energy range and source function in this work is a new application to the previous method, the expected signal-to-background count ratio (S/B) per energy bin is about the same as in the case for the 511 keV line, for example. The line shows an integrated flux of 10⁻³phcm⁻² s⁻¹ above an average background count rate of 0.6 cnts s⁻¹. Taking the background spectrum shown in Fig. 2 as representative for the whole mission, (S/B)511 is about 1.6 × 10⁻³. This value changes by about 50% over the course of the dataset (cf. the background variation in Fig. 1). While the expected signal in the continuum from 0.5 to 8 MeV decreases significantly from ~10⁻⁵ to ~10⁻⁷phcm^{−2 s−1} keV⁻¹ with a power-law index around −1.7, the background count rate also drops sharply with an index of −3 between 0.5 and 5 MeV. Depending on energy, the expected IC (S/B)_IC varies between 1 and 8 × 10⁻³. Judging from this ratio alone, using the Siegert et al. (2019) method seems justified.

In order to determine the temporal variability in the background per energy bin, we followed the same approach as in Siegert et al. (2019). Based on the fact that the germanium detector rate alone is insufficient to model or predict the pointing-to-pointing variation, an onboard radiation monitor is typically used to fix the background behaviour. We used the rate of saturating germanium detector events (GeDSat) to link the background amplitudes in time and employed it as a ‘tracer’ function. It was shown in several studies that used either this or previous methods that neither the GeDSat rate alone nor orthogonalised additional tracers can fully explain the background variability. To account for unexplained variance, we split the background tracer in time and set regularly spaced nodes to re-scale the background model. Because this choice is not unique, we attempted a tradeoff between the number of additional background parameters required and the likelihood. This is achieved by the use of the Akaike information criterion (AIC; Akaike 1974; Burnham & Anderson 2004), $$\text{AIC=2}\left(n_{\text{par}}-\text{ln}\left(\widehat{\mathcal{L}}\right)\right),$$(1)

with n_par being the number of fitted parameters and $\widehat{\mathcal{L}}$ the log-likelihood maximum, to determine which configuration of background variability is optimal for each energy bin. We tested background variability timescales between 0.19 d (1/16 of an orbit) and 30 d (ten orbits) for each energy bin and identified the optimum AIC (see Table 1).

The AIC has no absolute meaning, but its relative values can be used to identify similarly likely model configurations. We used the AIC optimum value and defined a threshold based on the required number of fitted parameters to select background variability timescales that also provide an adequate fit. Likewise, we can estimate a systematic uncertainty on the extracted flux per energy bin using the AIC (see Sect. 4.2).

Fig. 2 SPI data of one detector (00) between INTEGRAL revolutions 777 and 795 and spectral fits. Left: complete spectrum of SPI’s HE range between 2 and 8 MeV (top) and residuals (bottom). Right: zoomed-in view between 2.4 and 2.9 MeV, with individual lines indicated. The shown normalised residuals scatter around 0.04 with a standard deviation of 0.94 for the 6000 data points indicate a sufficiently well-described background spectrum fit.

3 Data and analysis

3.1 Filtered dataset

Based on the considerations in Sect. 2, we defined 12 logarithmic energy bins and considered INTEGRAL revolutions between 43 (February 2003) and 2047 (January 2019). Details about the included revolutions and the number of observations are provided in Appendix B. The dataset covers two independent SPI processing chains: pulse-shaped-discriminated (PSD) events between 0.5 and 2.0 MeV and HE between 2 and 8 MeV. To combine the extracted data points into one common spectrum, the PSD fluxes were scaled by the expected loss in efficiency of 1/0.85 due to increased dead time. The PSD range includes the Galactic diffuse ²⁶Al line at 1808.74 keV, which we included in a narrow bin between 1805 and 1813 keV.

We focused on a spatial region around the Galactic centre that is covered by targeted observations (pointings) falling into −40° ≤ℓ≤ 40°, −40° ≤ b ≤ 40°. Because of SPI’s fully coded 16° × 16° field of view, we considered diffuse emission out to |ℓ| ≤ 47.5° and |b| ≤ 47.5°, respectively. This avoids the partially coded field of view and its edge effects when the exposure is either very small (few pointings) or shows large gradients. Our flux estimates are therefore normalised to a spherical square with a side length of 95°, covering a solid angle of Ω = 2.43 sr.

Other data selections included radiation monitors and orbit parameters: we only chose pointings in the orbit phase between 0.15 and 0.85 to avoid residual activation by the Van Allen radiation belts. Whenever the running mean of the rate ratio between the anti-coincidence shield and the total rate of the Ge detectors exceeded a 3σ threshold, we excluded the observation. Pointings with a cooling plate difference of more than 0.8 K were also excluded. Finally, revolutions 1554–1558 were removed due to the outburst of the microquasar V404 Cygni.

These selections resulted in a total of 36103 pointings for the PSD and HE ranges. Based on a background-only fit to the selected data, we investigated the residuals as a function of pointing, detector, and energy, and removed individual observations whose deviations were larger than 7σ. Given the expectedly low signals, any diffuse emission contribution is about 0.1−1.0% of the total counts and would never distort the residuals in the broad logarithmic energy bins. The additional filter removed 0.6% of the PSD data, for a reduced dataset of 35 892 pointings. The HE range shows no such outliers. In total, the dead-time-corrected exposure time of our dataset is 68.5 Ms for a working detector.

The characteristics of our dataset are found in Table1, including the number of data points, the background variability timescale per energy bin, the number of degrees of freedom (d.o.f.) per energy bin, and a calculated goodness of fit criterion. As described in Sect. 2.2, the background variability is determined to first order by the GeDSat rate. Because this tracer is not sufficient to describe the true (measured) background variability, we inserted regularly spaced time nodes to capture the unexplained variance. As shown in Siegert et al. (2019), the number of time nodes, or in turn the background variability, depends on the energy, the bin width, and to some extent the source strength. With an optimisation to require the fewest number of parameters while at the same time obtaining the best likelihood (see the AIC approach in Sect. 2.2), this timescale was determined for each energy bin individually, always taking a baseline sky model into account (see Sect. 3.3). The background variability not explained by the tracer alone therefore changes between 0.75 and 6 days, increasing roughly with energy.

3.2 General method

SPI data analysis relies on a comparison between the raw count data per pointing, detector, and energy, with a combination of instrumental background and celestial emission. We modelled the data d_p per pointing p for each energy bin individually as $$m_p={\displaystyle \sum _t{\displaystyle \sum _j{R}_{jp}{\displaystyle \sum _{k=1}^{N_S}{\theta }_{k,t}{M}_{kj}+}{\displaystyle \sum _{t^{\prime }}{\displaystyle \sum _{\begin{array}{l}k=N_S+1\\ =\end{array}}^{N_S+N_B}{\theta }_{k,t^{\prime }}{B}_{kp}}}}},$$(2)

where the response R_jp is applied to each of the k = 1… N_S sky models M_kj pixelised by j. The N_B background models B_kp are independent of the response. The only free parameters of this model are the amplitudes θ_k,t and θ_k,t′ of the sky and background models, respectively. They were estimated through a maximum likelihood fit subject to the Poisson statistics $$\widehat{\mathcal{L}}(\theta |D)={\displaystyle \prod _{p=1}^{N_{\text{obs}}}\frac{m_p^{d_p}\exp (-m_p)}{d_p!}},$$(3)

where $D=\left\{d_1,\dots ,d_{N_{\text{obs}}}\right\}$ is the dataset of measured counts per pointing.

Both sky and background were allowed to change on different timescales, t and t′, respectively. Source variability above 500 keV is too faint to be detected in this dataset, and we assumed all sources as well as the diffuse emission to be constant in time. We followed the approach of Siegert et al. (2019) to model the instrumental background from the constructed line and continuum database (Sect. 2.2). We built two background models per analysis bin from the newly constructed HE background database, one for the instrumental lines and one for the instrumental continuum. The amplitudes of these models were fitted together with the flux(es) of expected emission model(s) (see Sect. 3.3). Any background variation that is not covered by this tracer was refined by additional time nodes to re-scale the GeDSat tracer function. The estimated background variability timescale, changing from 0.75 d (-1/4 of an orbit) between 0.5 and 1.4 MeV, up to 3–6 d above 4 MeV, is equivalent to ~2500 and ~500 fitted parameters, respectively.

The maximum likelihood fits to the raw data were performed with OSA/spimodfit (Halloin 2009; Strong et al. 2005). The extracted flux data points were governed by an energy redistribution matrix to take the instrument dispersion into account. Spectral fits were performed with 3ML (Vianello et al. 2015).

3.3 Emission templates

Four resolved point sources, 1E 1740.7–2942, GRS 1758–258, IGRJ17475–2822, and SWIFTJ1753.5−0127, are expected in addition to the diffuse emission (Bouchet et al. 2011). We modelled the point sources as constant in time up to an energy of 850 keV. The 1.8 MeV line from ²⁶Al was included as the SPI maximum entropy map by Bouchet et al. (2015). For the continuum, only the leptonic emission is relevant at the energies of interest, well below the so-called π⁰ bump. Previous analyses that included INTEGRAL, COMPTEL, and data from Fermi’s Large Area Telescope (LAT) suggested that electron bremsstrahlung is sub-leading by at least an order of magnitude in our range of energies (Strong 2011), and we neglect it in the following. We used energy-dependent IC scattering emission templates from the GALPROP (v56) CR propagation code (Strong et al. 2011). In particular, we used: (i) the model ^SL^Z4^R20^T150^C5 adopted in Ackermann et al. (2012), which reproduces Fermi/LAT gamma-ray observations well; and (ii) the model in Table 2 of Bisschoff et al. (2019), which adjusts primary CR spectra and propagation parameters, also accounting for Voyager 1 data. This latter model adopts different electron spectral indices, as well as diffusion scaling with rigidity below and above a reference rigidity of 4 GV, as encoded in the spectral indices δ₁ and δ₂, whose default values are 0.3 and 0.4, respectively.

By using the IC template as a tracer of the diffuse emission, we avoided adopting generic descriptions of the emission with, for instance, exponential discs or 2D Gaussians and/or continuum tracer maps such as the COMPTEL 1–30 MeV map (Strong et al. 1999). Such a model-inspired approach has a twofold advantage: first, our choice to resort to absolute models allowed us to gauge how far ‘off’ the fluxes are from typical expectations. Second, the predicted IC morphologies spanned by these models, together with the flexibility in the background following from the number of time nodes required, provide a measure of systematic uncertainties.

To this end, we tested different variants of the Voyager CR parameter configuration (Bisschoff et al. 2019) and assessed the magnitude of systematic uncertainties. We defined: (a) δ₁ = 0 to represent the possibly flatter behaviour of the diffusion scaling at low rigidities (Genolini et al. 2019); (b) δ₁ = δ₂ = 0.5 to test the effect of a single diffusion index closer to current best fits of the ratio of secondary to primary CR nuclei (Genolini et al. 2019; Weinrich et al. 2020); (c) 10×opt to account for a factor of 10 stronger optical ISRF, corresponding to a possible enhancement of this poorly known component of the ISM towards the inner Galaxy (Bouchet et al. 2011); and (d) thick halo − we adopted the extreme value L = 8 kpc halo half thickness, as opposed to the default 4 kpc, and we re-normalised the diffusion coefficient accordingly to account for their well-known degeneracy (Weinrich et al. 2020).

These variants affect both the spatial distribution of the IC photons (i.e. the morphology) and the spectral IC shape. It is suggested that the primary CR electron spectrum has a break around E_e = 2.2 GeV, changing from a power-law index of p₁ = −1.6 to p₂ = −2.4 (Fermi/LAT), or E_e = 4.0 GeV with power-law indices of p₁ = -1.9 and p₂ = −2.7 (Voyager). This implies a spectral break in the photon spectrum from ${\alpha }_1^{\text{IC}}=-1.3$ to ${\alpha }_2^{\text{IC}}=-1.7$ around 25 keV, 250 keV, and 25 MeV for the individual components of the ISRF (cosmic microwave background ~0.001 eV, dust ~0.01 eV, star light ~1 eV). The components were weighted with their respective spatial intensities, which results in a power-law-like spectrum with an index of α^IC = −1.4 to −1.3 up to a few MeV. We also notice that the photon spectrum curvature maximum can be shifted from around 20 MeV to around 3–5 MeV by increasing the optical ISRF by a factor of 10 (10×opt), which also increases the flux by at least a factor of 5. However, such a model might fall short in describing the absolute flux, especially below 4 MeV, and shows a steeper spectrum in our analysis band than what has previously been measured.

4 Fit quality and systematic uncertainties

4.1 Fit quality

We judged the adequacy of our maximum likelihood fits in each energy bin by the shape and distribution of the normalised residuals, $r=\left(d-m\right)/\sqrt{m},$, with data d and model m, as a function of time (pointing). To a lesser extent, mainly because the value has no proper meaning in this context but is frequently used in the literature, we considered the reduced χ² value of our fits, ${\chi }^2/\text{d}\text{.o}\text{.f}\text{.}={\displaystyle {\sum }_i{r}_i^2/\text{d}\text{.o}\text{.f}\text{.}}$ We refer the reader to Andrae et al. (2010) for why the use of χ² can be misleading in general, and in particular in the context of this work.

A ‘bad fit’ would be immediately seen in the residuals − even though χ²/d.o.f. might be close to the optimal value of 1.0. For example, individual outliers of even 50σ would still result in a reduced χ² value close to 1.0 but would distort the entire fit results. Likewise, an apparently large or small reduced χ² value must not be considered ‘bad’ in the case of Poisson statistics because it is only a calculated value and not related to the actual data generating process. Therefore, for a ‘good fit’ we demand the temporal sequence of residuals to show no individually strong outliers (≳10σ) and no clustered weak outliers (many neighbouring values above or below the mean).

In Fig. B.1 we show as an example the complete sequence of residuals of the camera combined and all individual detectors for the energy range 6–8 MeV. The reduced χ² value of this fit evaluates to 1.00383 with 581 390 d.o.f. Since we find no remaining structure in these residuals, we deem this fit adequate. This is also true for the remaining energy bins analysed in this work.

4.2 Background systematics

For each energy bin we calculated the standard deviation of flux estimates that follow ΔAIC ≤ n_par(ΔAIC = 0) to estimate our systematic uncertainties. This inequality still demands that the fit must be ‘good’ in the terms described above (Sect. 4.1) but allows the fitted parameter of interest (the amplitude of the sky model) to vary within a reasonable range. It does not describe another statistical uncertainty because the number of total parameters is increased when a new, smaller time variability scale is introduced. The likelihood would always increase towards a ‘better’ fit, which is why we used the AIC again to take the changing number of d.o.f. into account.

Since the pointing-to-pointing variation in our background model is fixed by an onboard monitor (Siegert et al. 2019) and consequently not entirely perfect, the background is re-scaled (fitted) according to the selected time nodes (Sect. 2.2). Because the sky components are either localised (point sources) or show gradients (diffuse emission), it is insufficient to scale this background model once for the entire dataset; it requires additional time nodes. This was realised in the maximum likelihood fit via the introduction of the background variability timescales, t′, or, equivalently, more background parameters, θ_k,t′ This re-scaling depends on energy, bin size, flux, and time (pointing) because different source strengths determine the total counts and because the INTEGRAL observation scheme is not a survey but pointed according to granted proposals.

We show two examples in Fig. 3 of how the AIC changes as a function of the background variability timescale and how the systematic uncertainties are estimated from this search. The energy bin from 2000–2440 keV is dominated by statistical uncertainties, meaning many, also unlikely, background model configurations result in the same flux estimate, given the same spatial template. From 4418 to 5945 keV, the systematic uncertainty is of the order of the statistical uncertainty because the background variability allows a larger range of flux estimates.

4.3 Source systematics

We used the different emission templates described in Sect. 3.3 to estimate another source of systematic uncertainties from the IC emission itself. Because the emission is expected to be weak, we refrained from extensive parameter scans and instead used a set of parameters that we explored within their uncertainties. Our results and extracted fluxes can thus be used in follow-up studies to constrain the CR propagation parameters. In total, we tested six different setups, one best-fit model from Fermi/LAT analyses, ^SL^Z4^R20^T150^C5 (Ackermann et al. 2012), and five variants of the combined study from Voyager, Fermi/LAT, and the Alpha Magnetic Spectrometer experiment (AMS-02) data from Bisschoff et al. (2019). We list the systematics according to different spatial models in Table 2.

5 Results

5.1 Spatial residuals

For a visual verification that we indeed measured emission from the Galactic plane, we fitted a background-only model to the data. The residuals of these fits are projected back onto the celestial sphere such that we obtain an image of where the residual counts are found. We caution that this is not an image reconstruction, nor should individual features be over-interpreted: the backward application of the coded-mask response to the residual counts is not unique and is limited by the source strength. If positive (or negative) regions consistently cluster in these residual images in the same areas as a function of energy, we can conclude that the measured fluxes are less likely to be an instrumental artefact. If instrumental background lines are not modelled properly, they will appear as residuals in these images but be restricted to one particular energy.

Figure 4 shows the residuals of a background-only fit and the changed appearance after including the IC template maps. We find positive residuals clustered in the region of the Galactic centre and disk for all energies. The magnitude of the residuals decreases with energy, as expected from the power-law behaviour of the IC emission. The residuals that include the IC template maps are devoid of the central enhancement and show a wreath-like pattern. This originates from large gradients in the exposure map, dropping from long observed regions to nearly zero within a few degrees.

We conducted an additional test for the detection of diffuse emission from 0.5 to 8.0 MeV by altering the IC sky model. If the emission is due to an instrumental effect and not from the Galactic plane, a similar spectrum (that is, similar to that of the background) will result if the template map used has no impact on the fit. We tested such a scenario by shifting the IC template maps for each energy bin by +20° in latitude and repeated the fit. The resulting spectra for both cases are shown in Fig. 5. Clearly, the spectrum follows a power-law shape for the template centred on the Galactic plane and is consistent with zero flux for the shifted template. We conclude that there is indeed diffuse emission detected by SPI in the Galactic plane up to 8 MeV.

Fig. 3 Background variability and systematics for two chosen energy bins, 2000–2440 keV (top, dominated by statistics) and 4418–5945 keV (bottom, statistics and systematics of the same magnitude). Left: scan for optimal background variability timescale as measured with the AIC (right axis), Eq. (1). The optimum is found at timescales of 6 and 3 d, respectively corresponding to two and one INTEGRAL orbits. For comparison, the calculated reduced Pearson χ2 is shown for each tested grid point. Right: corresponding flux estimates and statistical uncertainties (orange). The systematics are estimated from the standard deviation of fluxes whose ΔAIC values are below the number of fitted parameters at optimum AIC (shaded region).

Fig. 4 Count residuals projected back onto the sky as a function of energy. Top: background-only fits with consistent positive residuals along the Galactic plane and centre. Bottom: background plus IC template map fits, with the exposure map indicated. Shown are levels where the exposure drops to 50 and 10% of the maximum, respectively.

Fig. 5 Comparison of fitted IC source fluxes with GALPROP model δ1 =δ2 = 0.5 (left) and the same template maps but shifted in latitude by +25° (right). The expected spectrum is shown as a black curve and the extracted fluxes with crosses. The insets show a representative template map, with the unshifted map as contours in the right inset. No excess is found for the shifted templates, consolidating the signal found from the Galactic plane.

5.2 Spectrum

In Fig. 6 we show the extracted data points from our analysis of the IC emission. As expected, the²⁶ Al line at 1.8 MeV has no spatial component following the IC morphology, and we provide an upper limit. For a visual comparison to the 20 yr old COMPTEL data points (Strong et al. 1999), we binned our flux data points to a minimum signal-to-noise ratio of 6. We note that a comparison of ‘extracted fluxes’ from different instruments without taking the spectral response into account can be (and most of the time is) misleading. Nevertheless, it can provide a general overview of the consistency between the measurements. We find an excellent agreement between SPI and COMPTEL in the overlap region from 1–8 MeV and show that after 16 yr in space, SPI’s diffuse continuum measurements have smaller uncertainties than those of COMPTEL.

We fitted the spectrum phenomenologically with a power law, C₀(E/1 MeV)^α. Our best-fit parameters are a flux density of C₀ = (3.1 ± 0.3) × 10⁻⁶ ph cm⁻² s⁻¹ keV⁻¹ sr⁻¹ at 1 MeV and an index of α = -1.39 ± 0.09. The spectral index is consistent with the work by Bouchet et al. (2011), who found an index of 1.4–1.5 between 0.02 and 2.4 MeV. Extrapolating the fitted power law to the COMPTEL band up to 30 MeV and propagating the spectral uncertainties also shows a general agreement (violet band). Using instead a cutoff power law with a normal prior for the break energy of 4.0± 1.8 MeV (cf. Table 4 in Bouchet et al. 2011) leads to slightly larger flux values between 1 and 4 MeV and to slightly smaller fluxes (≲10% difference in both cases) elsewhere (red band). The resulting power-law index is then −0.95 ± 0.16 and the fitted break energy 4.9 ± 1.4 MeV.

We note that SPI also detects photons above 8 MeV; however, the official tools do not provide an imaging or spectral response at these energies. On the other hand, the SPI spectrum below 0.5 MeV is already well determined, and we refer the reader to Bouchet et al. (2011) and Siegert et al. (2022b) for details about this low-energy band. Extending the spectrum in either direction is beyond the scope of this paper.

As an alternative to a generic power law, we compare the extracted data points from each GALPROP IC morphology to the expected absolute model in Fig. 7. In this way, we can determine which propagation model provides the best absolute normalisation when compared to SPI data. The magnitudes of the systematic uncertainties were calculated as the mean absolute difference from the extracted flux values among the tested IC morphologies (thin error bars). The flux values (crosses) in Fig. 7 and their statistical uncertainties (thick error bars) are the means of the individually extracted fluxes (see also Table 2). At the spectral level, excessively extreme variations in the diffusive properties, as in the model δ₁ = 0, appear in tension with the data. All other models seem to lead to quasi-parallel spectra, in broad agreement with the deduced shape.

We note that the default predictions are about a factor of 2–3 below the measured fluxes, with increasing discrepancy towards higher energies. However, the plot also shows that it is hard to pin down the origin of the mismatch: Variations in the diffusion properties, variations in the photon targets by a factor of 3–5, or variations in the CR source spectra by a similar factor (not shown) could be involved in explaining the mismatch. Orlando (2018) argue, however, that this last option, also invoked in Bouchet et al. (2011), would lead to an overproduction of synchrotron emission, which disfavours such a hypothesis. Also, there is almost no sensitivity to the halo thickness (thick halo). The best match is found for the model variant with δ₁ = δ₂ = 0.5 that assumes a constant diffusion coefficient index for the entire CR electron spectrum. Finally, we note that part of the emission could be due to an unresolved population of Galactic sources, indistinguishable from a continuum emission. Such sources might show spectra similar to the ‘hard tails’ that have recently been detected in a few X-ray binaries (e.g. Cangemi et al. 2021a,b). Emission up to ~500 keV and beyond has been observed in individual sources, which could flatten out the cumulative spectrum of a population of weak sources. In term of the energy flux, E²F_E, this could lead to a peak in the unresolved point source spectrum around 0.5–3.0 MeV, depending on the objects’ properties and their luminosity function in the Milky Way.

Evaluating the different model variants, we find that the systematic uncertainties due to the background variability range between 5% (0.5–0.85 MeV) and 20% (3.3–8.0MeV). The systematic uncertainty from the IC morphology ranges between 20 and 30%.

Fig. 6 Spectrum of the analysed region between 0.5 and 8 MeV. The orange data points show the extracted fluxes from the energy-dependent IC template δ1 = δ2 = 0.5, and the fuchsia points a re-binning to a minimum signal-to-noise ratio of 6. The fitted power-law spectrum (F0.5–8.0 = (5.7 ± 0.8) × 10−8 ergcm2 s−1 α = −1.39 ± 0.09) is shown with its 68.3, 95.4, and 99.7 percentile bands in violet, and the fitted cutoff power law with EC = 4.9 ± 1.4 MeV in red. We compare the fluxes of this work with historic measurements by COMPTEL (green; Strong et al. 1999).

Fig. 7 Extracted spectrum (black crosses) with statistical (thick error bars) and systematic (thin error bars) uncertainties, as well as a generic power-law fit (band).

6 Summary, discussion, and conclusions

For the first time in 20 yr, we have provided a description of the Galactic diffuse γ-ray spectrum up to 8 MeV. Our results are compatible with previous estimates from COMPTEL and finally supersede its precision as measured by the signal-to-noise ratio. The spectrum is adequately described empirically by a power law with an index of −1.39 ± 0.09 and a flux of (5.7 ± 0.8_stat ± 1.7_syst) × 10⁻⁸ erg cm⁻² s⁻¹ between 0.5 and 8.0MeV. Our general finding is inline with Bouchet et al. (2011), showing the need for a continuum emission broadly peaking in the inner Galaxy and compatible in spectrum with the expected IC scattering of CR electrons onto the ISRF Such a model, however, overshoots baseline expectations of state-of-the-art models calibrated to local Voyager 1 and AMS-02 data by a factor of 2–3. With dedicated GALPROP runs, we discussed how enhanced ISRF in the Galactic centre or modified diffusion may be responsible for a similar discrepancy. A propagation model with a single diffusion index of δ₁ = δ₂ = 0.5 provides the best description of the SPI data in the photon energy range between 0.5 and 8.0MeV. Our analysis also includes an assessment of systematic uncertainties based on realistic morphologies of IC models, which lead to a systematic flux uncertainty from the IC spatial distribution of between 20 and 30%.

An O(10%) sub-leading bremsstrahlung component with a less steep electron spectrum (Strong et al. 2000,Strong et al. 2005; Bouchet et al. 2011; Ackermann et al. 2012) can further improve the agreement between CR propagation model expectations and data. Our improved estimates of the MeV spectrum in the Milky Way for broadband γ-ray analysis will provide more stringent estimates of the Galactic electron population at GeV energies.

Nonetheless, a better sensitivity in the MeV range, and therefore a future mission covering the MeV sensitivity gap, such as the recently selected small explorer mission COSI, the Compton Spectrometer and Imager² (Tomsick et al. 2019), will shed further light on the possibilities of additional continuum sources in lieu of true diffuse emission. This will be of relevance not only for the astrophysical study of Galactic CR populations, but also for searches of more exotic, beyond-the-standard-model emission processes, such as from dark matter candidates (Alves Batista et al. 2021).

The spectral data points and response are available in an online repository³. We encourage the use of this renewed dataset from INTEGRAL/SPI for comparisons to Galactic emission processes.

Acknowledgements

T.S. is supported by the German Research Foundation (DFG-Forschungsstipendium SI 2502/3-1) and acknowledges support by the Bundesministerium für Wirtschaft und Energie via the Deutsches Zentrum für Luft- und Raumfahrt (DLR) under contract number 50 OX 2201. F.C., J.B. and P.D.S. acknowledge support by the “Agence Nationale de la Recherche”, grant n ANR-19-CE31-0005-01 (PI: F. Calore).

Appendix A Spectral fits

The flux uncertainties are approximately Gaussian, so the log-likelihood for the following fits is proportional to χ². In Fig. A.1 we show as an example the spectral fit of a power law, C₀(E/1 MeV)^α, to the extracted data points for the morphology with diffusion indices δ₁ = δ₂ = 0.5, taking the energy redistribution matrix into account. At optimum, we find a χ² of 14.6 with ten dof. The spectral parameters are C₀ = (7.6 ± 0.6) × l0⁻⁶phcm⁻²s⁻¹ keV⁻¹ at 1 MeV and an index of α = −1.39 ± 0.09. The extracted data points for other morphologies are similar, differing by at most 1σ for individual energy bins. We show the fitted parameters for the alternative models in Table 2, from which we estimate systematic uncertainties. We chose the mean of all fits as a baseline to estimate systematic uncertainties and quoted the variant δ₁ = δ₂ = 0.5 for statistical uncertainties because it shows the individually highest likelihood across all analysed energy bins. We then picked the maximum differences to the mean values as systematics from the emission morphology for all parameters in Table 2. Using the above considerations from the background variability timescale, we also show the systematics from this component in the same table.

Fig. A.1 Spectral fit of Galactic diffuse emission between 0.5 and 8 MeV for morphology variant δ1 =δ2 = 0.5. Top: Extracted data points (orange) and fitted power-law model (purple; 68th and 95th percentiles) and cutoff power law (red). Bottom: Normalised residuals.

Appendix B Additional tables and figures

The selected observations from the filtered SPI database that excludes strong solar activity are listed in Table B.1. In Table 1 we summarise the energy bins used in this analysis and details about the instrumental background.

Fig. B.1 Fit residuals as a function of time (pointing). Shown are the 36103 data points, summed over the detector array (top left), and the 19 individual detectors for the energy bin 5945–8000 keV. The residuals are histogrammed in the panels on the right, together with the expected normal distribution of zero mean and unit standard deviation (red), indicating an adequate fit.

References

All Tables

All Figures

	Fig. 1 Background count rate of instrumental lines. The atmospheric 12C line at 4.4 MeV shows strong variations (one to three orders of magnitude) when a solar flare occurs. The instrumental 69Ge line at 882.5 keV is barely affected by solar events. We exclude all revolutions in which the 12C rate is 3σ above the running median.
In the text

Fig. 2 SPI data of one detector (00) between INTEGRAL revolutions 777 and 795 and spectral fits. Left: complete spectrum of SPI’s HE range between 2 and 8 MeV (top) and residuals (bottom). Right: zoomed-in view between 2.4 and 2.9 MeV, with individual lines indicated. The shown normalised residuals scatter around 0.04 with a standard deviation of 0.94 for the 6000 data points indicate a sufficiently well-described background spectrum fit.

In the text

Fig. 3 Background variability and systematics for two chosen energy bins, 2000–2440 keV (top, dominated by statistics) and 4418–5945 keV (bottom, statistics and systematics of the same magnitude). Left: scan for optimal background variability timescale as measured with the AIC (right axis), Eq. (1). The optimum is found at timescales of 6 and 3 d, respectively corresponding to two and one INTEGRAL orbits. For comparison, the calculated reduced Pearson χ2 is shown for each tested grid point. Right: corresponding flux estimates and statistical uncertainties (orange). The systematics are estimated from the standard deviation of fluxes whose ΔAIC values are below the number of fitted parameters at optimum AIC (shaded region).

In the text

	Fig. 4 Count residuals projected back onto the sky as a function of energy. Top: background-only fits with consistent positive residuals along the Galactic plane and centre. Bottom: background plus IC template map fits, with the exposure map indicated. Shown are levels where the exposure drops to 50 and 10% of the maximum, respectively.
In the text

Fig. 5 Comparison of fitted IC source fluxes with GALPROP model δ1 =δ2 = 0.5 (left) and the same template maps but shifted in latitude by +25° (right). The expected spectrum is shown as a black curve and the extracted fluxes with crosses. The insets show a representative template map, with the unshifted map as contours in the right inset. No excess is found for the shifted templates, consolidating the signal found from the Galactic plane.

In the text

Fig. 6 Spectrum of the analysed region between 0.5 and 8 MeV. The orange data points show the extracted fluxes from the energy-dependent IC template δ1 = δ2 = 0.5, and the fuchsia points a re-binning to a minimum signal-to-noise ratio of 6. The fitted power-law spectrum (F0.5–8.0 = (5.7 ± 0.8) × 10−8 ergcm2 s−1 α = −1.39 ± 0.09) is shown with its 68.3, 95.4, and 99.7 percentile bands in violet, and the fitted cutoff power law with EC = 4.9 ± 1.4 MeV in red. We compare the fluxes of this work with historic measurements by COMPTEL (green; Strong et al. 1999).

In the text

	Fig. 7 Extracted spectrum (black crosses) with statistical (thick error bars) and systematic (thin error bars) uncertainties, as well as a generic power-law fit (band).
In the text

	Fig. A.1 Spectral fit of Galactic diffuse emission between 0.5 and 8 MeV for morphology variant δ1 =δ2 = 0.5. Top: Extracted data points (orange) and fitted power-law model (purple; 68th and 95th percentiles) and cutoff power law (red). Bottom: Normalised residuals.
In the text

	Fig. B.1 Fit residuals as a function of time (pointing). Shown are the 36103 data points, summed over the detector array (top left), and the 19 individual detectors for the energy bin 5945–8000 keV. The residuals are histogrammed in the panels on the right, together with the expected normal distribution of zero mean and unit standard deviation (red), indicating an adequate fit.
In the text