Galactic Population Synthesis of Radioactive Nucleosynthesis Ejecta

Diffuse gamma-ray line emission traces freshly produced radioisotopes in the interstellar gas, providing a unique perspective on the entire Galactic cycle of matter from nucleosynthesis in massive stars to their ejection and mixing in the interstellar medium. We aim at constructing a model of nucleosynthesis ejecta on galactic scale which is specifically tailored to complement the physically most important and empirically accessible features of gamma-ray measurements in the MeV range, in particular for decay gamma-rays such as $^{26}$Al, $^{60}$Fe or $^{44}$Ti. Based on properties of massive star groups, we developed a Population Synthesis Code which can instantiate galaxy models quickly and based on many different parameter configurations, such as the star formation rate, density profiles, or stellar evolution models. As a result, we obtain model maps of nucleosynthesis ejecta in the Galaxy which incorporate the population synthesis calculations of individual massive star groups. Based on a variety of stellar evolution models, supernova explodabilities, and density distributions, we find that the measured $^{26}$Al distribution from INTEGRAL/SPI can be explained by a Galaxy-wide population synthesis model with a star formation rate of $4$-$8\,\mathrm{M_{\odot}\,yr^{-1}}$ and a spiral-arm dominated density profile with a scale height of at least 700 pc. Our model requires that most massive stars indeed undergo a supernova explosion. This corresponds to a supernova rate in the Milky Way of $1.8$-$2.8$ per century, with quasi-persistent $^{26}$Al and $^{60}$Fe masses of $1.2$-$2.4\,\mathrm{M_{\odot}}$ and $1$-$6\,\mathrm{M_{\odot}}$, respectively. Comparing the simulated morphologies to SPI data suggests that a frequent merging of superbubbles may take place in the Galaxy, and that an unknown but strong foreground emission at 1.8 MeV could be present.


