© M. M. M. Pleintinger et al. 2019

1. Introduction

²⁶Al is an ideal tracer of ongoing nucleosynthesis in the Galaxy. It is produced in massive stars and ejected into their surroundings via stellar winds and supernovae (SNe). It decays with a half-life time of ∼0.7 Myr to ²⁶Mg and emits a photon at 1809 keV, which can be measured by gamma-ray telescopes. The spatial distribution of ²⁶Al provides information about active sites of nucleosynthesis and galactic chemical enrichment, as well as dynamics and feedback processes in the interstellar medium (ISM) throughout the Milky Way (e.g. Diehl et al. 2006, 2010; Wang et al. 2009; Kretschmer et al. 2013; Bouchet et al. 2015; Siegert & Diehl 2017).

Hydrodynamic simulations are a crucial tool for understanding the dynamics and chemical enrichment of galaxies, and for interpreting observations. Empirical models that may be used to describe diffuse radioactivity in the ISM of the Galaxy have a long history (cf. Prantzos 1993; Prantzos & Diehl 1995; Knödlseder et al. 1996; Lentz et al. 1999; Sturner 2001; Drimmel 2002; Alexis et al. 2014). Yet it is more recently that Fujimoto et al. (2018) reported the first galactic-scale hydrodynamic simulation that starts from basic physical processes with the aim of tracking the synthesis and transport of radioactive isotopes, such as ²⁶Al and ⁶⁰Fe in a Milky Way-like galaxy.

Previous, heuristic model comparisons identify morphological similarities between ²⁶Al emission and multi-wavelength tracers or geometric emission models (e.g. Hartmann 1994; Prantzos & Diehl 1995; Diehl et al. 1997, 2004; Knödlseder et al. 1999; Kretschmer et al. 2013), leaving astrophysical implications to their interpretations. It is important to cross-check simulations that are based on astrophysical assumptions with observations. A fundamentally informative comparison of hydrodynamical simulations to actual measurements is challenging because the Milky Way is one particular realisation of a galaxy and any given hydrodynamic simulation will, even if it is intended to be similar, not necessarily match it perfectly. This is further complicated by the observational limitations of gamma-ray data due to the necessity of image reconstruction methods as compared with direct imaging.

In this paper, we investigate a range of methodological approaches for a generalised comparison of ²⁶Al full-sky emission maps from the simulation performed by Fujimoto et al. (2018) to gamma-ray data measured with the Spectrometer on INTEGRAL (SPI) (Vedrenne et al. 2003) aboard the International Gamma-Ray Astrophysics Laboratory (INTEGRAL) satellite (Winkler et al. 2003). In Sect. 2, we present the data analysis procedure for SPI observations, as well as the properties of the simulation by Fujimoto et al. (2018). Different comparison methods are described in Sect. 3, concluded by a discussion of the results in Sect. 4.

2. Observations and simulations

2.1. Gamma-ray measurements

Imaging with gamma-ray telescopes, either coded-mask based such as SPI or Compton telescopes, such as the Imaging Compton Telescope (COMPTEL) aboard the Compton Gamma Ray Observatory (CGRO) satellite, is impacted by the large instrumental background due to cosmic-ray bombardment. Typically, background amounts to 90–99% of measured events in SPI so that direct imaging is only possible for strong sources. On the other hand, maximum likelihood and maximum entropy approaches are well-tested and can be directly applied to the raw data, including an elaborate background model. The first full-sky image of the 1.8 MeV emission was obtained by Oberlack et al. (1996) using a maximum entropy deconvolution applied to data from 3.5 yr of observations with COMPTEL. In this work, we used the map from the entire nine-year CGRO mission obtained with the same method (Fig. 1, Plüschke et al. 2001). The COMPTEL telescope has an angular resolution of 3.8°. The map shows a large latitude extent on top of a clumpy structure concentrated to the inner galactic disk region and associated with spiral arm tangents as well as nearby massive star regions.

Fig. 1. Maximum entropy all-sky map of 1.8 MeV emission in the Milky Way obtained from 9 yr of observations with COMPTEL (Plüschke et al. 2001).

