Highlight
Open Access
Issue
A&A
Volume 672, April 2023
Article Number A54
Number of page(s) 19
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202244457
Published online 30 March 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 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 26Al and 60Fe 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. 2014, 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 26Al emission at 1809 keV (e.g. Diehl et al. 1995; Oberlack et al. 1996; Plüschke et al. 2001; Bouchet et al. 2015) and detailed γ-ray line spectroscopy (e.g. Diehl et al. 2006, 2010; Kretschmer et al. 2013; Siegert & Diehl 2017; Krause et al. 2018) provided observational insights into the relation between nucleosynthesis ejecta and the dynamics of massive star groups. Together with the detection of 60Fe emission lines from the ISM at 1173 keV and 1332 keV (Wang et al. 2007, 2020), this represents an indispensable astrophysical effort to understand the feedback cycle underlying Galactic chemical enrichment.

In order to follow and understand the complex dynamics, astrophysical models and simulations have to be utilised. Because directly mapping the diffuse emission of the 1.8 MeV line is difficult and image reconstructions can come with considerable bias, given the photon-counting nature of the measurements, empirical and descriptive models of 26Al in the Galaxy have been invoked in the past (e.g. Prantzos 1993; Prantzos & Diehl 1995; Knoedlseder et al. 1996; Lentz et al. 1999; Sturner 2001; Drimmel 2002; Alexis et al. 2014). Their scientific interpretations rely on comparisons between the measured 26Al emission and multi-wavelength ‘tracers’ or geometric (smooth) emission morphologies (e.g. Hartmann 1994; Prantzos & Diehl 1995; Diehl et al. 1997, 2004; Knödlseder et al. 1999; Kretschmer et al. 2013). Such heuristic approaches can come 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 26Al 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 26Al signal, distinguishing between massive stars, asymptotic giant branch (AGB) stars, classical novae, and cosmic-ray spal-lation. The ratio of 26Al 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 25Al(p, γ)26Si, suggest a nova contribution to the total26Al mass in the Milky Way of up to 30%. There is no resolved measurement of the 26Al 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 26Al 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 26Al 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 Galaxy-wide distribution of nucleosynthesis ejecta have been performed by Fujimoto et al. (2018, 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.

2 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).

2.1 Galactic scale

2.1.1 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 y 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 60Fe). 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 Ttot ≈ 50 Myr the galactic amount of the longer-lived 60Fe has decoupled from the initial conditions and reached a quasi-equilibrium. This defines the total model time and accordingly the total mass Mgal = SFR × Ttot, 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.

2.1.2 Spatial characteristics

In order to ensure fast calculation, we use the following assumptions about the overall morphology and the metallicity gradient at the Galactic level: During the expected visibility of an individual star group of ~ 25 Myr as defined by the decay of the 60Fe ejected from the last SN (see Sect. 3.1), the Milky Way rotates by about 38°. 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 , with the scale height z0 parametrising the Galactic disc thickness. The radial star formation density of the Galaxy can be approximated by a truncated Gaussian (Yusifov & Küçük 2004), (1)

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), (2)

along the radial variable R, with inner radius R0, 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 Gaussian-shaped 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.

thumbnail 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.

Table 1

Spiral centroid parameters (cf. Eq. (2)) following Faucher-Giguère & Kaspi (2006).

Table 2

Linear parameters for modelling the radial metallicity gradient in the Milky Way according to different heights above the Galactic plane (Cheng et al. 2012).

2.2 Stellar groups

2.2.1 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, MEC, which is empirically described by the power law (3)

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 ≤ MEC/M ≤ 107 (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 Ttot = 50 Myr.

thumbnail 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.

2.2.2 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 (4)

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), and Chabrier (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 x3 for stars above 1 M may also be uncertain (x3 ≈ 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, Mmax, seems to vary with the total cluster mass MEC (Pflamm-Altenburg et al. 2007; Kroupa et al. 2013): This relation can be approximated by (5)

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).

2.2.3 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 isotrop-ically 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 26Al. Due to this coincidence, the distribution of 26Al 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) (6)

as a function of the mechanical luminosity LW and time t. The mechanical luminosity was modelled in population synthesis calculations of Plüschke (2001) with a canonical value of LW = 1038erg 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 . Literature values for a 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. 2013, 2015) to 0.76 for the self-similar, analytical, adiabatic expanding solution (Weaver et al. 1977). We use a constant x of 4 × 103 kg−1/5 m3/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 26Al and 60Fe 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 MEC, age Tbubble, position (x, y, z)bubble, and metallicity Zbubble are assigned to each individual star group.

2.3 Stellar parameters

2.3.1 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.

thumbnail Fig. 3

26AI (blue) and 60Fe (green) yields from SNe by Limongi & 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.

2.3.2 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 26Al and 60Fe 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 60Fe is ejected only in SNe but 26Al also in winds, the 60Fe/26Al 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.