Introduction
The presence of radioisotopes in interstellar gas shows that the Milky Way is continuously evolving in its composition.This aspect of Galactic evolution proceeds in a cycle of star formation, nucleosynthesis feedback, and large-scale mixing.Radioisotopes trace this entire process through their γ-ray imprints, that is the production, ejection, and distribution of radionuclides.Thus, the investigation of radioactivity in the Galaxy provides an astrophysical key to the interlinkage of these fundamental processes.
The currently most frequently and thoroughly studied γ-ray tracers of nucleosynthesis feedback are the emission lines of 26 Al and 60 Fe ejecta from stellar winds and supernovae (SNe).Their half-life times of 0.7 Myr (Norris et al. 1983) and 2.6 Myr (Rugel et al. 2009), respectively, are comparable to the dynamic timescales of superbubble structures around stellar groups (de Avillez & Breitschwerdt 2005;Keller et al. 2014Keller et al. , 2016)).The nuclear decay radiation carries a direct signature of the physical connection between their production in massive stars and their distribution in the interstellar medium (ISM).Spatial mapping of the interstellar 26 Al emission at 1809 keV (e.g.Diehl et al. 1995;Oberlack et al. 1996;Plüschke et al. 2001;A&A 672, A54 (2023) with two major downsides: On the one hand, comparisons with tracer maps at other wavelengths include many astrophysical assumptions, such as the ionisation of massive stars that produce 26 Al that finally lead to free-free emission (Knödlseder 1999), which might put these comparisons on shaky grounds.On the other hand, descriptive models, such as doubly exponential discs, contain hardly any astrophysical input because only extents are determined.Heuristic comparisons therefore offer only limited potential for an astrophysical interpretation.These earlier studies focussed more on the general description of the overall γ-ray morphology because it was (and still is) unclear where exactly the emission originates.Prantzos & Diehl (1996) first discussed the different contributions to the Galactic 26 Al signal, distinguishing between massive stars, asymptotic giant branch (AGB) stars, classical novae, and cosmic-ray spallation.The ratio of 26 Al production from these earlier estimates is 1 : 2 : 3 for AGB, massive stars, and nova contributions in units of M ⊙ Myr −1 with a negligible cosmogenic production.From comparisons to COMPTEL data, Knödlseder (1999) found that the nova and AGB contribution may both only be up to 0.2 M ⊙ Myr −1 , whereas between 80-90% originate from massive stars.Recent modelling of classical novae by Bennett et al. (2013), taking into account the resonance interaction 25 Al(p, γ) 26 Si, suggest a nova contribution to the total 26 Al mass in the Milky Way of up to 30%.There is no resolved measurement of the 26 Al profile around Wolf-Rayet stars (Mowlavi & Meynet 2006), for example, and only a few measurements of the distribution inside superbubbles (Diehl et al. 2003;Krause et al. 2018), which can gauge the stellar evolution models.Owing to the instruments' capabilities with a typical angular resolution of 3 • , this is understandable.However, with more than 17 yr of data, the descriptive parameters, such as the total 26 Al line flux or the scale height of the disc (e.g.Pleintinger et al. 2019), are now precisely determined, which allows us to ask more fundamental questions.
In this study, we attempt to shift the focus from interpretations of descriptive parameters to astrophysical input parameters to describe the 26 Al sky.Modelling an entire galaxy for comparisons to γ-ray data appears intractable because hydrodynamics simulations suggest interdependencies that could hardly be modelled: First in-depth simulations regarding the Galaxywide distribution of nucleosynthesis ejecta have been performed by Fujimoto et al. (2018Fujimoto et al. ( , 2020) ) and Rodgers-Lee et al. (2019).These simulations started from the description of basic physical conditions, such as the gravitational potential of a galaxy and its temperature and density profile, to solve the hydrodynamics equations, and they followed the spread of freshly synthesised nuclei in a simulated galaxy.Comparing such simulations directly to γ-ray data has yielded insights into the astrophysical 3D modalities and dynamics underlying the measured radioactivity in the Milky Way (Pleintinger et al. 2019).A scientific exploitation of such simulations, however, faces major limitations: On the one hand, only a few realisations can be created because they are computationally expensive.Due to their particular characteristics (generalising the Galaxy, rather than accounting for the location of the Sun and nearby spiral arms, for example), these are only comparable to the Milky Way in a limited extent.On the other hand, current γ-ray instruments have poor sensitivities, so that the level of detail of hydrodynamic simulations cannot be covered by observations.
In this paper, we present an alternative approach to modelling the radioactive Galaxy that is specifically adapted to the empirical basis of γ-ray measurements.The essential scientific requirements here are that it can be repeated quickly and in many different parameter configurations and that, at the same time, the astrophysically most important and observationally accessible features can be addressed.For this purpose, we developed the Population SYnthesis COde (PSYCO), the general structure of which we outline in Sect. 2. It is designed as a bottom-up model of nucleosynthesis ejecta in the Galaxy, and based on population synthesis calculations of massive star groups, which is described in Sect.3. Our simulation results are shown in Sect. 4 and quantitative comparisons to the entire Galaxy are outlined in Sect.5, in particular summarising the astrophysical parameters of interest that can be obtained with γ-ray observations from the Spectrometer aboard the INTErnational Gamma-Ray Astrophysics Laboratory, INTEGRAL/SPI (Winkler et al. 2003;Vedrenne et al. 2003).We discuss our findings in Sect.6 and summarise in Sect.7.

Model structure and input parameters
The overall model structure of PSYCO is shown in Fig. 1 (Pleintinger 2020).While the input parameters are determined top-down from the galactic level to single star properties, the nucleosynthesis aspect is subsequently modelled bottom-up with population synthesis calculations of massive star groups (Sect.3).

Stellar mass and timescales
The overall SFR of the Milky Way remains debated, as it can be determined in many different ways and their results are almost as varied.For example, interpretations of HII regions obtain a value of (1.3 ± 0.2) M ⊙ yr −1 (Murray & Rahman 2010), of infrared measurements a value of 2.7 M ⊙ yr −1 (Misiriotis et al. 2006), of γ rays a value of 4 M ⊙ yr −1 (Diehl et al. 2006), or of HI gas a value of 8.25 M ⊙ yr −1 (Kennicutt & Evans 2012).Metaanalyses seem to settle on a value between 1.5 and 2.3 M ⊙ yr −1 (Chomiuk & Povich 2011;Licquia & Newman 2015).Likewise, a hydrodynamics-based approach by Rodgers-Lee et al. ( 2019) finds a range of 1.5-4.5 M ⊙ yr −1 .Because the SFR is an important model input variable, we use it as a free parameter.
The necessary model time is determined by the time scale of the longest radioactive decay (here 60 Fe).We start from an 'empty' galaxy which we gradually fill with mass in stars and radioactive ejecta.After the initial rise, owing to the constant SFR, a quasi-constant production rate balances a quasi-constant decay rate in the simulated galaxy.This leads to specific γ-ray line luminosities which can be compared to data.We exploit the fact that the galactic amount of species i with lifetime τ i will approach this equilibrium between production and decay after some time which is related to the SN rate which, in turn, is related to the stellar models used.In total, after T tot ≈ 50 Myr the galactic amount of the longer-lived 60 Fe has decoupled from the initial conditions and reached a quasi-equilibrium.This defines the total model time and accordingly the total mass M gal = SFR × T tot , which is processed into stars during that time.Treating the constant SFR as free parameter sets the main physical boundary on the stellar mass formed in the model.

Spatial characteristics
In order to ensure fast calculation, we use the following assumptions about the overall morphology and the metallicity gradient   .This is short enough that the position of the observer as well as gas and stars can be assumed as approximately co-rotating with respect to the spiral arms.Thus, a static galactic morphology divided in a radial and a vertical component is adopted.In the Galactic coordinate frame the Sun is positioned 8.5 kpc from the Galactic centre and 10 pc above the plane (Reid et al. 2016;Siegert 2019).We note that the Galactic centre distance may be uncertain by up to 8%, for which the resulting flux values and in turn luminosity, mass, and SFR may be uncertain by 14%.The vertical star formation density is chosen to follow an exponential distribution ρ(z; z 0 ) = z −1 0 exp(−|z|/z 0 ), with the scale height z 0 parametrising the Galactic disc thickness.The radial star formation density of the Galaxy can be approximated by a truncated Gaussian (Yusifov & Küçük 2004), with a maximum at radius R µ and width σ.Alternatively, an exponential profile can be chosen to obtain more centrally concentrated radial morphologies.The radial distribution is convolved with a 2D structure of four logarithmic spirals, approximating the expected spiral arm structure of the Milky Way.Each spiral centroid is defined by its rotation angle (Wainscoat et al. 1992), along the radial variable R, with inner radius R 0 , offset angle θ 0 , and widening k.In order to match the observed Milky Way spiral structure (Cordes & Lazio 2002;Vallée 2008) the values in Table 1 are adopted.Around the spiral centroids, a Gaussianshaped spread for star formation is chosen.In this work, it is generally assumed that the spatial distribution of star groups follows the Galactic-wide density distribution (Heitsch et al. 2008;Micic et al. 2013;Gong & Ostriker 2015;Krause et al. 2018).
Figure 2 depicts five different galactic morphologies that are implemented in the model.From pronounced spiral structures  to a smooth central peak, their nomenclature follows GM00-GM04, with GM00-GM02 using a Gaussian and GM03 and GM04 using an exponential density profile in radius.In the latter case, a scale radius of 5.5 kpc is adopted.In particular, GM00 describes a balance of a central density peak (for example from the central molecular zone) with the spiral arms.GM01 ignores the central peak, as found in Kretschmer et al. (2013) and modelled by Krause et al. (2015).GM02 enhances the peak central density by placing R µ of the Gaussian density closer to the centre.GM03 describes the distribution of pulsars (Faucher-Giguère & Kaspi 2006), and shows a large exponentially peaked maximum in the centre.Finally, GM04 represents an exponential radial profile, being completely agnostic about spiral features.The latter is typically used for inferences with γ-ray telescopes to determine the scale radius and height of the Galaxy in radioactivities.GM03 and GM04 could also mimic recent, several Myr year old, star burst activity near the Galactic centre as suggested from interpretations of the Fermi bubbles (e.g.Crocker 2012).We model the metallicity gradient in the Milky Way because it is measured to generally decrease radially (e.g.Cheng et al. 2012).
For the modelling of nucleosynthesis ejecta, this is an important factor because it determines the amount of seed nuclei and the opacity of stellar gas.To include this effect on nuclear yields and stellar wind strength, the Galactic metallicity gradients for different heights are included in the model following the measurements of Cheng et al. (2012).For different heights above the plane it can be approximated linearly with the parameters found in Table 2.

Distribution function of star formation events
Stars form in more extended groups as well as more concentrated clusters.Our model does not require information about the spatial distribution on cluster scale, so that we define 'star formation events'.Following Krause et al. (2015) we assume a single distribution function for all kinds of star-forming events, ξ EC , with the mass of the embedded cluster, M EC , which is empirically described by the power law with α EC = 2 in the Milky Way (Lada & Lada 2003;Kroupa et al. 2013;Krumholz et al. 2019).This applies to star formation events between 5 ≤ M EC / M ⊙ ≤ 10 7 (Pflamm-Altenburg et al. 2013;Yan et al. 2017).In order to approximate the physical properties of cluster formation, this relation is implemented in the model as probability distribution to describe the stochastic transfer of gas mass into star groups over the course of T tot = 50 Myr.

Initial mass function
Inside a star-forming cluster, the number of stars dN * that form in the mass interval between M and M + dM is empirically described by the stellar initial mass function (IMF) of the general form with power-law index α for stellar mass M * in units of M ⊙ and a normalisation constant k.The IMF itself is not directly measurable.Thus its determination relies on a variety of basic assumptions and it is subject to observational uncertainties and biases (Kroupa et al. 2013).Typically used IMFs are from Salpeter (1955, S55), Kroupa (2001, K01), andChabrier (2003, C03), which we include in our model.Details about the functional forms on these IMFs are found in Appendix A. These three variations of the IMF are implemented to describe the statistical manifestation of the physical star formation process in each star group.In all cases, we limit the IMF to masses below 150 M ⊙ since observationally higher stellar masses are, so far, unknown (Weidner & Kroupa 2004;Oey & Clarke 2005;Maíz-Apellániz 2008).We note that the yield models (Sect.2.3.3)typically only calculate up to 120 M ⊙ , so that extrapolations to higher masses are required, which might be unphysical.Given that the high stellar mass IMF index x 3 for stars above 1 M ⊙ may also be uncertain (x 3 ≈ 3 ± 1; Kroupa & Jerabkova 2021), but the fraction of stars above 100 M ⊙ is at most 0.05%, we consider the estimates robust against changes in the IMF above 100 M ⊙ .This mass limit, M max , seems to vary with the total cluster mass M EC (Pflamm-Altenburg et al. 2007;Kroupa et al. 2013): This relation can be approximated by log(M max ) = 2.56 log(M EC ) 3.
with λ = 9.17.The upper-mass limit is applied for each star group together with the lower-mass limit for deuterium burning of 0.012 M ⊙ (Luhman 2012).

Superbubbles
Radiative, thermal, and kinetic feedback mechanisms from individual stars in associations shape the surrounding ISM.This cumulative feedback produces superbubble structures around each stellar group.These structures are approximated as isotropically and spherically expanding bubbles (Castor et al. 1975;Krause & Diehl 2014).The crossing-time of stellar ejecta in such superbubbles is typically ∼ 1 Myr (Lada & Lada 2003), which corresponds to the lifetime of 26 Al.Due to this coincidence, the distribution of 26 Al is expected to follow the general dynamics of superbubbles (Krause et al. 2015).Therefore, the size scale of superbubbles is adopted as a basis for the spatial modelling of nucleosynthesis ejecta in our model.
Hydrodynamic ISM simulations by de Avillez & Breitschwerdt (2005), for example, show that after a similar timespan of ∼ 1 Myr the hot interior of a superbubble is homogenised.Thus, we model nucleosynthesis ejecta as homogeneously filled and expanding spheres with radius (Weaver et al. 1977;Kavanagh 2020) as a function of the mechanical luminosity L W and time t.The mechanical luminosity was modelled in population synthesis calculations of Plüschke (2001) with a canonical value of L W = 10 38 erg s −1 , which we adopt in this work.The free parameter x in Eq. ( 6) relates the ISM density, ρ 0 , and the constant for radial growth over time, α by x = αρ −1/5 0 .Literature values for α range between 0.51 for enhanced cooling in mixing regions (Krause et al. 2014) (see also Fierlinger et al. 2012;Fierlinger 2014;Krause et al. 2013Krause et al. , 2015) ) to 0.76 for the self-similar, analytical, adiabatic expanding solution (Weaver et al. 1977).We use a constant x of 4 × 10 3 kg −1/5 m 3/5 , which would correspond to a particle density around 100 cm −3 for the self-similar solution and a density around 20 cm −3 for enhanced cooling.This sets the boundary conditions for the temporal evolution of the individual superbubbles.Relevant bubbles sizes, that is when 26 Al and 60 Fe can still be found inside before decay, are therefore of the order of a few 100 pc.We discuss the possible biases in resulting fluxes and morphologies in Sect.6.1, and proceed with the assumptions presented above.
We assume cluster formation to be a stochastic process that occurs with a constant rate over time.Spatially, it is preferentially triggered, when gas is swept up by the gravitational potential of a spiral arm.Each newly formed stellar cluster is assigned a 3D position according to a galactic morphology as shown in Fig. 2. The 3D position then implies the metallicty as calculated from the Galactic metallicity gradient (Sect.2.1.2).In the model, the fundamental parameters of total mass M EC , age T bubble , position (x, y, z) bubble , and metallicity Z bubble are assigned to each individual star group.