For our analysis, we used observations of this ²⁶Al gamma-ray emission as were subsequently obtained with INTEGRAL, with a SPI dataset comprising ∼200 Ms exposure time from 13.5 yr of data in the energy range from 1.795 to 1.820 MeV. We excluded observations with high rates of saturated detector events as well as 20% of the orbital phase around the perigee in order to avoid background from solar flares and passages through the Van Allen radiation belt. The observational sky coverage of the dataset in time and space is shown in Fig. 2. The patchy structure comes from the observation strategy of INTEGRAL, observing regions of interest separately, rather than performing a uniform full-sky survey.

Fig. 2. SPI exposure map for 13.5 yr of the INTEGRAL mission with integrated observation time of ∼200 Ms. The black box encloses the inner Galactic region, which we treated separately in our analysis. The effective area of SPI at 1.8 MeV is 44 cm2, accounting for dead detectors (Attie et al. 2003).

In order to spatially resolve this emission, SPI uses a coded-mask technique. It measures gamma-ray energies between 20 and 8000 keV with an angular resolution of 2.7°. According to the morphology of a source or extended emission feature and the orientation of the telescope, a characteristic shadow pattern is projected onto an array of 19 high purity Ge-detectors.

The celestial signal is overlaid by an instrumental background, originating from nuclear excitations of the instrument and spacecraft material itself. To determine the celestial gamma-ray signal from the shadow pattern, simple background subtraction would lead to erroneous results as we expected only a small number of counts per energy bin per detector per second and 90% background. Instead, a simultaneous fit for celestial and background signals has to be applied (Strong et al. 2005; Diehl et al. 2018; Siegert et al. 2019). We used the model description

(1)

for the measured instrument counts in each particular energy bin k for N_I sky and N_B background model components (Siegert et al. 2019). For each image element j, the celestial source model intensities M_ij were convolved with the image response function R_jk, i.e. the mask pattern. The parameters θ_i scale the intensities for all model components i. The background contributions B_ik are independent of the mask and the spacecraft orientation but vary on different time scales and over the entire mission due to the solar cycle, nuclear build-up processes, solar flares, and radiation belt transits. A detailed background model was obtained from a spectral fitting of the entire mission data (Diehl et al. 2018), which was then adapted to the particular dataset subjected to a specific analysis (Siegert et al. 2019). In this adaption, the background component has to be rescaled on an adequate time scale to properly represent such temporal variations. For the ²⁶Al analysis, we chose half-year intervals for the background normalisation, in addition to detector failure times (cf. Siegert et al. 2019, Sect. 5.1.2). Short-term variations were taken into account and modelled according to the saturated detector events tracking these at high statistical precision (e.g. Jean et al. 2003).

The thickness of the SPI anticoincidence shield, made up of 91 individual BGO crystals, is 5 cm. The attenuation cross section for photon energies around 1.8 MeV in BGO is 4.7 × 10⁻² cm² g⁻¹ (Berger et al. 2010). With a density of the crystals of 7.13 g cm⁻³, we estimated that the transmission probability for perpendicularly-incident photons is about 18%. This would lead to an additional signal mainly when the Galactic disk emission is coming from the side. The transmission probability drops below 1% for incidence angles larger than 69°; only 32% of all pointings are oriented towards these higher latitudes of |b|> 21°. Thus, the additional signal due to shield transparency can be considered small. Additionally, this component imprints more signal in outer detectors compared to inner ones. As the orientation of the spacecraft usually remains rather constant during one orbit, this would lead to a characteristic and quasi-constant background detector pattern on that timescale. Such celestial detector count contributions were generally included in our modelled background because we determined the background and its detector pattern per orbit.

In contrast, the pattern for counts from a sky signal captured within the coded mask varied according to the spacecraft orientation due to the coded mask pattern. This enabled us to distinguish background and source contributions.

As the measured detector counts are Poisson-distributed, the likelihood of a set of model parameters θ given a dataset D with n data points was calculated by the full Poisson likelihood

(2)

where d_k is the measured number of instrument counts and m_k is the model predicted value as described in Eq. (1). To determine the parameter set θ that maximises the likelihood, we used the negative logarithm of the Poisson likelihood dropping the data-dependent term, which is commonly referred to as Cash statistic (Cash 1979)

(3)

In our case, full-sky maps were taken as emission models, which were fitted to the data in detector space for each energy bin separately. In order to compare non-nested models, we employed a likelihood-ratio test, which will be described in Sect. 3.1.