2.3.3 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 26Al and 60Fe 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 26Al yield is generally dominated by SN ejection in the lower-mass range, it is outweighed by wind ejection in the higher-mass regime. The 60Fe wind yields are negligible. On average, a massive star ejects of the order of 10−4 M of 26Al as well as 60Fe. 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 26Al 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 26Al to the ISM than more massive stars. Due to the sensitivity of 60Fe to explodability, the contribution to its overall amount by stars with Mzams ≥ 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) as they include 26Al and 60Fe production and ejection in time-resolved 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 59Fe(n, γ)60Fe has a linear impact on the 60Fe yields and suggest a smaller than previously accepted value to match γ-ray observations.

thumbnail Fig. 4

Nucleosynthesis yields (top 26Al, bottom 60Fe) for a selection of stellar evolution models, separated into the stellar wind and total contribution (wind + SN).

thumbnail Fig. 5

Population synthesis a canonical star group of 104 M including binary effects. Shown is the impact of 26Al binary wind yields from Brinkman 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.

2.4 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 26Al 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 104 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 26Al 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 26Al 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 overinter-preted. 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.

3 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, 1972; Cerviño 2013).

3.1 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. 2008, 2010; Diehl et al. 2010; Voss et al. 2010, 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*), (7)

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 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 26Al abundances and superbubble properties significantly. Details of the optimal sampling method is given in Appendix B. It is based on the total mass MEC of the cluster to be conserved, (8)

during the formation of stars with mass MminM*Mmax 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 26Al and 60Fe ejection with a population synthesis of a 104 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) are chosen to show explodability 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 26Al and 60Fe. 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 26Al or 60Fe.

The production of 26Al and 60Fe generally relies strongly on the initial presence of seed nuclei, mainly 25Mg and 59Fe, respectively. Thus, an overall metallicity reduction in the original star-forming gas ultimately also goes along with a decrease of ejection yields of 26Al and 60Fe (Fig. 6, top right) (Timmes et al. 1995; Limongi & Chieffi 2006). This effect is strong for 60Fe because 56Fe is produced only on a very short time scale in late evolutionary stages. Thus, only marginal amounts of processed 56Fe can reach the He- or C-burning shells where 60Fe 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 26Al 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 60Fe because this isotope is ejected only in SNe (Fig. 6, bottom left). The more extensive the inhibition of explosions, the stronger the reduction of 60Fe yields, especially at early times (higher initial masses). This effect is comparably weak for 26Al 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 26Al, this effect is mainly recognisable as a slightly earlier onset of winds. The changes in 60Fe 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 60Fe production by a factor of up to 25. Additionally, the duration of 60Fe ejection in a star group is prolonged by a factor of about 2.

Figure 6 shows that the mass ratio 60Fe/26Al 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 60Fe ejection, 60Fe/ 26Al is dominated by 60Fe after only ~ 10 Myr. If stellar rotation is not taken into account 26Al 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, 60Fe domination is delayed by ~ 13 Myr to about 16 Myr after cluster formation. This underlines that the 60Fe/26Al 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).

thumbnail Fig. 6