Stellar rotation
Stellar rotation creates additional advection, turbulent diffusion, and enhances convective regions inside the star (Endal & Sofia 1978;Heger et al. 2000), which increases mixing and transport of material inside the star.Stellar winds are amplified and also occur earlier due to the rotation.The wind phase as well as the entire lifetime of the stars is also extended.Stellar evolution models suggest that these effects have a significant impact on nucleosynthesis processes inside the stars as well as on their ejection yields (e.g.Limongi & Chieffi 2018;Prantzos et al. 2020;Banerjee et al. 2019;Choplin & Hirschi 2020).We implement stellar rotation for our galactic nucleosynthesis model by the following considerations: For each star that forms in the model, an additional step is included in the population synthesis process to randomly sample a rotation velocity according to measured distributions.Stellar rotation properties have been catalogued observationally for each spectral class by Glebocki et al. (2000); Glebocki & Gnacinski (2005).To include this information in population synthesis calculations, the observed rotation velocities are weighted with the average inclination angle of the stars.The resulting distributions of rotation velocities are then fitted for each spectral class individually by a Gaussian on top of a linear tail.In this context, the most relevant classes are the massive O-and B-type stars.They show a maximum at 100 km s −1 with a width of 60 km s −1 (O) and 0 km s −1 and width 180 km s −1 (B), respectively.

Explodability
At the end of the evolution of massive stars, a lot of processed material is ejected in SNe.However, this only applies if the stellar collapse is actually followed by an explosion, which is parametrised as 'explodability'.Explosions can be prevented under certain circumstances if the star collapses directly into a black hole instead.Depending on the complex pre-SN evolution, this naturally has a strong impact on nucleosynthesis ejecta.Different simulation approaches by Smartt (2009, S09), Janka et al. (2012, J12), Sukhbold et al. (2016, S+16), or Limongi & Chieffi (2018, LC18) provide strikingly different explodabilities.Figure 3 shows effects on 26 Al and 60 Fe ejection, respectively, over the entire stellar mass range with nuclear yield calculations by Limongi & Chieffi (2006, LC06).While the SN yields cease with suppression of the explosion, the wind ejecta remain unaffected by explodability.Because 60 Fe is ejected only in SNe but 26 Al also in winds, the 60 Fe/ 26 Al ratio is an important tracer of explodability effects on chemical enrichment.The explodability can thus be chosen in the model as input parameter in order to test their astrophysical impact on large-scale effects of nucleosynthesis ejecta.

Nucleosynthesis yields
The most fundamental input to the modelling of nucleosynthesis ejecta is the total mass produced by each star over its lifetime.This yield depends on many stellar factors, such as rotation, mixing, wind strength, and metallicity, as described above.It also involves detailed nuclear physics, which is represented in the nuclear reaction networks of the stellar evolution models.
Detailed yield calculations for 26 Al and 60 Fe in particular have been performed, for example, by Meynet et al. (1997), Limongi & Chieffi (2006), Woosley & Heger (2007), Ekström et al. (2012), Nomoto et al. (2013), Chieffi & Limongi (2013), or Limongi & Chieffi (2018).A comparison of the models is shown in Fig. 4.While the 26 Al yield is generally dominated by SN ejection in the lower-mass range, it is outweighed by wind ejection in the higher-mass regime.The 60 Fe wind yields are negligible.On average, a massive star ejects of the order of 10 −4 M ⊙ of 26 Al as well as 60 Fe.The predictions show a spread of about one to two orders of magnitude between models in the stellar mass range of 10-120 M ⊙ .Stars of lower mass, such as AGB stars with about 4 M ⊙ are expected to eject much less 26 Al around ∼5 × 10 −5 M ⊙ (Bazan et al. 1993).
For the overall contribution to nucleosynthesis feedback it is important to take the formation frequency of stars in a certain mass range into account.A convolution with the IMF shows that stars that form with a zero-age main sequence (ZAMS) mass of ≤30 M ⊙ contribute overall more 26 Al to the ISM than more massive stars.Due to the sensitivity of 60 Fe to explodability, the contribution to its overall amount by stars with M ZAMS ≥ 25 M ⊙ is negligible if SNe are not occurring.In our galactic nucleosynthesis model, we use the stellar evolution models by Limongi & Chieffi (2006, LC06) and Limongi & Chieffi (2018, LC18) A54, page 6 of 19 as they include 26 Al and 60 Fe production and ejection in timeresolved evolutionary tracks over the entire lifetime of the stars, thus avoiding the need for extensive extrapolations.We thus obtain population synthesis models of star groups that properly reflect the underlying physical feedback properties and their time variability.Other yield models can also be included, if a large enough grid in the required model parameters is available, so that interpolations can rely on a densely and regularly spaced input.We note again that individual yields of radioactive isotopes are quite uncertain as they depend on extrapolations from laboratory experiments towards different energy ranges.For example, Jones et al. (2019) found that the cross section 59 Fe(n, γ) 60 Fe has a linear impact on the 60 Fe yields and suggest a smaller than previously accepted value to match γ-ray observations.