We performed a spectral analysis of each fitted model to determine the 1.809 MeV line flux accurately above Galactic continuum emission, which contributes about 5% of the flux in the line band between 1805 and 1813 keV. We treated the spectral shape as a degraded Gaussian function which includes the effect of charge collection efficiency due to detector worsening over time (e.g. Kretschmer et al. 2013; Siegert 2017; Siegert et al. 2019). The average instrumental resolution at 1.8 MeV is 3.17 keV (Diehl et al. 2018).

2.2. Simulated maps

Fujimoto et al. (2018) present a high-resolution hydrodynamic simulation of a Milky Way-like spiral galaxy. They include the self-gravity of gas, fixed axisymmetric logarithmic potential to represent the gravity of old stars and dark matter, radiative cooling, photoelectric heating, as well as stellar feedback in the form of photoionisation, stellar winds, and SNe. The chemical enrichment of the galaxy was traced by following the dynamics of stellar ²⁶Al and ⁶⁰Fe ejecta, calculated via the Stochastically Lighting Up Galaxies (SLUG) population synthesis code (da Silva et al. 2012; Krumholz et al. 2015). Star-by-star yields of ²⁶Al were taken from Sukhbold et al. (2016). For their simulation, they assume an isolated gas disk orbiting in an otherwise static background potential representing dark matter and a stellar disk component.

The system evolved for t = 750 Myr with a maximum spatial resolution of 8 pc. After that running time, the physical and chemical structure of the ISM had reached a statistical equilibrium characterized by a steady large-scale structure of superbubbles filled with the freshly ejected nucleosynthesis products ²⁶Al and ⁶⁰Fe. These bubbles around massive star forming regions mainly followed the spiral arm pattern that developed spontaneously in the simulations, spatially exceeding the size of their host giant molecular clouds (GMCs). This is consistent with measurements from SPI determining the large scale gas dynamics of ²⁶Al (Kretschmer et al. 2013).

They derive full-sky flux maps for a hypothetical observer at the Solar Circle at r = 8 kpc by line-of-sight integration of the 1.8 MeV emission weighted with the distance squared. This was done for 36 different observer positions with 10° step size in the plane of the simulated galaxy. These maps cover the whole sky with a pixel size of 1° ×1°. An example of such an integrated map, for an observer at position 0°, is shown in Fig. 3.

Fig. 3. Simulated all-sky flux map of 1.8 MeV emission in Milky Way-like galaxy from hydrodynamic simulations by Fujimoto et al. (2018) plotted on same scale as COMPTEL map in Fig. 1.

3. Comparison of observations with simulations

3.1. Direct likelihood comparison

The most straightforward approach in comparing a simulation to observation is a direct evaluation of how well the simulation describes the data. We can expect, however, clear deviations due to the specific morphological features of the Milky Way, which are largely random results of the particular distribution of SN bubbles around the Sun, which we should not expect a simulation to match in detail. In this approach to comparison, we adopt the 36 flux maps obtained from the simulation by Fujimoto et al. (2018) as the emission models for the ²⁶Al sky. We then performed maximum likelihood fits to the observational datasets using proper instrument response and backgrounds (cf. Sect. 2.1) to determine the likelihood of the fitted simulation map.

Since our statistical method did not provide an absolute goodness of fit, we could not directly evaluate resulting likelihoods. As a criterion to rate the relative fit quality of different sky models, we applied the test statistic

(4)

which characterises a likelihood-ratio test of a sky model M₁ describing celestial emission on top of the background versus the null-model M₀ including only the background model. Thus, TS gives the likelihood of M₁ relative to the null-hypothesis of observing background only given the data D. With a sample of 1000 synthetic Monte Carlo datasets, we verified that TS is χ²/2-distributed and we can associate it with the probability of M₁ occurring by chance in D (cf. Appendix A).

We chose two different realms for this comparison: in the first case, we restricted our analysis to the inner galactic region with −10° < b <  10° and −30° < l <  30° (black box in Fig. 2) to treat the regions with highest intensity and longest exposure time separately. For the second case, we analysed the full sky emission. For each of the spatial realms, we determined the ²⁶Al signal in one energy band from 1805 keV to 1813 keV.