Population synthesis models of nucleosynthesis ejecta in a 104 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 26Al (upper), 60Fe (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 103 Monte Carlo runs.

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 mn,u distributed homogeneously in a bubble with radius RSB(t) at time t emit γ rays at energy Eγ from a nuclear transition with probability pEγ with a luminosity (9)

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 26Al, the 1.8 MeV luminosity per solar mass is . 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 pn 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 26Al 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 Ln,⊙ and lifetimes τn. In the case of short-lived isotopes such as 44Ti 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.

4 Simulation results

4.1 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 z0 ∊ {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 explod-abilities 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 estimate stochastic variations, which in total amounts to 30 000 simulated PSYCO maps.

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 26Al and 60Fe present in Galaxy: For LC18 (stellar evolution model and explodability), the total galactic 26Al mass follows roughly a linear trend M26/M ≈ 0.25 × SFR/(M yr−1); for other explodabilities, the SFR impact is larger M26/M ≈ 0.31–0.52 × SFR/(M yr−1). For 60Fe, the effects of explodability are reversed since 60Fe is only ejected in SNe. We find M60/M ≈ 1.72 × SFR/(M yr−1) for LC18, and M60/M ≈ 0.28–1.27 × SFR/(M yr−1) for the other explod-ability models. The resulting mass ratio 60Fe/26Al has therefore almost no SFR-dependence, and we find 60Fe/26Al of 0.9 for LC18, up to 7.1 for LC06. We note that there are crucial differences in the flux, mass, and isotopic 60Fe/26Al ratio: Given that the γ-ray flux Fn of an radioactive isotope n is proportional to (see Eq. (9)), the flux ratio of 60Fe to 26Al in the Galaxy as a whole is (10)

The conversions between flux ratio and isotopic ratio, mass ratio M60/M26, isotopic ratio N60/N26 = M60/M26 · , and production rate 60/26 = M60/M26 · , 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 explod-ability 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 (26Al) and 1.173 and 1.332 MeV (60Fe) fluxes are largely dependent on the chosen spiral-arm prominence. We show trends of these values as derived from PSYCO simulations in Appendix D.

Table 3

Conversions between 60Fe/26Al flux, mass, isotopic, and production ratios.

4.2 Overall appearance

Figure 7 shows the convergence of 60Fe and 26Al masses in the model within Ttot = 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 26Al and 60Fe 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 26Al γ-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 galactic-wide 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.

thumbnail Fig. 7

Steady-state settling time of the total 60Fe (blue) and 26Al (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.

5 Comparison to data

5.1 Galactic 26Al and 60Fe 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 F26 = (1.71 ± 0.06) × 10−3 ph cm−2 s−1 (Pleintinger et al. 2019) and F60 = (0.31 ± 0.06) × 10−3 ph cm−2 s−1 (Wang et al. 2020) are important model constraints. Indeed, the chosen density profile has a considerable impact on the total 26Al and 60Fe fluxes, changing by 50%, with GM03 typically showing the lowest fluxes and GM01 the highest. In addition to the total flux, the ‘Inner Galaxy’ (|| ≤ 30°, |b| ≤ 10°) is frequently used for comparisons in data analyses because in this range, the surface brightness of 26Al (and supposedly 60Fe) is particularly high. We summarise an evaluation of 26Al and 60Fe 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 60Fe than for 26Al 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 60Fe ejection compared to 26Al preferentially fills older and larger bubbles. The 60Fe 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 60Fe. We show the radial distributions of expected flux contributions for 26Al and 60Fe 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 26Al 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.

thumbnail Fig. 8

Radial distribution of modelled flux contributions for a theoretical observer (white circle) from 26Al (top) and 60Fe (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 26Al and 60Fe, 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 26Al emission at 1.809 MeV and of 60Fe 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.

5.2 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 26Al 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 26Al 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 goodness-of-fit value. Instead, a relative measure of the fit quality can be evaluated by a test statistic, (11)

with model M1 describing the general case of an image plus a background model, and M0 describing only the instrumental background model. The likelihood ℒ(D|M0) of the data D given the background model then describes the null-hypothesis, which is tested against the alternative ℒ(D|M1). 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 26Al 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 Myr−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.

thumbnail Fig. 9

Simulated full-sky γ-ray maps of the 1.8 MeV line from 26Al decay (left) and the 1.3 MeV line from 60Fe 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.

thumbnail 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).

thumbnail Fig. 11

Best fitting sky map with TS = 2061 out of 30 000 PSYCO models. It is based on GM00 with 700 pc scale height, SFR = 8 M yr−1, IMF K01, stellar models LC06, and explodability LC18.

6 Discussion

6.1 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 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 26Al 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, and discussion in Krause et al. 2013; Krause & Diehl 2014). Also, as discussed above, the wind-and SN-blown cavities may not be filled homogeneously with 26Al and 60Fe, 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.

thumbnail 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.

6.2 Star formation and supernova rate

The fluxes measured by SPI in the 1.8 MeV line from 26Al and the 1.3 MeV line from 60Fe, 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 26Al and 60Fe in the Galaxy to be 1.2–2.4 M and 1–6 M, respectively.

6.3 60Fe/26Al ratio

The galactic-wide mass ratio of 60Fe/26Al 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 F60/F26 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 26Al flux compared to the PSYCO values. Measurements of the 60Fe y-ray lines are currently difficult because radioactive buildup of 60Co in SPI leads to an ever-increasing background in these lines (Wang et al. 2020). Therefore, the information value of the current Galactic 60Fe/26Al should not be over-interpreted. In fact, Wang et al. (2020) also estimated systematic uncertainties of the 60Fe/26Al flux ratio to reach up to 0.4, which would be consistent with a range of parameter combinations of PSYCO. In fact, an increased 60Fe flux from measurements might also be possible considering current detection significances of the two 60Fe 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 26Al throughout the Galaxy compared to 60Fe might enhance the observed discrepancy even further.

The total 26Al flux is consistently underestimated in PSYCO compared to measurements. Either the total 26Al 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 26Al with 60Fe 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 26Al emission. Wang et al. (2020) found a 60Fe scale radius of kpc and a scale height of kpc. Using the measured 26Al values by Pleintinger et al. (2019); Wang et al. (2020) of kpc and kpc, respectively, it is clear that meaningful comparisons rely on the quality of 60Fe measurements. In our simulations, about 50% of the total 26Al line flux originates from within a radius between 2.8 and 6.0 kpc. The 60Fe 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 60Fe emission. If strong foreground sources are present, this fraction can even be as high as 90% for 26Al, however remains rather stable for 60Fe.

Our density model GM00 best describes the SPI 26Al 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.

6.4 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 26Al later. Thus, larger bubbles would appear brighter. The average diameter of 26Al-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 feedbackprocesses. 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. (2013, 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.

6.5 Superbubble blowout and Galactic wind

Using simulations similar to Rodgers-Lee et al. (2019), Krause et al. (2021) have characterised the vertical blowout of super-bubbles from the Galactic disc. In some of their simulations, the superbubbles tend to merge in the disc and create transsonic, 26Al 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 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).

7 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 26Al 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 26Al 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 26Al 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, 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 26Al emission will accumulate with the remaining INTEGRAL mission until 2029. The 60Fe γ-ray lines, however, have become difficult to analyse because radioactive buildup of 60Co 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 26Al and 60Fe 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.

Acknowledgements

Thomas Siegert acknowledges support by the Bundesministerium für Wirtschaft und Energie via the Deutsches Zentrum für Luft- und Raumfahrt (DLR) under contract number 50 OX 2201.

Appendix A Details on initial mass functions

Kroupa (2001) refined the IMF slope as multiple-broken power-law (A.1)

with a normalisation constant ki in each of the three mass regimes depending on the local SFR. This was adjusted by Chabrier (2003) to follow a smooth log-normal distribution (A.2)

with normalisation constant kC05 = 1/log(10)exp [−(log(µ)2/2σ2)]. The log-normal parameters are amplitude a = 0.086, µ = 0.2, and σ = 0.55 (Chabrier 2005).

Appendix B Optimal sampling

To achieve optimal sampling with a total mass MEC, the IMF has to be normalised by excluding unphysical objects above Mmax, which is determined by the relation (Weidner & Kroupa 2004) (B.1)

with a truncation condition Mtrunc = ∞. For a complete sampling of the IMF, it can then be discretised by the condition (B.2)

which ensures that within each mass interval mi > mi+1 exactly one object is formed. This yields the iterative formula for mass segments (B.3)

with the normalisation k of the IMF. The individual masses of each star in mass segment i is then given by (B.4)

This sets the discretisation of a total cluster mass into single star masses in our population synthesis model.

Appendix C Line of sight integration for spherical superbubbles

For an observer at position p0 = (0,0,0) with respect to the centre of the emitting sphere s0 = (x0, y0, z0), γ rays can be received along each line of sight (C.1)

for azimuth angle ϕ in the Galactic plane and zenith angle θ. Only along lines of sight intersecting the sphere photons are received. These intersections are calculated as (C.2)

for the nearby and distant points from the observer smin and smax, respectively. Here, the auxiliary variables p0 := x0 cos ϕ cos θ + y0 sin ϕ cos θ + z0 sin θ and were introduced. For an observed outside a sphere, the differential flux across the surface of the sky is therefore (C.3)

where is the emissivity (in units of s−1 cm−3) with the variables defined in Eq. (9), and V being the volume contributing to the total luminosity as (C.4)

For an observer inside a sphere, the lower integration limit smin equals zero, and equations Eq. (C.3) and (C.4) change accordingly. We note that the volume integral is equal to the volume of a sphere of .

Appendix D Trends of SFR, SN rate, and 60Fe/26Ai as a function of explodability

In Figs. D.1 and D.2, we show the effect of different explod-ability assumptions in the literature on the SFR, SN rate, and 60Fe/26Al mass ratio. It is clear that for the rigorous assumption of no SNe above an initial stellar mass of 25 M (LC18; Limongi & Chieffi 2018), both the SN rate as well as the 60Fe/26Al mass ratio are decreased compared to other assumptions that include more SNe also for higher initial masses. Especially for the measured 60Fe/26Al flux ratio of ~ 0.2 with systematic deviations up to ~ 0.4 (Wang et al. 2020), mass ratios of 1.7 up to 3.4 would be suggested. However, the flux in the Galaxy is not concentrated in one point and our PSYCO simulations suggest that 26Al and 60Fe are not co-spatial, so that the rigorous assumption by Limongi & Chieffi (2018) may still be valid.

thumbnail Fig. D.1

Supernova rate as a function of star formation rate for different explodability assumptions.

thumbnail Fig. D.2

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

References

  1. Alexis, A., Jean, P., Martin, P., & Ferriere, K. 2014, A&A, 564, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Bally, J. 2008, Handbook of Star Forming Regions, I, 459 [Google Scholar]
  3. Banerjee, P., Heger, A., & Qian, Y.-Z. 2019, ApJ, 887, 187 [Google Scholar]
  4. Bazan, G., Brown, L. E., Clayton, D. D., et al. 1993, AIP Conf. Proc., 280, 47 [NASA ADS] [CrossRef] [Google Scholar]
  5. Belokurov, V., Penoyre, Z., Oh, S., et al. 2020, MNRAS, 496, 1922 [Google Scholar]
  6. Bennett, M. B., Wrede, C., Chipps, K. A., et al. 2013, Phys. Rev. Lett., 111, 232503 [Google Scholar]
  7. Bouchet, L., Jourdain, E., & Roques, J.-P. 2015, ApJ, 801, 142 [NASA ADS] [CrossRef] [Google Scholar]
  8. Breitschwerdt, D., Feige, J., Schulreich, M. M., et al. 2016, Nature, 532, 73 [NASA ADS] [CrossRef] [Google Scholar]
  9. Brinkman, H. E., Doherty, C. L., Pols, O. R., et al. 2019, ApJ, 884, 38 [NASA ADS] [CrossRef] [Google Scholar]
  10. Brown, A. G. A., Hartmann, D., & Burton, W. B. 1995, A&A, 300, 903 [NASA ADS] [Google Scholar]
  11. Burgess, J., & Capel, F. 2021, J. Open Source Softw., 6, 3257 [NASA ADS] [CrossRef] [Google Scholar]
  12. Burrows, D. N., Singh, K. P., Nousek, J. A., Garmire, G. P., & Good, J. 1993, ApJ, 406, 97 [NASA ADS] [CrossRef] [Google Scholar]
  13. Cappellari, M., & Copin, Y. 2003, MNRAS, 342, 345 [Google Scholar]
  14. Castor, J., McCray, R., & Weaver, R. 1975, ApJ, 200, L107 [NASA ADS] [CrossRef] [Google Scholar]
  15. Cerviño, M. 2013, New Astron. Rev., 57, 123 [CrossRef] [Google Scholar]
  16. Cerviño, M., & Luridiana, V. 2006, A&A, 451, 475 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
  18. Chabrier, G. 2005, The Initial Mass Function 50 Years Later, 327, eds. E. Corbelli, & F. Palle, 41 [Google Scholar]
  19. Cheng, J. Y., Rockosi, C. M., Morrison, H. L., et al. 2012, ApJ, 746, 149 [NASA ADS] [CrossRef] [Google Scholar]
  20. Chieffi, A., & Limongi, M. 2013, ApJ, 764, 21 [Google Scholar]
  21. Chomiuk, L., & Povich, M. S. 2011, AJ, 142, 197 [Google Scholar]
  22. Choplin, A., & Hirschi, R. 2020, J. Phys.: Conf. Ser., 1668, 012006 [CrossRef] [Google Scholar]
  23. Cordes, J. M., & Lazio, T. J. W. 2002, ArXiv e-prints [arXiv:astro-ph/0207156] [Google Scholar]
  24. Crocker, R. M. 2012, MNRAS, 423, 3512 [Google Scholar]
  25. da Silva, R. L., Fumagalli, M., & Krumholz, M. R. 2012, ApJ, 745, 145 [CrossRef] [Google Scholar]
  26. de Avillez, M. A., & Breitschwerdt, D. 2005, A&A, 436, 585 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Diehl, R., Dupraz, C., Bennett, K., et al. 1995, A&A, 298, 445 [NASA ADS] [Google Scholar]
  28. Diehl, R., Oberlack, U., Knödlseder, J., et al. 1997, in Proceedings of the fourth compton symposium (AIP), 1114 [CrossRef] [Google Scholar]
  29. Diehl, R., Knödlseder, J., Lichti, G. G., et al. 2003, A&A, 411, L451 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Diehl, R., Kretschmer, K., Lichti, G., et al. 2004, in Proceedings of the th INTEGRAL Workshop on the INTEGRAL Universe, 27 [Google Scholar]
  31. Diehl, R., Halloin, H., Kretschmer, K., et al. 2006, Nature, 439, 45 [CrossRef] [PubMed] [Google Scholar]
  32. Diehl, R., Lang, M. G., Martin, P., et al. 2010, A&A, 522, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Diehl, R., Siegert, T., Greiner, J., et al. 2018, A&A, 611, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Drimmel, R. 2002, New Astron. Rev., 46, 585 [CrossRef] [Google Scholar]
  35. Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146 [Google Scholar]
  36. Endal, A. S., & Sofia, S. 1978, ApJ, 220, 279 [NASA ADS] [CrossRef] [Google Scholar]
  37. Faucher-Giguere, C.-A., & Kaspi, V. M. 2006, ApJ, 643, 332 [Google Scholar]
  38. Fierlinger, K. M. 2014, PhD Thesis, LMU Munich, Germany [Google Scholar]
  39. Fierlinger, K. M., Burkert, A., Diehl, R., et al. 2012, Adv. Comput. Astrophys.: Methods, 453, 25 [NASA ADS] [Google Scholar]
  40. Fujimoto, Y., Krumholz, M. R., & Tachibana, S. 2018, MNRAS, 480, 4025 [NASA ADS] [CrossRef] [Google Scholar]
  41. Fujimoto, Y., Krumholz, M. R., & Inutsuka, S.-I.. 2020, MNRAS, 497, 2442 [NASA ADS] [CrossRef] [Google Scholar]
  42. Glebocki, R., & Gnacinski, P. 2005, VizieR On-line Data Catalog, III/244 [Google Scholar]
  43. Glebocki, R., Gnacinski, P., & Stawikowski, A. 2000, Acta Astron., 50, 509 [Google Scholar]
  44. Gong, M., & Ostriker, E. C. 2015, ApJ, 806, 31 [NASA ADS] [CrossRef] [Google Scholar]
  45. Hartmann, D. H. 1994, AIP Conf. Proc., 304, 176 [NASA ADS] [CrossRef] [Google Scholar]
  46. Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368 [NASA ADS] [CrossRef] [Google Scholar]
  47. Heitsch, F., Hartmann, L. W., Slyz, A. D., Devriendt, J. E. G., & Burkert, A. 2008, ApJ, 674, 316 [NASA ADS] [CrossRef] [Google Scholar]
  48. Janka, H.-T. 2012, Annu. Rev. Nuclear Particle Sci., 62, 407 [NASA ADS] [CrossRef] [Google Scholar]
  49. Janka, H.-T., Hanke, F., Hüdepohl, L., et al. 2012, Progr. Theor. Exp. Phys., 2012, 01A309 [Google Scholar]
  50. Jones, S. W., Möller, H., Fryer, C. L., et al. 2019, MNRAS, 485, 4287 [NASA ADS] [CrossRef] [Google Scholar]
  51. Kavanagh, P. J. 2020, Astrophys. Space Sci., 365, 6 [NASA ADS] [CrossRef] [Google Scholar]
  52. Keller, B. W., Wadsley, J., Benincasa, S. M., & Couchman, H. M. P. 2014, MNRAS, 442, 3013 [NASA ADS] [CrossRef] [Google Scholar]
  53. Keller, B. W., Wadsley, J., & Couchman, H. M. P. 2016, MNRAS, 463, 1431 [Google Scholar]
  54. Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531 [NASA ADS] [CrossRef] [Google Scholar]
  55. Knödlseder, J. 1999, ApJ, 510, 915 [CrossRef] [Google Scholar]
  56. Knödlseder, J. 2000, A&A, 360, 539 [NASA ADS] [Google Scholar]
  57. Knoedlseder, J., Prantzos, N., Bennett, K., et al. 1996, A&AS, 120, 335 [NASA ADS] [Google Scholar]
  58. Knödlseder, J., Bennett, K., Bloemen, H., et al. 1999, A&A, 344, 68 [NASA ADS] [Google Scholar]
  59. Krause, M., & Diehl, R. 2014, ApJ, 794, L21 [NASA ADS] [CrossRef] [Google Scholar]
  60. Krause, M., Fierlinger, K., Diehl, R., et al. 2013, A&A, 550, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Krause, M., Diehl, R., Böhringer, H., Freyberg, M., & Lubos, D. 2014, A&A, 566, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Krause, M., Diehl, R., Bagetakos, Y., et al. 2015, A&A, 578, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Krause, M., Burkert, A., Diehl, R., et al. 2018, A&A, 619, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Krause, M. G. H., Rodgers-Lee, D., Dale, J. E., Diehl, R., & Kobayashi, C. 2021, MNRAS, 501, 210 [Google Scholar]
  65. Kretschmer, K., Plüschke, S., Diehl, R., & Hartmann, D. H. 2000, Astron. Gesellsch. Meet. Abstr., 17, 77 [Google Scholar]
  66. Kretschmer, K., Diehl, R., Krause, M., et al. 2013, A&A, 559, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  68. Kroupa, P., & Jerabkova, T. 2021, ArXiv e-prints [arXiv:2112.10788] [Google Scholar]
  69. Kroupa, P., Weidner, C., Pflamm-Altenburg, J., et al. 2013, in Planets, Stars and Stellar Systems (Dordrecht: Springer, Dordrecht), 115 [CrossRef] [Google Scholar]
  70. Krumholz, M. R., McKee, C. F., & Bland-Hawthorn, J. 2019, ARA&A, 57, 227 [NASA ADS] [CrossRef] [Google Scholar]
  71. Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57 [Google Scholar]
  72. Lentz, E. J., Branch, D., & Baron, E. 1999, ApJ, 512, 678 [NASA ADS] [CrossRef] [Google Scholar]
  73. Licquia, T. C., & Newman, J. A. 2015, ApJ, 806, 96 [NASA ADS] [CrossRef] [Google Scholar]
  74. Limongi, M., & Chieffi, A. 2006, ApJ, 647, 483 [NASA ADS] [CrossRef] [Google Scholar]
  75. Limongi, M., & Chieffi, A. 2018, ApJS, 237, 13 [NASA ADS] [CrossRef] [Google Scholar]
  76. Luhman, K. L. 2012, ARA&A, 50, 65 [CrossRef] [Google Scholar]
  77. Maíz-Apellániz, J. 2008, ApJ, 677, 1278 [CrossRef] [Google Scholar]
  78. Martin, P., Knödlseder, J., & Meynet, G. 2008, New Astron. Rev., 52, 445 [CrossRef] [Google Scholar]
  79. Martin, P., Knoedlseder, J., Diehl, R., & Meynet, G. 2009, A&A, 506, 703 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Martin, P., Knödlseder, J., Meynet, G., & Diehl, R. 2010, A&A, 511, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Meynet, G., Arnould, M., Prantzos, N., & Paulus, G. 1997, A&A, 320, 460 [NASA ADS] [Google Scholar]
  82. Micic, M., Glover, S. C. O., Banerjee, R., & Klessen, R. S. 2013, MNRAS, 432, 626 [NASA ADS] [CrossRef] [Google Scholar]
  83. Misiriotis, A., Xilouris, E. M., Papamastorakis, J., Boumis, P., & Goudis, C. D. 2006, A&A, 459, 113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Mowlavi, N., & Meynet, G. 2006, New Astron. Rev., 50, 484 [CrossRef] [Google Scholar]
  85. Murray, N., & Rahman, M. 2010, ApJ, 709, 424 [NASA ADS] [CrossRef] [Google Scholar]
  86. Nomoto, K., Kobayashi, C., & Tominaga, N. 2013, ARA&A, 51, 457 [CrossRef] [Google Scholar]
  87. Norris, T. L., Gancarz, A. J., Rokop, D. J., & Thomas, K. W. 1983, Am. Geophys. Union NASA, 14, B331 [Google Scholar]
  88. Oberlack, U., Bennett, K., Bloemen, H., et al. 1995, Int. Cosmic Ray Conf., 2, 207 [NASA ADS] [Google Scholar]
  89. Oberlack, U., Bennett, K., Bloemen, H., et al. 1996, A&A Suppl. Ser., 120, 311 [NASA ADS] [Google Scholar]
  90. Oey, M. S., & Clarke, C. J. 2005, ApJ, 620, L43 [NASA ADS] [CrossRef] [Google Scholar]
  91. Pflamm-Altenburg, J., Weidner, C., & Kroupa, P. 2007, ApJ, 671, 1550 [NASA ADS] [CrossRef] [Google Scholar]
  92. Pflamm-Altenburg, J., González-Lópezlira, R. A., & Kroupa, P. 2013, MNRAS, 435, 2604 [CrossRef] [Google Scholar]
  93. Piskunov, A. E., Kharchenko, N. V., Schilbach, E., et al. 2011, A&A, 525, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Pleintinger, M. M. M. 2020, PhD Dissertation, Technical University of Munich, Germany [Google Scholar]
  95. Pleintinger, M. M. M., Siegert, T., Diehl, R., et al. 2019, A&A, 632, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Plüschke, S. 2001, PhD thesis, Technische Universität München, Germany [Google Scholar]
  97. Plüschke, S., Diehl, R., Schönfelder, V., et al. 2001, in Proceedings of the Fourth INTEGRAL Workshop, eds. B. Battrick, A. Giménez, V. Reglero, & C. Winkler, ESA Publications Division, ESA SP-459, Alicante, Spain, 55–58 [Google Scholar]
  98. Popescu, B., & Hanson, M. M. 2009, AJ, 138, 1724 [NASA ADS] [CrossRef] [Google Scholar]
  99. Prantzos, N. 1993, AIP Conf. Proc., 280, 52 [NASA ADS] [CrossRef] [Google Scholar]
  100. Prantzos, N., & Diehl, R. 1995, Adv. Space Res., 15, 99 [NASA ADS] [CrossRef] [Google Scholar]
  101. Prantzos, N., & Diehl, R. 1996, Phys. Rep., 267, 1 [NASA ADS] [CrossRef] [Google Scholar]
  102. Prantzos, N., Abia, C., Cristallo, S., Limongi, M., & Chieffi, A. 2020, MNRAS, 491, 1832 [Google Scholar]
  103. Predehl, P., Sunyaev, R. A., Becker, W., et al. 2020, Nature, 588, 227 [CrossRef] [Google Scholar]
  104. Reid, M. J., Dame, T. M., Menten, K. M., & Brunthaler, A. 2016, ApJ, 823, 77 [Google Scholar]
  105. Renzo, M., Zapartas, E., de Mink, S. E., et al. 2019, A&A, 624, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Rodgers-Lee, D., Krause, M., Dale, J., & Diehl, R. 2019, MNRAS, 490, 1894 [NASA ADS] [CrossRef] [Google Scholar]
  107. Rugel, G., Faestermann, T., Knie, K., et al. 2009, PRL, 103, 072502 [CrossRef] [Google Scholar]
  108. Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
  109. Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [Google Scholar]
  110. Schulz, C., Pflamm-Altenburg, J., & Kroupa, P. 2015, A&A, 582, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  111. Siegert, T. 2019, A&A, 632, A1 [Google Scholar]
  112. Siegert, T., & Diehl, R. 2017, Proc. 14th International Symposium on Nuclei in the Cosmos (NIC2016), 020305 [Google Scholar]
  113. Siegert, T., Diehl, R., Weinberger, C., et al. 2019, A&A, 626, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  114. Smartt, S. J. 2009, ARA&A, 47, 63 [NASA ADS] [CrossRef] [Google Scholar]
  115. Sofue, Y. 2000, ApJ, 540, 224 [NASA ADS] [CrossRef] [Google Scholar]
  116. Sturner, S. J. 2001, in Proceedings of the Fourth Integral Workshop, 101 [Google Scholar]
  117. Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H.-T. 2016, ApJ, 821, 38 [Google Scholar]
  118. Timmes, F. X., Woosley, S. E., Hartmann, D. H., et al. 1995, ApJ, 449, 204 [NASA ADS] [CrossRef] [Google Scholar]
  119. Tinsley, B. M. 1968, ApJ, 151, 547 [Google Scholar]
  120. Tinsley, B. M. 1972, A&A, 20, 383 [NASA ADS] [Google Scholar]
  121. Tomsick, J. A., Zoglauer, A., Sleator, C., et al. 2019, Bull. Am. Astron. Soc., 51, 98 [Google Scholar]
  122. Tur, C., Heger, A., & Austin, S. M. 2010, ApJ, 718, 357 [NASA ADS] [CrossRef] [Google Scholar]
  123. Uberseder, E., Adachi, T., Aumann, T., et al. 2014, PRL, 112, 211101 [NASA ADS] [CrossRef] [Google Scholar]
  124. Vallée, J. P. 2008, AJ, 135, 1301 [CrossRef] [Google Scholar]
  125. Vedrenne, G., Roques, J.-P., Schönfelder, V., et al. 2003, A&A, 411, L63 [CrossRef] [EDP Sciences] [Google Scholar]
  126. Voss, R., Diehl, R., Hartmann, D. H., et al. 2009, A&A, 504, 531 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  127. Voss, R., Diehl, R., Vink, J. S., & Hartmann, D. H. 2010, A&A, 520, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  128. Voss, R., Martin, P., Diehl, R., et al. 2012, A&A, 539, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  129. Wainscoat, R. J., Cohen, M., Volk, K., Walker, H. J., & Schwartz, D. E. 1992, ApJS, 83, 111 [NASA ADS] [CrossRef] [Google Scholar]
  130. Wang, W., Harris, M. J., Diehl, R., et al. 2007, A&A, 469, 1005 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  131. Wang, W., Siegert, T., Dai, Z. G., et al. 2020, ApJ, 889, 169 [NASA ADS] [CrossRef] [Google Scholar]
  132. Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377 [Google Scholar]
  133. Weidner, C., & Kroupa, P. 2004, MNRAS, 348, 187 [NASA ADS] [CrossRef] [Google Scholar]
  134. Winkler, C., Courvoisier, T. J. L., Di Cocco, G., et al. 2003, A&A, 411, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  135. Woosley, S. E., & Heger, A. 2007, Phys. Rep., 442, 269 [Google Scholar]
  136. Yan, Z., Jerabkova, T., & Kroupa, P. 2017, A&A, 607, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  137. Yang, H. Y. K., Ruszkowski, M., & Zweibel, E. G. 2022, Nat. Astron., 6, 584 [NASA ADS] [CrossRef] [Google Scholar]
  138. Yusifov, I., & Küçük, I. 2004, A&A, 422, 545 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Tables

Table 1

Spiral centroid parameters (cf. Eq. (2)) following Faucher-Giguère & Kaspi (2006).

Table 2

Linear parameters for modelling the radial metallicity gradient in the Milky Way according to different heights above the Galactic plane (Cheng et al. 2012).

Table 3

Conversions between 60Fe/26Al flux, mass, isotopic, and production ratios.

Table 4

Fluxes of 26Al emission at 1.809 MeV and of 60Fe 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.

All Figures

thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail Fig. 3

26AI (blue) and 60Fe (green) yields from SNe by Limongi & 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.

In the text
thumbnail Fig. 4

Nucleosynthesis yields (top 26Al, bottom 60Fe) for a selection of stellar evolution models, separated into the stellar wind and total contribution (wind + SN).

In the text
thumbnail Fig. 5

Population synthesis a canonical star group of 104 M including binary effects. Shown is the impact of 26Al binary wind yields from Brinkman 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.

In the text
thumbnail Fig. 6

Population synthesis models of nucleosynthesis ejecta in a 104 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 26Al (upper), 60Fe (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 103 Monte Carlo runs.

In the text
thumbnail Fig. 7

Steady-state settling time of the total 60Fe (blue) and 26Al (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.

In the text
thumbnail Fig. 8

Radial distribution of modelled flux contributions for a theoretical observer (white circle) from 26Al (top) and 60Fe (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 26Al and 60Fe, 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.

In the text
thumbnail Fig. 9

Simulated full-sky γ-ray maps of the 1.8 MeV line from 26Al decay (left) and the 1.3 MeV line from 60Fe 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.

In the text
thumbnail 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).

In the text
thumbnail Fig. 11

Best fitting sky map with TS = 2061 out of 30 000 PSYCO models. It is based on GM00 with 700 pc scale height, SFR = 8 M yr−1, IMF K01, stellar models LC06, and explodability LC18.

In the text
thumbnail 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.

In the text
thumbnail Fig. D.1

Supernova rate as a function of star formation rate for different explodability assumptions.

In the text
thumbnail Fig. D.2

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

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.