Stellar binaries
The impact of binary star evolution, especially in terms of nucleosynthesis, is a heavily debated field.Roche-lobe overflows and tidal interactions can change the composition of binary stars as a whole and also enhance the ejection of material.The effects of binary star evolution are highly complex because of the unknown influence of the many stellar evolution parameters.First binary yield calculations for 26 Al have been provided by Brinkman et al. (2019).We perform a quantitative check of binarity impacts on the stellar group scale by including these yields in PSYCO.
For a canonical stellar group of 10 4 M ⊙ , we assume an overall binary fractions of 70% (Sana et al. 2012;Renzo et al. 2019;Belokurov et al. 2020), as well as extreme values of 0% and 90%.If a star has a companion or not is sampled randomly according to the selected fraction as they emerge from the same IMF.Brinkman et al. (2019) restricted the stellar evolution to that of the primary star so that we treat the companion as a single star.In addition, Brinkman et al. (2019) only considers wind yields so that we have to assume SN ejecta to follow other models.We inter-and extrapolate the parameter grid from Brinkman et al. (2019), including orbital periods, masses, and orbit separations to a similar grid as described above.The resulting population synthesis for a single stellar group is shown in Fig. 5.The two extreme explodability assumptions S09 and LC18 are shown for comparison.Independent of the SN model choice, the effects from binary wind yields appear rather marginal.In particular for low-mass stars, being the dominant 26 Al producers after 10 Myr, considering binaries has no impact.The reduction of wind yields in binary systems with large separation and primary stars of 25-30 M ⊙ leads to less 26 Al ejection after 15 Myr.The increased binary wind yields for stars with ≲20 M ⊙ results in a slightly enhanced ejection after 17 Myr.In addition, at early stages, the ejection from very massive stars tends to be reduced due to binary interactions.It is important to note here that these extra-and interpolations come with large uncertainties so that the binary yield considerations here should not be overinterpreted.Given the wind yield models by Brinkman et al. (2019), the variations for a canonical star group all lie within the 68th percentile of a single star population synthesis results.

Population synthesis
The cumulative outputs of the PSYCO model is built up step by step using the input parameters outlined in Sect. 2 as depicted on the left in Fig. 1 from the yields of single stars to properties of massive star groups to the entire Galaxy.The underlying method is a population synthesis approach, which relates the integrated signal of a composite system with the evolutionary properties of its constituents (Tinsley 1968(Tinsley , 1972;;Cerviño 2013).

Star group
Individual sources of interstellar radioactivity remain unresolved by current γ-ray instruments.However, integrated cumulative signals from stellar associations can be observed (e.g.Oberlack et al. 1995;Knoedlseder et al. 1996;Kretschmer et al. 2000;Knödlseder 2000;Martin et al. 2008Martin et al. , 2010;;Diehl et al. 2010;Voss et al. 2010Voss et al. , 2012)).This describes star groups as the fundamental scale on which the Galactic model of nucleosynthesis ejecta is based upon.
In massive star population synthesis, time profiles ψ(M * , t) of stellar properties (e.g.ejecta masses, UV brightness) are integrated over the entire mass range of single stars, weighted with the IMF ξ(M * ), with the normalisation A according to the total mass of the population (group or cluster).It has been shown that continuous integration of the IMF within population synthesis calculations can result in a considerable bias for cluster properties such as its luminosity (Piskunov et al. 2011).This is related to the discrete nature of the IMF, resulting in a counting experiment with seen and unseen objects of a larger population, which is difficult to properly treat without knowing the selection effects.In order to take stochastic effects in smaller populations into account, we therefore use a discrete population synthesis by Monte Carlo (MC) sampling (Cerviño & Luridiana 2006).Publicly accessible population synthesis codes have been developed and applied successfully to astrophysical questions (e.g.Popescu & Hanson 2009;da Silva et al. 2012).A generic population synthesis code including selection effects, biases, and overarching distributions has recently been developed by Burgess & Capel (2021).
Focussing on the specific case of radionuclei and interstellar γ-ray emission, we base our population synthesis approach on Fig. 6.Population synthesis models of nucleosynthesis ejecta in a 10 4 M ⊙ star group based on yield models by Limongi & Chieffi (2006, LC06) or Limongi & Chieffi (2018, LC18) as indicated in subcaptions.Different panels illustrate the impact of the physical parameters (a) IMF, (b) metallicity, (c) explodability, and (d) stellar rotation.Each panel shows the time evolution of the ejection of 26 Al (upper), 60 Fe (middle), and their average mass ratio (bottom).In each case, coeval star formation is assumed and the shaded regions indicate the 68th percentiles derived from 10 3 Monte Carlo runs.
the work by Voss et al. (2009), who also included kinetic energy and UV luminosity evolution of superbubbles.
In a first step, individual initial mass values are sampled according to the IMF.In order to assure a discretisation and to reproduce the shape of the IMF as close as possible, we apply the optimal sampling technique, which was developed by Kroupa et al. (2013) and revised as well as laid out in detail by Schulz et al. (2015).Krause et al. (2015) have shown that the details of the sampling method do not influence 26 Al abundances and superbubble properties significantly.Details of the optimal sampling method is given in Appendix B. It is based on the total mass M EC of the cluster to be conserved, during the formation of stars with mass M min ≤ M * ≤ M max according to an IMF ξ(M * ).
In the next step, each single star is assigned a stellar rotation velocity according to its spectral class (Sect.2.3.1).In addition to the total mass, each star group is assigned a position drawn from the Galactic density distribution (Sect.2.1.2) and accordingly a metallicity (Table 2).
By adding the respective stellar isochrones, we thus obtain cumulative properties of star groups based on the described parameters and assuming coeval star formation for any given event/group.The original evolutionary tracks by Limongi & Chieffi (2006, 2018) cover only a few ZAMS masses and irregular time steps.For the population synthesis, a uniform and closely meshed (fine) grid of stellar masses in 0.1 M ⊙ steps and evolution times in 0.01 Myr is created by interpolations.We use linear extrapolation to include stellar masses above 120 M ⊙ and below 13 M ⊙ .Figure 6 shows the effect of the main physical input parameters on 26 Al and 60 Fe ejection with a population synthesis of a 10 4 M ⊙ cluster for different assumptions of the IMF, metallicity, explodability, and stellar rotation.They are mainly based on models by Limongi & Chieffi (2018), as they cover this whole range of physical parameters.Models by Limongi & Chieffi (2006)  effects, as they do not include intrinsic assumptions about this parameter.For better cross-comparison, stellar masses have been obtained by random sampling in this case to determine stochastic uncertainty regions.
The choice of the IMF (Fig. 6, top left) affects the amplitude of nucleosynthesis feedback for both 26 Al and 60 Fe.While K01 and C05 yield similar results within the statistical uncertainties, S55 shows a strong reduction of nucleosynthesis ejecta.This is readily understood because S55 continues unbroken towards low-mass stars which is unphysical.This distributes a great amount of mass into a large number of low-mass stars -that is, those that do not produce major amounts of 26 Al or 60 Fe.
The production of 26 Al and 60 Fe generally relies strongly on the initial presence of seed nuclei, mainly 25 Mg and 59 Fe, respectively.Thus, an overall metallicity reduction in the original star-forming gas ultimately also goes along with a decrease of ejection yields of 26 Al and 60 Fe (Fig. 6, top right) (Timmes et al. 1995;Limongi & Chieffi 2006).This effect is strong for 60 Fe because 56 Fe is produced only on a very short time scale in late evolutionary stages.Thus, only marginal amounts of processed 56 Fe can reach the He-or C-burning shells where 60 Fe production can occur (Tur et al. 2010;Uberseder et al. 2014).Reduced metallicity also decreases the opacity of stellar material (Limongi & Chieffi 2018).As a consequence, convective zones shrink and stellar winds decrease with lower radiative pressure.Both effects reduce 26 Al yields because it resides in hot regions where it is destroyed and the wind component ceases.
The impact of explodability is shown to be strongest for 60 Fe because this isotope is ejected only in SNe (Fig. 6, bottom left).The more extensive the inhibition of explosions, the stronger the reduction of 60 Fe yields, especially at early times (higher initial masses).This effect is comparably weak for 26 Al and accounts for only a factor of 2 less ejection from the most massive stars because the wind component remains unaffected.
Due to the increased centrifugal forces in fast rotating stars (Fig. 6, bottom right), the core pressure is reduced and its overall lifetime extended (Limongi & Chieffi 2018).Thus, nucleosynthesis feedback is stretched in time when stellar rotation is included.For 26 Al, this effect is mainly recognisable as a slightly earlier onset of winds.The changes in 60 Fe are much stronger.The enhancement of convection zones increases the neutron-rich C-and He-burning shells significantly and enhanced material from even deeper layers can be mixed into these regions.This leads to a boost in 60 Fe production by a factor of up to 25.Additionally, the duration of 60 Fe ejection in a star group is prolonged by a factor of about 2.
Figure 6 shows that the mass ratio 60 Fe/ 26 Al is particularly sensitive to changes in metallicity, explodability and rotational velocity.A one order of magnitude reduction in metallicity leads to a decrease of the same magnitude in this mass ratio.The most significant impacts on the temporal behaviour have rotation effects.Due to the prolongation of stellar evolution and the drastic increase in 60 Fe ejection, 60 Fe/ 26 Al is dominated by 60 Fe after only ∼ 10 Myr.If stellar rotation is not taken into account 26 Al dominance lasts for ∼ 18 Myr.Changes in explodability show also a shift in this time profile.If explosions of massive stars with M * > 25 M ⊙ are excluded, 60 Fe domination is delayed by ∼13 Myr to about 16 Myr after cluster formation.This underlines that the 60 Fe/ 26 Al ratio is an important observational parameter that provides crucial information about detailed stellar physics.
Due to the reproducibility and uniqueness of optimal sampling, time profiles of star groups can be calculated in advance.We take advantage of this fact and create a database that covers a broad parameter space of combinations of cluster masses, explodability models, yield models and IMF shapes.This procedure drastically reduces the computing time of the overall galaxy model (Sect.3.2).