The model comparison is shown in Fig. 4. As observation-based reference points we include model comparisons with the full-sky maximum entropy map obtained from COMPTEL (Plüschke et al. 2001) and a map obtained from SPI measurements (Bouchet et al. 2015) as sky models. This gives test statistic values for maps representing observations of the actual Milky Way emission.

Fig. 4. Likelihood ratio of different sky maps relative to likelihood of a background-only fit (see Eq. (4)). This is obtained by fitting INTEGRAL data in the energy band 1805 keV–1813 keV for the inner Galaxy (solid lines) and the full sky (dashed lines). The simulated sky models (green) correspond to 36 observer positions in the simulation by Fujimoto et al. (2018). Values of TS for the COMPTEL (Plüschke et al. 2001) and SPI (Bouchet et al. 2015) maximum entropy maps are given as observation-based reference points. TS values for the COMPTEL and SPI map differ only by 2 (full sky) and 7 (inner galaxy) and overlap in the figure; therefore they are shown combined in single lines (blue/purple). Higher values of TS indicate a lower probability of occurrence by chance in our set of 13.5 years of SPI data (cf. Appendix A).

For the inner Galaxy, the simulated maps show values of TS mostly below the observation-based maps which are consistent with one another. Nevertheless, there are maps from sight-lines in the simulation which are as unlikely to occur by chance in the data as the observation-based maximum entropy reconstructions. This indicates that the observed structure in this region is, by and large, dominated by the overall Galactic morphology and therefore well-described by the generic hydrodynamic simulation. There are surprising variations among the different simulation samples, however. It is particularly striking that there are some observer positions where the simulations are actually less likely to be found by chance in the SPI data than the maximum entropy reconstruction. It is difficult to identify distinct morphological features, to which the striking improvement for some of the observer positions may be attributed. There is indication that this may be mainly due to the geometric configuration of superbubbles in the direction of the Galactic centre. We return to this issue in Sect. 3.2.1.

Taking the full sky into account, the observation-based maps show overall larger and also consistent values of the test statistic TS, while values for the simulated maps fall below, with a larger scatter than for the inner galaxy case. This indicates less-matching models (through a higher probability of chance coincidence for the simulated maps, compared to the observation-based maps). Thus, the simulations diverge significantly from the maximum entropy reconstruction over the full sky, particularly when we include higher latitudes.

Because the simulation primarily gives relative intensity variations on the sky, we fit the line flux to estimate the absolute changes. The obtained fluxes are given in Table 1 in comparison with the observational maximum entropy maps. For the isolated treatment of the inner Galaxy, the line fluxes derived from all maps match within the uncertainties. The flux values obtained from the full-sky analysis are significantly (∼7σ) lower for the simulated maps than for the COMPTEL map. Conversely, the SPI map gives a larger flux (∼4σ). As Fujimoto et al. (2018) state, prominent 1.8 MeV emission regions such as Cygnus, Carina, or Sco-Cen, as well as the characteristic spiral arm structure of the Miky Way, have been omitted in their simulation. Thus, it is reasonable that especially the characteristic foreground emission in the Milky Way, which originates from regions relatively close to the Earth and extends to higher latitudes, is missing. The COMPTEL map, as well as the SPI map, include such regions and features, thus providing a better fit to INTEGRAL/SPI data. Nevertheless, the good fitting results for the inner Galaxy indicate that this region is less influenced by characteristic foreground features but by the galaxy-wide emission.

3.2. Scale height analysis

Given our finding that numerical simulations fail to mimic the particular structure of the sky as seen from Earth, we now seek a more general method of comparing simulations and observations. Thus, we analyse their morphological features in a generalized way which can be applied to SPI measurements as well. Common morphological analyses like expansion in spherical harmonics or wavelets are very sensitive to the assumptions made for image reconstructions, for example, starting points of maximum entropy deconvolution. Thus, we stick to maximum likelihood estimations of chosen models, and evaluate characteristic distributions that can be inferred from those models. A representation of the sky in spherical harmonics basically carries an infinite number of realisations for each combination of degree and mode already as a result of a simple rotation on the sky. As we can only test each realisation of such an analytic model individually, this is not a basis for a viable approach. Thus, we have to choose a simpler model which contains basic morphological information and for which we can fit individual realisations to the data.

3.2.1. Galaxy-wide scale height and scale radius