Galaxy
We extend the population synthesis to the galactic level by calculating a total galactic mass that is processed into stars with a constant star formation rate (SFR) over 50 Myr as described in Sect.2.1.1.Because the embedded cluster mass function (ECMF) behaves similarly to the IMF, the optimal sampling approach is also used here.In addition to the variables of mass and time, the spatial dimension is added at this level.Star groups form at different positions and times in the Galaxy and their spatial extents evolve individually over time (Eq.( 6)).In order to transfer this 3D information into an all-sky map that can be compared to actual γ-ray measurements, we perform a line-of-sight integration (for details, see Appendix C).
The ejecta of an isotope n with atomic mass m n,u distributed homogeneously in a bubble with radius R SB (t) at time t emit γ rays at energy E γ from a nuclear transition with probability p E γ with a luminosity which is normalised to a unit Solar mass of isotope n and is expressed in units of ph s −1 M ⊙ −1 .For example in the case of 26 Al, the 1.8 MeV luminosity per solar mass is 1.4 × 10 42 ph s The amount of isotope n present in the superbubble is predetermined for each point in time by the massive star population synthesis.This determines an isotope density ρ n which is constant inside and zero outside the bubble.Burrows et al. (1993) and Diehl et al. (2004) suggest that ejecta remain 'inside' the bubble.There could be some mixing of 26 Al with the HI walls.Hydrodynamic simulations find a varying degree of concentration of nucleosynthesis ejecta in the supershells (Breitschwerdt et al. 2016;Krause et al. 2018).In any case, the superbubble crossing time of about 1 Myr again would make the ejecta appear homogenised.Given the angular resolution of current γ-ray instruments of a few degrees, we therefore find that a homogeneous density inside the bubbles is a good first-order approximation (see Sect. 6.1 for a discussion).
By line-of-sight integration of these homogeneously filled spheres, a spatial γ-ray emission model for each superbubble is created and added onto the current model map.Their cumulative effect finally gives a complete galactic picture.This formulation is easily adaptable to arbitrary isotopes by scaling L n,⊙ and lifetimes τ n .In the case of short-lived isotopes such as 44 Ti with a lifetime of 89 yr, for example, the spatial modelling reduces to point sources for SPI because the ejecta do not travel far from their production site before decaying.

Evaluation of PSYCO models
We have evaluated a grid of models, varying our input parameters with SFR ∈ {2, 4, 8} M ⊙ yr −1 , scale height z 0 ∈ {0.1, 0.2, 0.3, 0.5, 0.7} kpc, density profiles GM00-GM04, the two stellar evolution models LC06 and LC18, and the explodabilities S09, and S+16 (and LC18 to match the LC18 stellar evolution model).We chose to use only the IMF K01.For each parameter value combination, 100 MC runs are performed to From this number of simulations, we can explore links (correlations) between the parameters and assign some uncertainty to those.Naturally, the SFR and explodability have an impact on the total amount of 26 Al and 60 Fe present in Galaxy: For LC18 (stellar evolution model and explodability), the total galactic 26 Al mass follows roughly a linear trend M 26 /M ⊙ ≈ 0.25 × SFR/(M ⊙ yr −1 ); for other explodabilities, the SFR impact is larger M 26 /M ⊙ ≈ 0.31-0.52× SFR/(M ⊙ yr −1 ).For 60 Fe, the effects of explodability are reversed since 60 Fe is only ejected in SNe.We find M 60 /M ⊙ ≈ 1.72 × SFR/(M ⊙ yr −1 ) for LC18, and M 60 /M ⊙ ≈ 0.28-1.27× SFR/(M ⊙ yr −1 ) for the other explodability models.The resulting mass ratio 60 Fe/ 26 Al has therefore almost no SFR-dependence, and we find 60 Fe/ 26 Al of 0.9 for LC18, up to 7.1 for LC06.We note that there are crucial differences in the flux, mass, and isotopic 60 Fe/ 26 Al ratio: Given that the γ-ray flux F n of an radioactive isotope n is proportional to M n p γ,n m −1 n τ −1 n (see Eq. ( 9)), the flux ratio of 60 Fe to 26 Al in the Galaxy as a whole is The conversions between flux ratio and isotopic ratio, mass ratio M 60 /M 26 , isotopic ratio N 60 /N 26 = M 60 /M 26 • m 26 m 60 , and production rate Ṁ60 / Ṁ26 = M 60 /M 26 • τ 26 τ 60 , are given in Table 3.The SN rates from these model configurations are directly proportional to the SFR, as expected, and follow the trend SN rate/century −1 ≈ 0.37-0.56× SFR/(M ⊙ yr −1 ), with the explodability LC18 giving the lowest SN rate and S09 the highest.The values above are independent of the chosen density profiles.By contrast, the 1.809 MeV ( 26 Al) and 1.173 and 1.332 MeV ( 60 Fe) fluxes are largely dependent on the chosen spiral-arm prominence.We show trends of these values as derived from PSYCO simulations in Appendix D.

Overall appearance
Figure 7 shows the convergence of 60 Fe and 26 Al masses in the model within T tot = 50 Myr.This is an artificial diagnostic of radioactive masses to reach a constant value in a steady state.The stellar evolution models LC06 and LC18 both lead to an equilibrium between production and decay within this time span owing to the constant star formation.A change in SFR alters the overall amplitude, leaving the general convergence behaviour unaffected.This means that after a modelling time of 50 Myr, the distribution of the isotopes 26 Al and 60 Fe is determined only by the assumed distribution of the star groups in space and time, and no longer by the initial conditions.Due to the specific scientific focus on nucleosynthesis ejecta, we therefore choose the snapshot at 50 Myr to evaluate all our models (see also Sect. 4.1).About 50% of the total 26 Al γ-ray flux is received from within 6 kpc in GM03, that is the most shallow profile.For the spiral-arm-dominated profiles, GM00 and GM01, half the flux is already contained within 2.8 kpc.It is important to note here that most of the flux received excludes the Galactic centre with a distance of 8.5 kpc.Until the distance of the Local Arm tangent at about 2 kpc, on average 30% of the flux is enclosed (see also Fig. 8).In addition, it is interesting that in 0.3% of all cases (i.e. 90 out of 30 000 models, see Sect.4.1), about 90% of the total flux comes from a region of only 6 kpc around the observer.This means that the local components outweigh the overall Galactic emission by far in these cases.More centrally dominated morphologies (especially GM04) show a flatter slope than spirally dominated ones (GM00-GM02).Our models show similar flux profiles compared to the simulation by Rodgers-Lee et al. (2019), for example.With respect to this hydro-simulation, the best agreement is found with the centrally weighted spiral morphology GM03.Particular observer positions in the hydrodynamic simulation can show a strong contribution from the Local Arm.Such a behaviour is reflected by spiral-arm-dominated morphologies such as GM01 in our model.Based on these general agreements between galacticwide population synthesis and hydrodynamics simulations, it is suggested that some important properties of the Galaxy can be transferred to the simpler-structured population synthesis model to test stochastic effects and a variety of parameters.