Extragalactic studies show an exponential decrease of young massive stars with radius. As the 1.8 MeV emission from ²⁶Al traces such massive star groups, we expect it to follow a similar trend. Therefore, we assume a doubly exponential disk model

(5)

with the galactocentric radius r² = x² + y², the height above the disk z, scale radius r₀, scale height z₀, and amplitude of the disk A₀, as a first-order model for the galactic 3D distribution of the ²⁶Al emissivity in units of ph cm⁻³ s⁻¹. The 2D flux projection of this emission model onto the celestial sphere in galactic longitude l and latitude b is obtained by line-of-sight integration

(6)

from the relative position of the Sun with respect to the Galactic centre p_s = (x_s, y_s, z_s) = (8.5, 0, 0) kpc along the line-of-sight vector u = (u_x, u_y, u_z) = (cos(l)cos(b), sin(l)cos(b), sin(b)). The integration boundaries s_min and s_max confine the emission model to the volume of the Milky Way. Earlier studies of the inner Galaxy with SPI find a ²⁶Al scale height between 60 and 250 pc (Diehl et al. 2006; Wang et al. 2009) for a fixed scale radius of 4 kpc. COMPTEL all-sky datasets also suggested a scale height of 150−170 pc (Oberlack 1997; Diehl et al. 1997). In order to obtain a galaxy-wide evaluation of the ²⁶Al emission with SPI, we fit a grid of 32 × 64 combinations of different r₀ and z₀ to our set of 13.5 yr of SPI data. The scale radius r₀ ranges from 0.50 to 8.25 kpc in 250 pc steps. The scale height z₀ ranges from 10 to 475 pc in 15 pc steps and from 500 to 2050 pc in 50 pc steps (Siegert 2017). From the resulting Poisson likelihood values we calculate the probability density distribution for combinations of scale height and scale radius of the galaxy-wide ²⁶Al emission. The results are shown in Fig. 5. From the likelihood profiles of both parameters, we find a best-fitting exponential disk model with r₀ = 5.81 ± 0.64 kpc and z₀ = 0.77 ± 0.17 kpc for the Milky Way emission. The smaller scale height of 180 pc found with SPI by Wang et al. (2009) is due to a spatial restriction to the inner Galaxy.

Fig. 5. Fit results for 26Al data described by exponential disk model. The colour-coded surface gives the probability distribution function (PDF) of scale height and scale radius. The contour plot shows the 1, 2, and 3σ uncertainty contours around the best fit values (yellow dot) and the shaded histograms give the likelihood profile for scale height (green) and scale radius (blue).

We fit the same grid of exponential disk emission models to the 36 simulated ²⁶Al flux maps. This enables a comparison of the galaxy-wide morphology of simulations and gamma-ray measurements. We evaluate the best-fitting amplitude for each exponential disk model and calculate the difference between each simulated map and all exponential disk maps. In order to capture the same dynamical range as in SPI observations, we multiply the simulated sky maps with the SPI exposure map. This transfers the sky maps into count space projected onto the sky with effective area and sky coverage of SPI taken into account. We perform a Pearson’s χ² fit to estimate the best fitting amplitudes of all exponential disk models for each simulated map, treating the maps from the simulation as synthetic data and the exponential disk maps as model prediction. From the minimum χ² values we find the mean values of r₀ = 3.02 kpc and z₀ = 0.07 kpc for all synthetic sky maps. While the overall scale radius in the simulation is close to what we observe in the Milky Way, the scale height appears to be one order of magnitude smaller. However, we find one outlier at an observer position of 100° in the simulated galaxy, which shows the maximum overall scale height of z₀ = 0.7 kpc in agreement with observations. In the simulation, this is a unique spot where the observer is placed directly inside a superbubble of ∼1 kpc in size around two high star formation clumps, located in the direction of the galactic centre. This could indicate that, from a nucleosynthesis point of view, the Sun inside the Local Bubble is located at a rather exceptional location in the Galaxy.

3.2.2. Scale height frequency spectrum