Galactic 26 Al and 60 Fe fluxes
The total flux as well as the flux distribution of these models, either across the celestial sphere or along different lines of sight, play a major role in the interpretation of the γ-ray signals.The total flux in γ-ray measurements is nearly independent of the chosen morphology (e.g. the SPI or COMPTEL maps, or an exponential disc lead to the same fluxes within 5%), so that the absolute measurements of F 26 = (1.71± 0.06) × 10 −3 ph cm −2 s −1 (Pleintinger et al. 2019) and F 60 = (0.31 ± 0.06) × 10 −3 ph cm −2 s −1 (Wang et al. 2020)
constraints.Indeed, the chosen density profile has a considerable impact on the total 26 Al and 60 Fe fluxes, changing by 50%, with GM03 typically showing the lowest fluxes and GM01 the highest.In addition to the total flux, the 'Inner Galaxy' (|l| ≤ 30 • , |b| ≤ 10 • ) is frequently used for comparisons in data analyses because in this range, the surface brightness of 26 Al (and supposedly 60 Fe) is particularly high.We summarise an evaluation of 26 Al and 60 Fe fluxes in Table 4 for both the entire sky and the Inner Galaxy.This is obtained by an average across density profiles and therefore includes an intrinsic scatter of 25%.
The scale height z of the density profile also influences the total fluxes, with smaller scale heights typically resulting in larger fluxes for otherwise identical parameter sets.The effect is stronger for 60 Fe than for 26 Al because a larger scale height leads to a larger average distance of sources to the observer (the galaxy is 'larger').For density profiles which show an enhancement closer to the observer (GM00, GM01, GM04), this effect is stronger than for more centrally peaked profiles.In addition, the later 60 Fe ejection compared to 26 Al preferentially fills older and larger bubbles.The 60 Fe emission is intrinsically more diffuse so that an additional vertical spread enhances the r −2 -dependence of the flux which amplifies the scale height effect for 60 Fe.We show the radial distributions of expected flux contributions for 26 Al and 60 Fe in Fig. 8. Bright regions indicate 'where' the most measured flux would originate in.Clearly in all profiles, the local flux contributions, and especially the spiral arms (GM00-GM03), shape the resulting images (Fig. 9).Even though there is a density enhancement in GM00, for example, the resulting image would appear devoid of such a feature because the r −2 effect lets the Galactic centre feature appear washed out.Interestingly, the exponential profiles (GM03, GM04) show both, a central flux enhancement as well as a local contribution, which might be closer to the expected profile of classical novae in the Milky Way plus massive star 26 Al emission.However, such a strong central enhancement is not seen in either the COMPTEL map nor the SPI map, even though exponential discs nicely fit the raw data of the two instruments.Interpretations are discussed in Sect.6.

Likelihood comparisons
MeV γ-ray telescopes cannot directly 'image' the sky -the typically shown all-sky maps are individual reconstructions (realisations) of a dataset projected back to the celestial sphere, assuming boundary conditions and an instrumental background model.As soon as those datasets change or another instrument with another aperture and exposure is used, also the reconstructions can vary significantly.Likewise, different realisations of PSYCO, even with the same input parameters of density profile, scale height, SFR, yield model, and explodability, will always look different due to the stochastic approach.Therefore, a comparison of individual models in 'image space' can -and will most of the time -be flawed.
In order to alleviate this problem, the comparisons should happen in the instruments' native data spaces.This means that any type of image is to be convolved with the imaging response functions and compared in the native data space.By taking into account an instrumental background model, then, the likelihoods of different images (all-sky maps) can be calculated and set in relation to each other.As an absolute reference point (likelihood maximum), we use the 26 Al all-sky maps from COMPTEL (Plüschke et al. 2001) and SPI (Bouchet et al. 2015).
Clearly, individual morphological features of the Milky Way may not be mapped in all realisations of PSYCO, which will  result in 'bad fits'.These discrepancies are expected as they are mostly dominated by random effects of the particular distribution of superbubbles in the Galaxy and in the MC simulations.One particular realisation is therefore not expected to match all (relevant) data structures.
An almost direct comparison is nevertheless possible to some extent: We use the INTEGRAL/SPI 26 Al dataset from Pleintinger et al. (2019) to investigate which of the large-scale parameters are required to maximise the likelihood.Ultimately, this results in an all-sky map (or many realisations thereof) which can be compared to SPI (and other) data.In particular, we describe our sky maps as standalone models M, which are converted to the SPI instrument space by applying the coded-mask pattern for each of the individual observations in the dataset (see Pleintinger et al. 2019, for details on the dataset as well as the general analysis method).An optimisation of the full model, including a flexible background model (Diehl et al. 2018;Siegert et al. 2019) and a scaling parameter for the image as a whole, results in a likelihood value for each model emission map.We note that this does not necessarily provide an absolute goodnessof-fit value.Instead, a relative measure of the fit quality can be evaluated by a test statistic, TS = 2(log(L (D|M 0 )) − log(L (D|M 1 ))), (11) with model M 1 describing the general case of an image plus a background model, and M 0 describing only the instrumental background model.The likelihood L (D|M 0 ) of the data D given the background model then describes the null-hypothesis, which is tested against the alternative L (D|M 1 ).TS values can hence be associated to the probability of occurrence by chance of a certain emission map in the SPI dataset.The higher the TS values are, the 'better fitting' an image therefore is.We show a summary of the TS values for GM00 (best cases) and GM04 (typically used model) in Fig. 10.Independent of the actual SFR, the fit quality is almost the same, showing that the SFR has no morphological impact between 2 and 8 M ⊙ yr −1 with GM00 resulting in a slightly higher SFR.This is understood because the amplitude of the emission model, that is one scaling parameter for the entire image, is optimised during the maximum likelihood fit.Only models with SFR ≥ 4 M ⊙ yr −1 show a fitted amplitude near 1.0 which means that SFRs of less than 4 M ⊙ yr −1 mostly underpredict the total 26 Al flux.The density profiles GM03 and GM04 (exponential profiles) show on average lower TS values than Gaussian profiles, independent of the chosen spiral structure.The reason for these apparently worse fits compared to GM00-GM02 is the central peakedness of the exponential profiles, which is absent in actual data.The highest average TS values are found with GM02, which can be described by a protruding and rather homogeneous emission from the Inner Galaxy (Fig. 9, middle), mainly originating from the nearest spiral arm.However, high-latitude emission is barely present in GM02 because the density maximum is about 2.5 kpc away from the observer.In general, a flatter gradient between the two nearest spiral arms describes the data better than a steep gradient.The scale height affects the quality of the global fit with a trend of a higher scale height generally fitting the data worse, except for the density profile GM00.A large-scale height increases the average distance to the supperbubbles, which decreases their apparent sizes and the general impact of foreground emission.However, since GM00 shows the largest individual TS values (i.e.fits the SPI data best among all combinations considered), which also improves further with larger scale heights, this indicates that SPI data includes strong contributions from high latitudes in the direction of the galactic anticentre.In fact, GM00 is the only density profile that shows improved fits with both increasing SFR and larger scale height.Compared to the other density profiles, GM00 requires a large SFR to explain the total flux and to fully develop its characteristics on the galactic scale because it is also the most radially extended profile reaching up to 15 kpc.
The individually best model is found for GM00 with a scale height of 700 pc, a SFR of 8 M ⊙ yr −1 , IMF K01, stellar evolution model LC06, and explodability L18, and shown in Fig. 11.With a TS value of 2061, it is about ∆TS = 100 away from the maximum of the SPI map itself with (TS = 2166).The characteristics of this map is a rather bimodal distribution, peaking towards the Inner Galaxy and the Galactic anticentre.The spiral arm 'gaps' (tangents) are clearly seen in this representation, and the map appears rather homogeneous (many bubbles overlapping to smear out hard gradients).The map in Fig. 11, however, would result in an 'unfair' comparison to the reconstructed maps of COMPTEL or SPI with an intrinsic resolution of 3 • , and a sensitivity of much more (i.e.worse) than the minimum depicted one of 10 −6 ph s −1 cm −2 sr −1 .We therefore convolve the image with a 2D-Gaussian of width 2 • , and set the minimum intensity of the maps in Fig. 12 to 5 × 10 −5 ph s −1 cm −2 sr −1 .Clearly, many features of the actual map disappear because the sensitivities of the instruments are not good enough so that especially high-latitude features beyond |b| ≳ 30 • would drown in the background.