In order to spatially resolve how the overall scale height is composed of certain features with different latitude extent, we investigate separate rectangular regions of interest (ROIs) of 12° in longitude and 180° in latitude each. The width in longitude was chosen to achieve a compromise between spatial resolution and intensity of the ²⁶Al signal per bin (Kretschmer et al. 2013). We also apply synthetic noise reflecting the dynamical range of the observations by a weighting with the SPI exposure as it was done in Sect. 3.2.1. For simplicity, we assume a fixed scale radius of r₀ = 5.5 kpc for the exponential disk models and retain only z₀ as free parameter. We determine the fit quality for a set of 123 full-sky exponential disk models with different scale heights between 10 pc and 5 kpc for each ROI separately. The step size is 15 pc from 10 to 475 pc and 50 pc from 0.5 to 5 kpc. The scale height of the best-fitting model in each longitude bin gives a measure of the extent of ²⁶Al emission in latitude. Since SPI has a field of view of 16 × 16 deg², two adjacent ROIs are overlapping. To account for the variability arising from the arbitrary placement of bins and their overlap, the evaluation was done for 12 different sets of bins shifted in 1° longitude steps.

This allows for a generalized comparison of the variances present in synthetic and observed maps of the galactic ²⁶Al emission that is particularly independent of the specific morphology seen in the Milky Way. The scale height frequencies derived from the simulation and SPI data are shown in Fig. 6. From SPI measurements we see a major contribution in small scale height features of the order of z₀ ∼ 10 pc. Additionally, a lower and fairly constant contribution of larger latitude extents is measured with a slight increase even up to a few kpc. This indicates the presence of two quite distinct emission features rather than a uniform exponential morphology with a consistent scale height. This implies that the actual Milky Way emission deviates substantially from an exponential disk model.

Fig. 6. Frequency spectrum of best-fitting scale heights in 12° longitude bins. The green line shows the median of all 36 flux maps from the simulation by Fujimoto et al. (2018), of which 68% lie inside the shaded green region. The purple line and shaded region give the median and the 68th percentile of SPI observations for 12 sets of longitude bins shifted by 1° each.

The simulation shows a similar excess for small-scale heights and a second component of scale heights around z₀ ∼ 100 pc. Small scale heights indicate dominant galactic disk emission and large-scale heights correspond to dominant nearby regions extending towards higher latitudes or even local emission around the Sun. While the frequency of small-scale heights in the simulation matches well with observations, the nearby large scale height component seems to be under-represented. This finding is consistent with the results from the direct comparison in Sect. 3.1.

3.2.3. Scale height vs. longitude

To see how the different scale height components are linked to the galactic morphology, we break down the results according to galactic longitude. This is depicted in Fig. 7. The uncertainties in scale height represent the average width of the log(ℒ)-profile in each longitude bin. Due to overlapping ROIs neighbouring data points are not strictly independent and grouped together alternately. The error bars in longitude account for the overlap.

Fig. 7. Best fitting scale height for individual 12° longitude bins along galactic plane. The grey lines show the results for all 36 observer positions in the simulation by Fujimoto et al. (2018) individually. Their median is given by the green line with the 68th percentile as shaded green region. The purple crosses represent the maximum likelihood scale heights obtained from SPI data with errors in y-direction depicting the average width of the log(ℒ)-profile in each longitude bin. Since adjacent ROIs are overlapping, data points are shown alternately with bold and thin markers to group mostly-independent data points together. The shaded region follows only the bold data points to guide the eye for a possible trend. The fits are restricted to 5 kpc scale heights.

In regions that are located away from the galactic centre at −60° > l >  60°, SPI observations are dominated by emission with small scale heights of ∼10 pc with a few contributions at intermediate large scale heights. On average, this general structure is also seen in the simulation, indicating that these regions in the Milky Way are shaped according to the stochastic galaxy-wide star formation processes. On the other hand, around the galactic centre at −60° < l <  60°, the measurements show intermediate-to-very-large-scale heights of the order of kpc. The flux maps from the simulation show overall smaller-scale heights between 10 pc and 100 pc in this region. In contrast to the previously considered outer regions, this implies a strong contribution of characteristic nearby foreground emission in the direction around l = 0°, e.g. from the Sco-Cen region at l = 350° (Krause et al. 2018). The structural difference could indicate that the run of star formation rate density with galactocentric radius differs in the simulation from the true conditions in the Milky Way, meaning that there may be more star formation locally than the simulation assumes. Another possible reason could be that the clustering of star formation in the simulation could differ from that of the observation. The clustering of superbubbles affects the vertical flow of nucleosynthesis material, with bigger superbubbles causing stronger vertical flows (von Glasow et al. 2013). Additionally, the spiral arms in the Milky Way could be attributed more to classical density waves than the spontaneously-formed structures in the simulation, related to external perturbations or the bar, perhaps. In reality, this may concentrate star formation in fewer, more prominent spiral arms.