Degeneracy between superbubble physics, yields, and star formation
The PSYCO model includes assumptions, which cannot be verified individually, yet they interplay to create the structures of the Fig. 12. Compilation of observational maps (top: COMPTEL; middle: SPI) compared to our best-fitting PSYCO simulation, adopted to match the instrument resolution of 3 • .The minimum intensity in the maps is set to 5 × 10 −5 ph s −1 cm −2 sr −1 to mimic potentially observable structures.
output maps.In particular, nucleosynthesis yields and SFRs both scale the total flux predicted by the model.Observable structures depend on the size adopted to be filled with nucleosynthesis ejecta: we use one location centred on the massive-star group as a whole, and distribute the cumulative ejecta from the group within a spherical radius determined by the group age.Projecting this onto the sky with the adopted distance of the group to the observer, a circular patch of 26 Al emission results on the map.We know that real cavities around massive stars are not spherical, but of irregular shapes (see the example of Orion-Eridanus Brown et al. 1995;Bally 2008, anddiscussion in Krause et al. 2013;Krause & Diehl 2014).Also, as discussed above, the windand SN-blown cavities may not be filled homogeneously with 26 Al and 60 Fe, respectively.Therefore, we caution that values of the SFR and SN rate discussed in the following should be viewed as dependent in detail on the validity of these assumptions.

Star formation and supernova rate
The fluxes measured by SPI in the 1.8 MeV line from 26 Al and the 1.3 MeV line from 60 Fe, when fitted to the PSYCO bottom-up model, lead to SFRs of ≳4 M ⊙ yr −1 when considering the entire Galaxy, and ∼ 5 M ⊙ yr −1 when considering the Inner Galaxy.These values are obtained for the stellar evolution models by Limongi & Chieffi (2006, 2018), explodability models of Smartt (2009); Janka (2012); Sukhbold et al. (2016), and the range of velocties, metallicities and scale heights tested.
Considering uncertainties and variations with such assumptions, the SN rate in the Milky Way is estimated to be 1.8-2.8 per century.This range is purely determined by systematic uncertainties from model calculations; the statistical uncertainties would be on the order of 2-4%.Furthermore, this SFR sets the total mass of 26 Al and 60 Fe in the Galaxy to be 1.2-2.4M ⊙ and 1-6 M ⊙ , respectively.

60 Fe/ 26 Al ratio
The galactic-wide mass ratio of 60 Fe/ 26 Al mainly depends on the stellar evolution model and the assigned explodabilities.We find a mass ratio of 2.3 for model & explodability combination LC06 & S09 and 7.1 for LC18 & LC18.This converts to expected flux ratios F 60 /F 26 of 0.29 and 0.86, respectively, independent of the region in the sky chosen, that is either full sky or Inner Galaxy.These values are higher compared to measured value of 0.18 ± 0.04 (Wang et al. 2020), which is mainly due to the larger measured 26 Al flux compared to the PSYCO values.Measurements of the 60 Fe γ-ray lines are currently difficult because radioactive buildup of 60 Co in SPI leads to an ever-increasing background in these lines (Wang et al. 2020).Therefore, the information value of the current Galactic 60 Fe/ 26 Al should not be over-interpreted.In fact, Wang et al. (2020) also estimated systematic uncertainties of the 60 Fe/ 26 Al flux ratio to reach up to 0.4, which would be consistent with a range of parameter combinations of PSYCO.In fact, an increased 60 Fe flux from measurements might also be possible considering current detection significances of the two 60 Fe lines combined of 5σ.Similar to the differences in Fig. 11, showing the total emission, and Fig. 12, showing only 'detectable' features, the generally higher flux of 26 Al throughout the Galaxy compared to 60 Fe might enhance the observed discrepancy even further.The total 26 Al flux is consistently underestimated in PSYCO compared to measurements.Either the total 26 Al mass is underestimated (also ejected per star, for example), or there is a general mismatch in the density profiles (spiral arms, foreground emission).Only some extreme model configurations, for example with SFR = 8 M ⊙ yr −1 and no explodability constraints (S09), can reach the high observed fluxes.The density profiles have a 50% impact on the total Galactic flux.The Inner Galaxy only contributes to 16% on average, but can be measured more accurately with SPI because of the enhanced exposure time along the Galactic plane and bulge.The Inner Galaxy contribution to the total sky varies in PSYCO between 23% for centrally peaked profiles and 40% for spiral-arm-dominated profiles.This can be interpreted as requiring an additional strong local component at high latitudes, or a particularly enhanced Local Arm towards the Galactic anticentre.A bright spiral arm component in the Milky Way is supported by a general latitudinal asymmetry from the fourth quadrant, (−90 • ≤ l ≤ 0 • ) (Kretschmer et al. 2013;Bouchet et al. 2015).For example, a more prominent Sagittarius-Carina arm could achieve such an asymmetry.
Comparing the scale sizes of PSYCO 26 Al with 60 Fe simulations (Fig. 9), it appears that scale radius and scale height should differ. 60Fe γ-ray emission rarely appears towards the Galactic anticentre, while it is even required in 26 Al emission.Wang et al. (2020)  −1.0 kpc and 0.8 +0.3 −0.2 kpc, respectively, it is clear that meaningful comparisons rely on the quality of 60 Fe measurements.In our simulations, about 50% of the total 26 Al line flux originates from within a radius between 2.8 and 6.0 kpc.The 60 Fe fluxes originate from an even smaller region up to only 4 kpc radius, with a tendency towards the Galactic centre, so that only the Inner Galaxy appears bright in 60 Fe emission.If strong foreground sources are present, this fraction can even be as high as 90% for 26 Al, however remains rather stable for 60 Fe.
Our density model GM00 best describes the SPI 26 Al data, also putting a large emphasis on the Local Arm as well as the next nearest arms.With a large-scale height of 700 pc also nearby sources can be mimicked -underlining the importance of foreground emission once more.

Foreground emission and superbubble merging
The explodability of LC18 generally describes the 1.8 MeV SPI data better than other combinations.This is related to a trend to fill superbubbles with the majority of 26 Al later.Thus, larger bubbles would appear brighter.The average diameter of 26 Al-filled superbubbles in the Galaxy is about 300 pc.The underestimation in models of both the 1.8 MeV local foreground components and the average bubble size indicates that the actual local star formation density in the vicinity of the Solar System is larger than the average, and that it occurs overall more clustered than currently assumed.Our PSYCO full-sky images, when compared to SPI and COMPTEL, might support this as they appear more structured than in the simulations by Fujimoto et al. (2018), for example.Clustered star formation releases energy more concentrated through stellar feedback processes.As a result, the average size of the superbubbles would be larger and such a mechanism could also account for the more salient granularity in the observed scale height distribution of the Milky Way (Pleintinger et al. 2019).An increased bubble size might point to frequent superbubble merging, as suggested by Krause et al. (2013Krause et al. ( , 2015) ) or Rodgers-Lee et al. (2019).Here, HI shells break up frequently and open up into previously blown cavities, which lets them grow larger as a consequence of the feedback contributions from multiple star groups.