4. Conclusions

In this paper we present methodological routes for making comparisons between observed maps of ²⁶Al decay in the gamma-ray sky and hydrodynamic simulations of a galactic ²⁶Al distribution. A direct comparison of a generic simulation with observations necessarily yields relatively poor fits because the simulation does not match particular prominent foreground features of the sky as seen from Earth. Nevertheless, we find that the portion of the sky around the Milky Way centre seems to be dominated by the overall galactic emission, rather than characteristic foreground structures, and that the simulation provides a reasonably good match to the observed gamma-ray sky in this region. This provides information about the 3D distribution of the ²⁶Al emission.

For a more generalized approach in terms of morphology, we investigate the latitude extent of the 1.8 MeV emission as parametrized by the exponential scale height. The rather broad emission with z₀ = 0.77 ± 0.17 kpc from SPI observations is only seen in one configuration in the simulation. In this case, the observer is placed inside a large superbubble structure filled with freshly produced ²⁶Al. This may imply that the Sun is placed in a similarly exceptional location in the Milky Way.

In order to spatially resolve certain emission features, we evaluated the scale height extent in 12° longitude bins. Analysis of SPI data reveals an almost bimodal distribution, with most longitudes showing a small scale height ≲50 pc, but with a small number showing extremely large scale heights of a few kpc. These large-scale height bins are mainly seen around the direction towards the Galactic centre. The simulation lacks such features on average. Such a large latitude extent of ²⁶Al must be associated with nearby regions. Thus, we confirm that the 1.8 MeV emission contains a significant contribution from nearby superbubbles, that is, from Sco-Cen or regions along the spiral arm tangents (e.g. del Rio et al. 1996; Diehl et al. 2010; Krause et al. 2018). We find indications that one of the most-nearby ²⁶Al-filled superbubbles may have overrun the Solar System already and, thus, they contribute an omnidirectional emission component (cf. Krause et al. 2018). This characteristic seen in the observations, and its relative infrequency in simulated sky maps, suggests that the Milky Way may have a more coarse-grained superbubble structure than that which has been modelled in the simulation. Indeed, Fujimoto et al. (2019) find that the pre-supernovae stage feedback implemented in Fujimoto et al. (2018) is not strong enough to disperse the surrounding gas completely, leaving star formation tracer emission too strongly associated with molecular gas tracer emission, which is inconsistent with observations of nearby galaxies. Thus, the 1.8 MeV map contains important information about the detailed geometry of massive star feedback in the Milky Way. The bimodal distribution of scale heights apparent in SPI measurements along different longitudes indicates that Galactic 1.8 MeV emission deviates significantly from a simple exponential disk model. This implies that previously determined scale heights to describe the 1.8 MeV emission of the entire Galaxy do not correspond to a physical scale height of ²⁶Al, and they are, instead, likely to be significantly biased by local foregrounds. Thus, an exponential disk model is insufficient to fit the local scale heights and a purely phenomenological model should be used for this kind of analysis in further studies. This adds systematic uncertainties to the current total ²⁶Al mass estimate in the Milky Way, because it is based on the assumption of a consistent Galactic scale height (Diehl et al. 2006).

Acknowledgments

We thank the anonymous referee for the constructive feedback and helpful suggestions on this work. We thank J. Michael Burgess for many useful discussions and comments on the statistics part. This research was supported by the Deutsches Zentrum für Luft- und Raumfahrt (DLR). Thomas Siegert is supported by the German Research Society (DFG-Forschungsstipedium SI 2502/1-1). The INTEGRAL/SPI project has been completed under the responsibility and leadership of CNES; we are grateful to ASI, CEA, CNES, DLR, ESA, INTA, NASA and OSTC for support of this ESA space science mission. This research was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government.

References

Appendix A: Distribution of TS

In order to judge the values of the likelihood ratio TS (cf. Eq. (4)) given in Fig. 4 with respect to the absolute level of comparison, we investigated its distribution. We constructed 1000 synthetic Poisson datasets containing background only. Fitting a sky model M₁ as well as a background only null-model M₀ to those datasets gives TS for occurrence of M₁ by chance due to Poisson fluctuations in each dataset. We fitted the COMPTEL map (Plüschke et al. 2001), the SPI map (Bouchet et al. 2015), as well as six cases of outstanding maps by Fujimoto et al. (2018, observer positions 10°, 70°, 110°, 120°, 250°, and 260°). The results are displayed in Fig. A.1, with the Fujimoto et al. (2018) results combined into one average curve. We confirm that TS is -distributed for all the maps. The factor 2 is due to the positivity of the Poisson distribution, whereby the negative half of symmetric normal distribution is left out (e.g. Mattox et al. 1996). Thus, we can associate TS with the probability for occurrence by chance of a sky model M₁ in our dataset, with an upper limit of 10⁻³ and evaluate absolute values of the test statistic TS.

Fig. A.1. Distribution of likelihood-ratio test statistic for simulated sky maps by Fujimoto et al. (2018), COMPTEL map (Plüschke et al. 2001), and SPI map (Bouchet et al. 2015) according to Eq. (4) using 1000 Poisson datasets sampled from background model only.

All Tables

All Figures

	Fig. 1. Maximum entropy all-sky map of 1.8 MeV emission in the Milky Way obtained from 9 yr of observations with COMPTEL (Plüschke et al. 2001).
In the text

	Fig. 2. SPI exposure map for 13.5 yr of the INTEGRAL mission with integrated observation time of ∼200 Ms. The black box encloses the inner Galactic region, which we treated separately in our analysis. The effective area of SPI at 1.8 MeV is 44 cm2, accounting for dead detectors (Attie et al. 2003).
In the text

	Fig. 3. Simulated all-sky flux map of 1.8 MeV emission in Milky Way-like galaxy from hydrodynamic simulations by Fujimoto et al. (2018) plotted on same scale as COMPTEL map in Fig. 1.
In the text

Fig. 4. Likelihood ratio of different sky maps relative to likelihood of a background-only fit (see Eq. (4)). This is obtained by fitting INTEGRAL data in the energy band 1805 keV–1813 keV for the inner Galaxy (solid lines) and the full sky (dashed lines). The simulated sky models (green) correspond to 36 observer positions in the simulation by Fujimoto et al. (2018). Values of TS for the COMPTEL (Plüschke et al. 2001) and SPI (Bouchet et al. 2015) maximum entropy maps are given as observation-based reference points. TS values for the COMPTEL and SPI map differ only by 2 (full sky) and 7 (inner galaxy) and overlap in the figure; therefore they are shown combined in single lines (blue/purple). Higher values of TS indicate a lower probability of occurrence by chance in our set of 13.5 years of SPI data (cf. Appendix A).

In the text

	Fig. 5. Fit results for 26Al data described by exponential disk model. The colour-coded surface gives the probability distribution function (PDF) of scale height and scale radius. The contour plot shows the 1, 2, and 3σ uncertainty contours around the best fit values (yellow dot) and the shaded histograms give the likelihood profile for scale height (green) and scale radius (blue).
In the text

	Fig. 6. Frequency spectrum of best-fitting scale heights in 12° longitude bins. The green line shows the median of all 36 flux maps from the simulation by Fujimoto et al. (2018), of which 68% lie inside the shaded green region. The purple line and shaded region give the median and the 68th percentile of SPI observations for 12 sets of longitude bins shifted by 1° each.
In the text

Fig. 7. Best fitting scale height for individual 12° longitude bins along galactic plane. The grey lines show the results for all 36 observer positions in the simulation by Fujimoto et al. (2018) individually. Their median is given by the green line with the 68th percentile as shaded green region. The purple crosses represent the maximum likelihood scale heights obtained from SPI data with errors in y-direction depicting the average width of the log(ℒ)-profile in each longitude bin. Since adjacent ROIs are overlapping, data points are shown alternately with bold and thin markers to group mostly-independent data points together. The shaded region follows only the bold data points to guide the eye for a possible trend. The fits are restricted to 5 kpc scale heights.

In the text

	Fig. A.1. Distribution of likelihood-ratio test statistic for simulated sky maps by Fujimoto et al. (2018), COMPTEL map (Plüschke et al. 2001), and SPI map (Bouchet et al. 2015) according to Eq. (4) using 1000 Poisson datasets sampled from background model only.
In the text