Superbubble blowout and Galactic wind
Using simulations similar to Rodgers-Lee et al. (2019), Krause et al. (2021) have characterised the vertical blowout of superbubbles from the Galactic disc.In some of their simulations, the superbubbles tend to merge in the disc and create transsonic, 26 Al carrying outflows into the halo.With typical velocities of 1000 km s −1 and a half-life of 0.7 Myr, kpc scale heights can be expected.Given that massive stars are formed within the Milky Way disc within a typical scale height of about 100 pc or less (Reid et al. 2016), the blowout interpretation appears to be a likely explanation for the large-scale heights we find.
Rodgers-Lee et al. ( 2019) also show that the halo density constrains such vertical blowouts, and due to the higher halo density in the Galactic centre for a hydrostatic equilibrium model the A54, page 15 of 19 A&A 672, A54 (2023) scale height should be smallest there.The apparent bimodal scale height distribution in our best-fitting models (higher towards Galactic centre and anti-centre) might hence point to a significant temporal reduction of the halo density via the Fermi bubbles (e.g.Sofue 2000;Predehl et al. 2020;Yang et al. 2022).

Summary and outlook
We have established a population-synthesis-based bottom-up model for the appearance of the sky in γ-ray emission from radioactive ejecta of massive star groups.This is based on stellar-evolution models with their nucleosynthesis yields, as well as representations of the spatial distribution of sources in the Galaxy.Parameters allow for these components to be adjusted, and thus they provide direct feedback of varying model parameters on the appearance of the sky.We parametrised, specifically, the explodability of massive stars, the contributions from binary evolution, the density profile and spiral-arm structure of the Galaxy, and the overall SFR.PSYCO can be easily adapted to other galaxies and, for example, model radioactivity in the Large Magellanic Cloud.
Application of the PSYCO approach to Galactic 26 Al finds agreement of all major features of the observed sky.This suggests that, on the large scale, such a bottom-up model captures the sources of 26 Al on the large scale of the Galaxy with their ingredients.Yet, quantitatively, the PSYCO model as based on current best knowledge fails to reproduce the all-sky γ-ray flux as observed.This suggests that nucleosynthesis yields from current models may be underestimated.On the other hand, mismatches as to the detail appear particularly at higher latitudes, and indicate that nearby sources of 26 Al with their specific locations play a significant role in the real appearance of the sky, and also in the total flux observed from the sky.We know that several such associations, for example Cygnus OB2 (Martin et al. 2009), Scorpius Centaurus (Diehl et al. 2010), or Orion OB1 (Siegert & Diehl 2017), should be included for a more-realistic representation.We note, however, that here details of superbubble cavity morphologies will be important (Krause et al. 2013(Krause et al. , 2015)), and our spherical volume approximation, while adequate for more distant sources on average, will be inadequate.Such refinements, and inclusions of very nearby cavities from the Scorpius-Centaurus association and possibly the Local Bubble, are beyond the scope of this paper.
Measurements of the 26 Al emission will accumulate with the remaining INTEGRAL mission until 2029.The 60 Fe γ-ray lines, however, have become difficult to analyse because radioactive buildup of 60 Co in SPI leads to an ever-increasing background in these lines (Wang et al. 2020).With the COSI instrument (Tomsick et al. 2019) on a SMEX mission planned for launch in 2027, a one order of magnitude better sensitivity after 2 yr could be achieved, so that also weaker structures, as predicted from our PSYCO models, could be identified in both the 26 Al and 60 Fe lines.Also, the Large Magellanic Cloud with an expected 1.8 MeV flux of 2 × 10 −6 ph cm −2 s −1 would be within reach for the COSI mission.

Fig. 1 .
Fig. 1.Structure of the PSYCO model.Using the model input on the right, model parameters were accumulated from the top down.The nucleosynthesis aspect on the left was finally built from the bottom up to construct all-sky γ-ray maps.

Fig. 2 .
Fig. 2. Radial probability density distributions for star formation.The morphologies GM00-GM04 are ordered from top to bottom by decreasing dominance of the Galactic spiral structure.GM04 is a purely exponential disc.

Fig. 3 .
Fig.3. 26Al (blue) and Fe (green) yields from SNe byLimongi & Chieffi (2006) for different explodability models.Stars with an initial mass inside the grey shaded regions eject no material during the SN.Islands of explodability following each other closely appear as green regions.
Fig. 4. Nucleosynthesis yields (top 26 Al, bottom 60 Fe) for a selection of stellar evolution models, separated into the stellar wind and total contribution (wind + SN).

Fig. 5 .
Fig.5.Population synthesis a canonical star group of 10 4 M ⊙ including binary effects.Shown is the impact of 26 Al binary wind yields fromBrinkman et al. (2019, B+19).Binary systems are included with orbital periods between 3 and 100 days and with an overall fraction of 70% (dashed line) and 90% (dotted line), which is contrasted with an association of only single stars (solid line).Binaries with longer periods, i.e. wider separations are considered to evolve as single stars.Lines indicate the average of 1000 Monte Carlo runs each.
Difference when including rotation.

Fig. 7 .
Fig. 7. Steady-state settling time of the total 60 Fe (blue) and 26 Al (green) in PSYCO galaxy models.Shaded regions denote the 68th percentile of 100 MC model runs.All models are based on evolutionary tracks LC06 (dashed lines) or LC18 (solid lines) and explodability models S09 and LC18, respectively, for SFR = 4 M ⊙ yr −1 and the K01 IMF.

Fig. 8 .
Fig. 8. Radial distribution of modelled flux contributions for a theoretical observer (white circle) from 26 Al (top) and 60 Fe (bottom) decay.Each column represents average results of 500 model instantiations based on the density profiles GM00-GM04 (grey boxes).The models shown are stellar evolution models LC06, explodability S09, IMF K01, and SFR = 4 M ⊙ yr −1 .The latter corresponds to a total mass of 1.8 ± 0.2 M ⊙ and 4.2 ± 0.2 M ⊙ of 26 Al and 60 Fe, respectively.Adaptive spatial binning (Cappellari & Copin 2003) is used to obtain Voronoi tessellations as spatial bins, each of which contribute a flux of 10 −6 ph cm −2 s −1 for the observer.The colour scale refers to the collecting area covered by each such pixel.

Table 4 .
Fluxes of 26 Al emission at 1.809 MeV and of 60 Fe emission at 1.173 or 1.332 MeV, respectively, in units of 10 −4 ph cm −2 s −1 from PYSCO simulations for the entire sky or the Inner Galaxy (|l| ≤ 30 • , |b| ≤ 10 • ), as a function of SFR (in units of M ⊙ yr −1 ) and different stellar evolution models.

Fig. 9 .
Fig. 9. Simulated full-sky γ-ray maps of the 1.8 MeV line from 26 Al decay (left) and the 1.3 MeV line from 60 Fe decay (right) constructed with PSYCO.Each row represents an individual MC run based on a different density profile (grey boxes) with a scale height of 300 pc.Nucleosynthesis yields are based on LC06 with the S09 explodability and the K01 IMF.

Fig. 10 .
Fig. 10.Likelihood ratio of 6000 sky maps modelled with PSYCO relative to the likelihood of a background-only fit with the SPI.Dots and solid lines denote the average values from 100 MC runs as a function of scale height.The colours correspond to stellar model configuration as noted in the legend.Triangles mark the maximum TS value obtained from the 100 MC samples for each model configuration.The thick grey lines denote the reference value obtained with COMPTEL (TS = 2160) and SPI (TS = 2166).

Fig
Fig. D.1: Supernova rate as a function of star formation rate for different explodability assumptions.

Fig
Fig. D.2: 60 Fe/ 26 Al mass ratio as a function of star formation rate for different explodability assumptions.

Table 2 .
(Cheng et al. 2012for modelling the radial metallicity gradient in the Milky Way according to different heights above the Galactic plane(Cheng et al. 2012).
Notes.The Intersect is in units of the metallicity, [Fe/H], and the Slope in units of [Fe/H] kpc −1 .