Chemical evolution models: the role of type Ia supernovae in the α-elements over iron relative abundances and their variations in time and space

ABSTRACT

The role of type Ia supernovae (SN Ia), mainly the delay time distributions (DTDs) determined by the binary systems, and the yields of elements created by different explosion mechanisms, are studied by using the MulChem chemical evolution model applied to our Galaxy. We explored 15 DTDs and 12 tables of elemental yields produced by different SN Ia explosion mechanisms, doing a total of 180 models. Chemical abundances for |$\alpha$|-elements (O, Mg, Si, and Ca) and Fe derived from these models are compared with recent solar region observational data of |$\alpha$|-elements over Fe relative abundances, [X/Fe], as a function of [Fe/H] and age. A multidimensional maximum-likelihood analysis shows that 52 models are able to fit all these data sets simultaneously, considering the 1|$\sigma$| confidence level. The combination of STROLG1 DTD from Strolger et al. (2020) and LN20181 SN Ia yields from Leung & Nomoto (2018) provides the best fit. The exponential model with very prompt events is a possible DTD, but a combination of several channels is more probable. The SN Ia yields that include MCh or Near MCh correspond to 39 (75 per cent) of the 52 best models. Regarding the DTD, 31 (60 per cent) of the 52 most probable models correspond to the SD scenario, while the remaining 21 (40 per cent) are based on the DD scenario. Our results also show that the relatively large dispersion of the observational data may be explained by the stellar migration from other radial regions, and/or perhaps a combination of DTDs and explosion channels.

1 INTRODUCTION

Supernovae (SNe) explosions are one of the key ingredients in chemical evolution models, since they produce most of elements heavier than N. However, each SN type produces different proportions of elements. As an example, a typical type Ia supernova (SN Ia) produces around 0.01 M|$_{\odot }$| of Ca (Seitenzahl et al. 2013a), while a core-collapse (CC) SN progenitor produces between 0.05 to 0.10 M|$_{\odot }$| (Nomoto, Kobayashi & Tominaga 2013). Massive stars, through CC SNe, pollute the interstellar medium (ISM) both with |$\alpha$|-elements, as, for example, O, Mg, Si, S, and Ca (Weinberg, Andrews & Freudenburg 2017), produced during the stellar evolution, and also iron peak elements, such as Cr, Mn, Fe, Co, and Ni, created in the CC SNe explosion. Instead, type Ia SNe mainly pollutes the ISM with the latter, leading to a decrease of the [|$\alpha$|/Fe] with increasing metallicity when they appear.

The importance of SNe Ia in the chemical evolution of galaxies resides in this amount of Fe released in these events. The mean Fe yield for a SN Ia is |$y_{\rm Fe,Ia} \sim 0.7$| M|$_{\odot }$| (Iwamoto et al. 1999; Mazzali et al. 2007; Howell et al. 2009). Recent works give values in the range |$[0.622-0.790]$| M|$_{\odot }$| per event (see Section 3.3), while for a typical CC SNe, it is a factor of 10 lower. When considering the percentage of explosions of all types of CC SNe (II, Ib, Ic, etc.) (Li et al. 2011b; Shivvers et al. 2017; Graur et al. 2017a, b), an average yield |$\lt y_{\rm Fe,CC}\gt = 0.066$| M|$_{\odot }$| is obtained.

Besides the yield of each individual SNe, it is necessary to consider their rates of production: there are between 7 and 8 times more CC SNe than SNe Ia, given their Hubble-time integrated production efficiency |$N/M_{*}$|⁠, that is, the expected SN number by solar mass formed. For CC SNe, this number is |$N_{\mathrm{ CC}}/M_{*} \sim 0.01$| M|$_{\odot }$||$^{-1}$|⁠, depending on the assumed initial mass function (IMF), while for SNe Ia it is in the range |$N_{\mathrm{ Ia}}/M_{*} \sim [0.0016-0.054]$| M|$_{\odot }$||$^{-1}$|⁠, with field galaxies displaying the lowest value (Maoz & Graur, 2017). This way, SNe Ia and CC SNe contribute approximately to half of the Fe that is formed in a typical galaxy.

The stellar abundances may reflect this production of Fe through the relative abundance [|$\alpha$|/Fe], which retains the CC SN and massive stars production plus the contribution of SNe Ia. In other words, for any |$\alpha$|-element as O, Mg, Si, S, Ca, or Ne, the relative abundances [|$\alpha$|/Fe] depend on the number of SNe Ia that have exploded compared with the number of CC SNe. The ratio [|$\alpha$|/Fe] is, therefore, an indicator of the evolution of the region, mostly how the different SN rates (SNRs) have varied over time. If CC SN events were the only contributor to the Fe yield, there would exist a plateau in approximately [|$\alpha$|/Fe] |$\sim 0.3$|– 0.5 (Maoz & Graur 2017), defined by the corresponding yields of |$\alpha$|-elements and Fe from massive stars and CC SNe. In fact, this is the observed value for the metal-poorest stars of the Milky Way Galaxy (MWG).

In turn, SN Ia progenitors are expected to reside in binary stellar systems, involving at least one white dwarf (WD) that accretes mass from a companion star increasing its mass, which approaches to the Chandrasekhar mass (MCh). The WD will explode after the evolution of the most massive component; that is, there is a delay time after the creation of the binary system before the Fe created in that SNe Ia is ejected into the ISM. This delay depends on the scenario involving the evolution of the binary system to reach the explosion. The final consequence is a decreasing in [|$\alpha$|/Fe] from the initial ratio 0.3–0.5 from massive stars, once the Fe from SN Ia appears in the ISM. Thus, the chemical evolution may give clues about the proportions of both types of contributors.

The time distributions of SNRs are evaluated by using the so-called delay time distribution (DTD) function, which is the expected SNR after an instantaneous star formation burst or simple stellar population (SSP). A minimum time of delay is necessary to allow the most massive star to evolve, which would be a few tens of Myr for a star of 8 M|$_{\odot }$|⁠. Then, the DTD must be modulated by the star formation rate (SFR) of each region or galaxy in order to compute the total SNR within it. In recent years, improved estimates for DTDs of different SN species (CC and Ia) have been available. There are different methods to empirically infer the SN Ia DTD (see Maoz, Mannucci & Nelemans 2014, for a review) from the existing data. From the theoretical side, the SN Ia DTD was computed some years ago using binary star population synthesis (bps) models (Yungelson & Livio 2000; Greggio 2005). These simulations computed the DTD of the different scenarios by using probability distribution functions for stellar orbital parameters and following the mass transfer in the system.

There are important differences in the evolution after the birth of the binary system until the explosion, depending on the assumed scenarios for that system. The two common ones are the single degenerate (SD) systems, where a single WD accretes mass from a non-WD star, such as a main-sequence star, a red super giant, or a He-rich star; and the double degenerate (DD) systems, where both stars are WD, with at least one CO WD.

For the DD case, the time to reach the merger strongly depends on the separation of the WDs in the post-common envelope phase (Totani et al. 2008). The bps models give for this case a DTD |$\propto t^{-1}$| (Katz & Dong 2012; Toonen, Nelemans & Portegies Zwart 2012) after a certain delay. There exist models showing a great variety of results depending on their parameters, as, for example, the exponential model (Strolger et al. 2020) but, in general, they predict a time dependence compatible with a power law, even with short delays (Hachisu, Kato & Nomoto 2008; Meng & Yang 2010; Mennekens et al. 2010; Bours, Toonen & Nelemans 2013). For times shorter than 30–40 Myr, which is the mean lifetime of stars with 8 M|$_{\odot }$|⁠, we do not expect any SN Ia event, since this is the necessary time for the binary system, assuming it has a total mass in the range [3—16] M|$_{\odot }$|⁠, evolves.

For the SD scenario, the mass of the non-degenerate donor star was very restricted in the early models, only 2–3 M|$_{\odot }$|⁠, by the accretion conditions (Langer et al. 2000). Recent studies, however, include different interactions in their bps codes such as Roche lobe overflows or wind ejecta interactions, expanding the mass range up to 8 M|$_{\odot }$| stars and some red giant stars. Depending on this mass, the resulting delay may be as short as |$\sim$| 100 Myr or even 40 Myr (Greggio 2005).

Other studies use different approaches to recover the DTD. Most of these observational DTDs find a similar dependence |$\propto t^{-1}$| in a continuous way with no cut-off times, which supports the DD scenario. Sometimes a prompt population of the order of the 50  per cent of SNe Ia taking place with delay times shorter than 500 Myr is found, which seems to indicate the SD scenario as probable. There exist, however, some tension about the existence or not of such short delays, without a consensus about the dominant channel for the SNe Ia.

On the other hand, the nucleosynthetic yields produced by SNe Ia are dependent on the explosion mechanisms, and finally determined by the density of the matter at which the thermonuclear burning starts. According to this, the yields of elements ejected by SNe Ia are usually divided in two categories depending on the mass of the WD: some SNe Ia explode when the Chandrasekhar mass is reached – or almost reached – (MCh explosions), but there also exists the possibility that the WD explodes before to reach that MCh, giving place to a sub-MChandrasekhar (sub-MCh) explosion.¹

The relationship between the SN Ia formation scenarios (SD or DD) and the WD masses involved in the explosions (sub-MCh and MCh) is a complex question and a current open debate in astronomy (e.g Ruiter 2020; Pakmor et al. 2022; Liu, Röpke & Han 2023. There is no clear consensus either about the SN Ia progenitor system and the actual explosion mechanisms, that is, the two unknown points: which is the companion star donating the mass to the WD, (and the process of mass transfer and its time-scale) and which is the WD mass or density when the SN Ia explodes are still unclear. Therefore, this dichotomy SD/DD or sub-MCh/MCh has changed very much in the last years. The idea that SN Ia originating in a given scenario tend to be associated with WD reaching a certain mass, is beginning to be considered too simple. In practice, it is not always easy to determine the exact mass of a WD before its explosion as an SN Ia. Actually, the diversity of the observed SN Ia is large, as it is shown in fig. 1 from Taubenberger (2017), what makes it difficult to classify them in only two categories.

In last years, there have been a lot of works developed around the explosion mechanisms besides the classic near-MCh explosion of a WD in a SD binary system (Gilfanov & Bogdán 2010; Sim et al. 2013), or the violent merger of two WDs in a DD systems (Iben & Tutukov 1984; Guillochon et al. 2010; Pakmor et al. 2013), the two old basic channels. The mass accretion process and conditions in binary systems can vary widely. In some SD systems, the WD can accumulate enough mass to reach MCh, while in other cases, it can explode before reaching that critical mass (see, for instance, the SD models for both sub-MCh and MCh from Greggio 2005). In the DD scenario, although the WDs involved are expected to be less massive than MCh, not all mergers will result in sub-MCh explosions, as the details of the merger are complex and may vary. There are some studies that support that a large fraction of SNe Ia comes from this sub-MCh channel, either by SD or DD He mass transfer scenario or through DD violent mergers (Ruiter, Belczynski & Fryer 2009; Gilfanov & Bogdán 2010; Shappee, Kochanek & Stanek 2013; Goldstein & Kasen 2018; Kuuttila et al. 2019; Flörs et al. 2020). Calculations of DTD with a DD scenario for MCh are also available (Ruiter et al. 2011).

Soker (2019, see his table 1) claims that now there are until five or more different explosions channels, such as: (1) single degenerate with a time of delay until the explosion from the merger (SD with MED); (2) core degenerate (CD); (3) double degenerate (DD with MED); (4) double detonation (DDet); and (5) WD–WD collision. The CD and DD-MED models explode with MCh, while DD, SD-MED, and DDet do that with sub-MCh. In turn, Ruiter (2020, see table 1) established that an SN Ia may take place by: (i) an MCh WD that accretes mass in an SD scenario by a delayed detonation (or failed detonation) also named deflagration-to-denotation transition – DDT – model; (ii) a sub-MCh WD accreting helium-rich material that explodes by a double detonation in an SD scenario (DDet), through the He mass accretion from a companion (Nomoto 1982; Greggio 2005; Bildsten et al. 2007; Kromer et al. 2010; Ruiter et al. 2011; Woosley & Kasen 2011); or (iii) a double WD merger with MCh or sub-MCh through a delayed detonation or DDet. Nomoto & Leung (2019) clarified that ‘The sub-MCh mass explosions could occur in both SD and DD detonations’.² In their recent review, Liu et al. (2023, see their fig. 3) give a detailed explanation of these mechanisms about the different channels, where the authors specify that it is possible to have sub-MCh explosions in both SD and DD scenarios, and there are new scenarios with either one or two degenerate WD that allow explosions with MCh or sub-MCh. Several of these research groups have calculated different explosions with sophisticated burning and explosion simulations, obtaining the corresponding SN Ia yields (Iwamoto et al. 1999; Leung & Nomoto 2018, 2020; Shen et al. 2018; Eitner et al. 2020; Gronow et al. 2020, see details in Section 3.3).

Thus, our idea is to use the fact that the [|$\alpha$|/Fe] versus age and/or [|$\alpha$|/Fe] versus [Fe/H] relations are able to trace the chemical enrichment produced by SNe Ia along the evolution, seeking to distinguish among the described scenarios, in order to clarify which DTD is the most probable. Furthermore, we will vary the sets of yields for SN Ia coming from different authors by assuming different explosion mechanisms. Therefore, the basic objective of this work is to compute the same chemical evolution model for the MWG varying the DTD prescriptions as well as the SNe Ia yields, in order to, after comparing results with observations of the solar neighbourhood, to search for the best model, in particular which DTD and explosion mechanism are the most probable. We assume in this work that both inputs are independent, since most of works devoted to DTD or SN Ia yields give this information separately, by computing all possible combinations DTD|$+$|yields; however, we are conscious that some models/combinations could be less realistic than others. We will comment this problem when analysing the results.

The relative abundance [|$\alpha$|/Fe], sensitive to the SNRs and this way to the DTD convolved with the SFR, has been previously used in other chemical evolution models (Tinsley 1979; Matteucci & Greggio 1986; De Donder & Vanbeveren 2004; Matteucci et al. 2006; Kobayashi & Nomoto 2009; Tsujimoto & Shigeyama 2012) in MWG or in other nearby galaxies (Walcher et al. 2016). As an example, Calura, Matteucci & Tozzi (2007) find some evidence of SN Ia with short delays from the abundance–age relation. More recently, Palla (2021) has contributed to this subject by computing chemical evolution models with a set of SN Ia yields from different authors and simulating the three more probable explosion mechanisms within a model with only one DTD. They found that the classical W7 + WDD models (two channels) produce similar patterns that near-MCh and sub-MCh mass models.

Our compilation of observational data is presented in Section 2. Section 3 describes our MWG chemical evolution model, in particular the different DTD prescriptions under study, discussed in Section 3.2, and the yields used for SNe Ia described in Section 3.3. The results of the comparison model versus data are discussed in Section 4. Our main conclusions are summarized in Section 5. Furthermore, we give as Supporting Information the normalization process of the observational data sets in Appendix  A; and the calibration of the MWG model in Appendix  B.

2 OBSERVATIONAL DATA: STELLAR CATALOGUES

Stellar data will be used to compare with the [X/Fe] relative abundance from models and to be able to determine which of them is the best one to fit these data, the most ones affected by variations in DTDs and SN Ia yields. Here, we describe these data sets.

To compare the model predictions for the time evolution of the disc, we searched for observations of stars with available elemental abundances for Fe and for |$\alpha$|-elements (as many of them as possible) and also age estimates. To conduct the statistical analysis, we only select data sets that include all the analysed chemical abundances and ages. This allows us to construct a complete sample, meaning that all the plots analysed in this work contain the same stars. We used here the stellar data from several authors/surveys, as given in Table 1, in particular those from Chen et al. (2000), the HARPS-Guaranteed Time Observations (GTO) and the GALAH (Galactic Archaeology with HERMES)-DR3 samples, which fulfill these conditions. This table describes the references for each data set and the names adopted through this work to refer to the data sets.

Chen et al. (2000) gathered high-resolution and high signal-to-noise (S/N) spectra for a sample of 90 F and G main-sequence disc stars. The observations were performed with the coudé Echelle spectrograph attached to the 2.16 m telescope at Beijing Astronomical Observatory (Xinglong, PR China). They provide metallicities in the range |$-1.0 \lt $| [Fe/H] |$\lt +0.1$|⁠, with approximately the same number of stars in each metallicity bin of 0.1 dex. Stellar ages are also provided from interpolated evolutionary tracks. Their |$\alpha$|-chemical abundances for O, Mg, Si, and Ca were used in our analysis.

The HARPS (HARPS-GTO, Mayor et al. 2003; Lo Curto et al. 2010; Santos et al. 2011) sample is composed by the stellar data given in Delgado Mena et al. (2010, 2015, 2017, 2019), Bertran de Lis et al. (2015), and Suárez-Andrés et al. (2016, 2017), which are compilations from the HARPS planet search programme, HARPS-1, observed with the Ultraviolet and Visual Echelle Spectrograph (UVES) installed at the VLT/UT2 Kueyen telescope (Paranal Observatory, ESO, Chile) during several campaigns. We have used the catalogues³ taken from the CDS-Strasbourg. Bertran de Lis et al. (2015) provide stellar parameters and abundances for 762 stars observed at high resolution and high S/N. We prefer to use their oxygen abundances derived from O i 6158 Å line, since these authors suggest that this line can be used as a reliable indicator of oxygen abundance in solar-type stars. We selected [Fe/H] abundances and stellar ages for 65 disc dwarf stars from Suárez-Andrés et al. (2016). In Suárez-Andrés et al. (2017), the [C/H] and [Fe/H] abundances and stellar ages for 1058 FGK solar-type disc stars were selected. Delgado Mena et al. (2017, 2019) give abundances for Mg, Si, and Ca. The HARPS data from these different authors are, in principle, compatible, since the abundances are calculated with the same methods and codes. Using all sets together, we will simulate the real dispersion of data due to different methods and instruments. Finally, in Delgado Mena et al. (2019), the stellar ages are obtained, allowing the analysis of the time evolution in the solar vicinity. Additionally, we apply a quality selection criteria by choosing stars with precision better than 0.2 dex in the abundances and also those with the flag thin in the catalogues, to select only stars that pertain to the thin disc component.

We have also used the information from the second Gaia data release already included in the high-resolution GALAH spectroscopic survey (Buder et al., 2021). This survey aims at to analyse the structure of the Galactic disc components using a large number of stars in the MWG by studying the kinematics and chemical abundances. The total number of stars from this DR3 is |$\sim$| 588 571. However, not all of them have similar precision for stellar ages and abundances. We have selected stars with error in ages smaller than 15 per cent. As GALAH-DR3 also provides astrometric parameters from Gaia DR2 (Gaia Collaboration 2016, 2018), we have also selected stars with height over the Galactic plane |$|z| \lt 0.2$| kpc, and within the solar distance, with Galactocentric distances in the range |$7 \le R_{\rm gal} \le 9$| kpc, and we have restricted to good parallaxes with relative uncertainties lower than 10 per cent. The geometric estimated distances propagated from Bailer-Jones et al. (2018) were used and the heliocentric distances were converted to Galactocentric distances by adopting the Sun’s distance to the Galactic Centre |$R_{0} = 8.1$| kpc (GRAVITY Collaboration 2018). These cuts implies errors in the ages for the oldest stars of |$\sim 1.5$| Gyr, although the average value is slightly higher than 0.5 Gyr. We have also selected stars for which the abundances of the elements Fe, O, Mg, Si, and Ca are available, and also limiting their errors to 0.1 dex. This subsample, which consists of 5098 stars, is used through this work and referred as GALAH. More details about this subsample can be seen in Appendix  A.

All these data are given using different solar abundances, as shown in Appendix  A (Supporting Information, Table A1). Therefore, in order to use all of them together, we have normalized the data to the same solar abundances scale from Lodders (2019). After this, differences in their scale could still be present, since we consider that all abundances [Fe/H] and [X/Fe] must be zero for |${\rm [Fe/H]=0}$| and at the Sun age of 4.6 Gyr (when the time since the Galaxy formation is |$t=8.6$| Gyr). That is, the data must be around the solar values at the solar metallicity and age. We have therefore shifted each cloud of data to locate the solar value in its centre. Further details about the normalization process can be found in Appendix  A. In Table 1, we provide the sets of data used and the chemical abundances and ages used for the analysis. The number of stars containing in each data set is indicated in the second to last line and in the last line the references used for each data set. The resulting combined data, considering all the data sets together, are represented in Fig. 1 for the age–metallicity relation and in Figs 2 and 3 for the [X/Fe] relative abundances as a function of age and of [Fe/H], respectively.

3 SUMMARY DESCRIPTION OF GALAXY EVOLUTIONARY MODELS

Our MulChem chemical evolution models, particularly the one applied to the MWG, have been described in Mollá et al. (2015, 2016, 2017, 2019, hereinafter MOL15, MOL16, MOL17, and MOL19, respectively), updating the models from Ferrini et al. (1994) and Mollá & Díaz (2005) for what refers to the inputs for modelling a given galaxy. We summarize here the changes done in this new version of the Galaxy model, related with stellar yields and SN Ia DTDs. The other model hypotheses and equations are widely described in the cited references.

3.1 Stellar yields

To compute the elemental abundances, we use the technique based on the Q-matrix formalism (Talbot & Arnett 1973; Ferrini et al. 1992; Portinari, Chiosi & Bressan 1998). Each element |$(i,j)$| of a matrix, |$Q_{i,j}$| gives the proportion of a star which was initially element j, ejected as i when the star dies. Thus,

where |$m_{i,j,\mathrm{ exp}}$| is the ejected mass of an element i initially being other j. These numbers must be multiplied by the number of stars of mass m, given by the used IMF, here being the one from Kroupa (2001, hereinafter KRO) with limits |$m_{\rm low}=0.15$| M|$_{\odot }$| and |$m_{\rm up}=40$| M|$_{\odot }$|

In this work, we use for the ejected mass from stars the stellar yield sets from Cristallo et al. (2011, 2015), the fruity set, for low- and intermediate-mass stars, and from Limongi & Chieffi (2018) for massive stars. The first ones produced a first set of stellar yields for stars between 1 and 3 M|$_{\odot }$| (Cristallo et al. 2011), which is complemented with other set for masses in the range from 4 to 6 M|$_{\odot }$| (Cristallo et al., 2015). They follow carefully the evolution of each star, particularly the AGB (asymptotic giant branch) phases with successive thermal pulses and dredge-ups. They give in their fruity web page (http://fruity.oa-teramo.inaf.it/) the complete information, stellar masses and core and surface abundances in each step, as well as the net and total yields for 400 isotopes and 10 metallicities. We have used their so-called total yields, which actually are the total ejected mass (new and already existing in the star when formed) of each isotope to include in our code, and this way to obtain the Q’s matrices for the eight stellar masses given by the authors. Moreover, using the information of the successive steps of the AGB stars evolution, that is, each remaining mass and the corresponding surface abundance, we have computed the secondary and primary components for |$^{14}$|N and |$^{13}$|C, by assuming that the production during the third dredge-up may be considered as primary and the rest is secondary. For massive stars, Limongi & Chieffi (2018) give the stellar yields for four metallicities ([Fe/H] = 0, −1, −2, and −3) and for nine stellar masses, from 13 to 120  M|$_{\odot }$| (http://orfeo.iaps.inaf.it). They calculated the same sets for three different values of the stellar rotation velocity: 0, 150-, and 300 km s|$^{-1}$|⁠, finding, in agreement with previous works (Meynet, Ekström & Maeder 2006; Chiappini et al. 2011), that the rotation modifies the stellar yields, mainly the one for |$^{14}$|N, which appears as primary for low metallicities. In fact, since not all stars rotate at the same velocity, it is necessary to use a distribution of them for each metallicity. Following the recommendations of these authors, we have used the distribution shown in fig. 4 from Prantzos et al. (2018). Therefore, we have finally a mix of stellar yields at each rotation velocity for each one of the four metallicities. Then, we have interpolated at the same 10 metallicities as fruity for low- and intermediate-mass stars, to obtain a set of ejected masses for the whole mass range from 1 to 120  M|$_{\odot }$| and 10 metallicities. The stellar yields given by the authors for |$^{24}$|Mg, |$^{28}$|Si, and |$^{40}$|Ca have been multiplied by a factor of 2, 0.7, and 1.6, respectively, to better reproduce the observations. Finally, we added the stellar yields for SNe-Ia, see Section 3.3. These new yields were already included in Cavichia, Mollá & Bazán (2023), successfully reproducing the Galactic bulge abundances ratios and their dependence along the evolutionary time. In this work, we use starmatrix (Bazán & Mollá, 2022), an Open Source python code⁴ that compute the new ejected elements by each SSP, then included in the MulChem model by a convolution with the SFR.

3.2 Supernova rates and delay time distributions for SNe Ia

We consider that all stars with |$m \gt 8\, \mathrm{ and}\, \le 40$| M|$_{\odot }$| will end their life as CC SNe (including all types together), although the limit of mass of stars that explode as CC SNe is under debate (Doherty et al., 2017). For an SSP, the total number of such events, |$N_{\mathrm{ CC}}$|⁠, per unit stellar mass formed is determined by the IMF |$\phi (m)$| as:

which gives a value of |$\sim 0.0052$| M|$_{\odot }$||$^{-1}$| for a KRO IMF, defined within the range 0.15 and 40 M|$_{\odot }$| (⁠|$\sim 0.0054$| M|$_{\odot }$||$^{-1}$| if m|$_{\mathrm{ up}}=100$| M|$_{\odot }$|⁠), and 0.0083 M|$_{\odot }$||$^{-1}$| for a Salpeter (1955) IMF, within the same mass range. The IMF from Chabrier (2003) produces two times the number of the CC SNe created compared with a KRO IMF.

For an arbitrary star formation history, |$\Psi (t)$|⁠, the instantaneous rate of CC SNe is given by:

where |$\phi ^{\prime }(m)$| denotes the adopted IMF subtracting the fraction of stars in binary systems that would yield SN Ia.

In turn, the rate of SNe Ia explosions is given by the equation:

which involves the total number of stars

created per unit time at time |$t-\tau$|⁠, the fraction of them (⁠|$A_{Ia}$|⁠) that eventually produce an SN Ia, and the delay time distribution |$\mathrm{ DTD}(\tau)$| in terms of the age |$\tau$|⁠.

If the IMF is normalized in mass, |$\int _{m_{\rm low}}^{m_{\rm up}} m\ \phi (m)\ \mathrm{ d}m=1$|⁠, as usual, then |$N_{*} = \int _{m_{\rm low}}^{m_{\rm up}} \phi (m)\ \mathrm{ d}m$| is the number of stars per unit stellar mass of a generation or single stellar population. This value depends again on the IMF, being 1.71 M|$_{\odot }$||$^{-1}$| for the KRO IMF used here (2.02 M|$_{\odot }$||$^{-1}$| for Salpeter). In that case, |$N_{\mathrm{ Ia}}=N_{*}\times A_{\mathrm{ Ia}}$| gives the number of SN Ia produced by each stellar generation. Using observational data of SN Ia rates in galaxies, Greggio (2005) found a value |$A_{\mathrm{ Ia}\sim 10^{-3}}$|⁠, although this value may vary (Bonaparte et al. 2013).

This way, equation (5) may be written as:

The function |$\mathrm{ DTD}(\tau)$| is the delay time distribution that describes how many type Ia SN progenitors die after a delay time |$\tau$| for a SSP of mass 1 M|$_{\odot }$|⁠. This function is usually normalized to 1 when integrated along the age of a galaxy:⁵

where |$\tau _{i}$| is the minimum delay, or minimum lifetime for one star precursor of an SN Ia. Since it is usually assumed that the SNe Ia occur in binary systems where the mass limits are assumed to be between 3 and 16 M|$_{\odot }$|⁠, |$\tau _i$| would correspond to the age of the most massive component possible, that is 8 M|$_{\odot }$|⁠. For the upper limit, we need to choose between the longer age corresponding to the lower mass that could exist in a binary system, |$\sim$| 1.5 M|$_{\odot }$| and, if the conditions allow it, the time in which the explosion may occur. In practice, this |$\tau _{x}$| may be as long as the Hubble time.

Different progenitor scenarios for SNe Ia, and therefore different possible DTDs are at present supported by data. The MulChem model in its original form (Ferrini et al. 1992, 1994) used the technique explained in Greggio & Renzini (1983) to compute the SNe Ia rate from a secondary mass distribution. In more recent versions, from Mollá & Díaz (2005) until MOL19, the code incorporated the rates given by Ruiz-Lapuente, Canal & Burkert (1997) and Ruiz-Lapuente et al. (2000), who provided a numerical table (private communication) with the time evolution of the SNRs for an SSP, computed under different assumptions or scenarios and their probabilities of occurrence. The final function, called RLPPC00, composed by several Gaussian functions, each one describing a given channel and showing a maximum at a different delay time, is shown in Fig. 4(a). In addition to the RLPPC00 DTD, which has been consistently used for the verification of our results, we have employed 14 additional DTDs, as listed in Table 2. These were obtained for both DD and SD scenarios by various authors, following either theoretical (Greggio 2005) or observational prescriptions.

As explained in Introduction, there also exist some empirical works that deduce the SNR or DTD necessary to obtain the observed SN frequency. A review about this technique is given in Maoz et al. (2011). From an observational point of view, the first DTD studies from Mannucci et al. (2005) and Mannucci, Della Valle & Panagia (2006) established a relationship between the SNR and the colour and Hubble type of the SN host galaxy, both measurements considered approximations to the SFR. After these first attempts, Sullivan et al. (2006), Li et al. (2011a, b), and Smith et al. (2012) found two populations of SNe Ia: one associated to the SFR with short delays of 100–500 Myr, the prompt population; and a tardy population associated to the galaxy mass and delays of up to 5 Gyr. However, more recent studies have found a general relation between the SNR and the age of the underlying stellar population, implying a continuous DTD proportional to |$t^{-1}$| (Maoz & Mannucci 2012), a result supported by theoretical models for DD scenarios.

Measuring the SNR in galaxy clusters, a DTD also proportional to |$t^{-1}$|⁠, but with delays longer than 2 Gyr up to 10 Gyr, was found (Maoz, Sharon & Gal-Yam 2010; Barbary et al. 2012). The same result was obtained for field elliptical galaxies (Totani et al. 2008). Aubourg et al. (2008), Raskin et al. (2009), and Thomson & Chary (2011) also establish the existence of a prompt population that can dominate the SNR, even if the current SFR is low compared to the present value in the MWG (Mannucci 2009; Maoz & Mannucci 2012). More recently some works computed the DTD from the volumetric rate evolution of SN Ia. Large SN surveys are then necessary to infer the SNR, while the SFR comes from cosmological simulations (Graur & Maoz 2013; Rodney et al. 2014; Friedmann & Maoz 2018), from the modelling of the colours for the galaxy (Heringer et al. 2017; Heringer, Pritchet & van Kerkwijk 2019), from the use of stellar population synthesis codes for individual galaxies (Maoz, Mannucci & Brandt 2012), or even from integral field spectroscopy analysis of the host galaxies (Castrillo et al. 2021).

This way, some power-law functions that vary with time from Maoz & Graur (2017), Castrillo et al. (2021), and Chen et al. (2021), and are valid for the DD scenario, are used here:

where A, B, and C are normalization constants. We have also added another function, STROLG5, as given by Strolger et al. (2020, their fig. 5, right panel), obtained with parameters given in that figure. All these functions are represented in Fig. 4(b), with different colour lines, as labelled. In most of the empirical works, the exponent of the function is in the range |$\sim -1.1$| to |$-1.35$|⁠: |$-1.1$| for MAOZ017, |$-1.30$| for Friedmann & Maoz (2018), |$-1.34$| for Heringer et al. (2019), and |$-1.1$| for CASTRIL; while the delay time |$\Delta \tau$|⁠, as given in equation (9), is more variable: 40 Myr for CASTRIL, 50 Myr for MAOZ017, and 120 Myr for CHE2021. These latter authors established that although their results indicate a high confidence of |$\Delta \tau = 120$| Myr, and their slope is |$-1$| consistent with the DD scenario, they cannot firmly discard other channels of SNe Ia with other times of appearance.

Some theoretical functions from Greggio (2005) are also used and included in Fig. 4(b) (GRCDD04 and GRCDD10, and GRWDD04 and GRWDD10), which correspond to the close and wide DD scenarios with time parameters 0.4 and 1.0 Gyr. They are compatible with models represented by power-law functions for times longer than some hundreds of Myr, except the one from STROLG5.

Moreover, we have also tried other DTDs that would correspond to the SD scenario. In that case, we have used the curves given by Greggio (2005) for her models SD with MCh and with sub-MCh (GRESDCH and GRSDSCH). To these theoretical prescriptions, we have added other four functions as given by Strolger et al. (2020, their fig. 5, left panel), obtained by modifying the parameters of their equation to fit the observational data from Nelemans et al. (2001). All these functions and the one from RLPPC00, which is a mix of channels, are drawn Fig. 4(c) with different colour lines as labelled.

Some years ago, Mennekens et al. (2010) performed synthesis of populations to produce models of SD and DD scenarios. They compare the resulting DTD curves with observational data (Mannucci et al. 2005; Totani et al. 2008), as seen in their fig. 2. They show that the SD scenario produces a maximum in 1.4 Gyr and then disappears after 8 Gyr. The DD scenario creates a curve proportional to |$t^{-1}$| and it is the dominant channel. However, the SD alone does not fit the observed distributions, being DD or a mix of both channels SD|$+$|DD the most probable scenario. In fact, in that figure they presented a mix of both channels that reproduce well the data. Similarly, Ruiter et al. (2009, 2011) did similar models for four scenarios: DD, SD and HeR with MCh, and one DDet sub-MCh scenario involving helium-rich donors. Computing the DTDs for the four cases and comparing with the observations, they show that the sub-MCh may create both, a prompt (⁠|$\lt 500$| Myr) channel and a more delayed (⁠|$\gt 500$| Myr) one.

Such as Nelemans, Toonen & Bours (2013) claimed, ‘the observed SN populations probably originated from a combination of multiple channels’. That is, although the SNe Ia scenario defines mainly the delay time for the explosion will occur, a diversity of models is possible due to the different mass accretion rates, the existence or not of winds, the stripping rate of the companion, or the variations to radiate away the angular momentum, which change the time necessary to have a explosion in the DD scenario. Thus, theoretical models give delay times between 0.5 to 3 Gyr for SD scenarios and a mean time larger than 0.3 Gyr for DD scenarios. In fact, there is not a unique model able to reproduce satisfactorily the diversity of observations of SNe Ia. It is, therefore, quite probable that several channels are mixed to produce the SN Ia observed distribution (Mernier et al. 2016; Hitomi Collaboration 2017). Although we have mainly used only one DTD in each model, except when using RLPPC00 that is by itself a mix of scenarios, it is necessary to take into account that the most probable solution will need two or more DTDs, either simultaneously or varying along the time, in the same model.

3.3 The SN Ia element yields

As explained above, there is a variety of possible thermonuclear explosion models, which produce elements in different quantities (Arcones & Thielemann, 2023). The explosion mechanisms may be due to different situations, which are widely explained in Liu et al. (2023, see their table 1 and fig. 3) and summarized as: (1) the burning of C in the centre of the star, when the WD reaches the MCh (the classical case): the simmering phase starts, with the consequent compression heating and a central burning that produces the thermonuclear runaway; (2) the burning of C in the outer envelope, outer ignition that occurs with a merger of two WD (DD scenario); and (3) the burning of He, which starts in the surface of the star, when He is accreted from a companion rich in helium (He-WD or CO-WD with He envelope). This burning provokes a shock wave towards the centre of the WD and produces a WD detonation.

There exist differences in the described mechanisms, since the flame propagation may occur with different velocities: (i) detonation, when the burning occurs at supersonic velocities; (ii) deflagration, when the burning occurs at subsonic velocity; and (iii) when a deflagration to detonation transition (DDT) occurs.

For this particular subsection we are interested in works where the nucleosynthesis yield sets are calculated for the different models. We refer to Table 3 for the details about the computed models in this work (Iwamoto et al. 1999; Seitenzahl et al. 2013a; Leung & Nomoto 2018, 2020; Mori et al. 2018; Bravo et al. 2019; Gronow et al. 2021). These models are divided into two categories: (i) those where the explosion occurs at MCh (or near-Mch), and (ii) those where it occurs at sub-MCh. The first ones are usually computed for a DDT channel in an SD or DD scenario; while the second ones are calculated for the DDet case in an SD scenario, or for the merging of two WD (violent mergers) in a DD scenario. Therefore, although a priori there is no direct correlation between the SN Ia scenarios and the explosion mechanisms, there are some limitations: the SD scenario cannot produce an outer ignition, while triple systems occur basically with outer C ignition. The intermediate case CD scenario does not occur with surface He-burning. But both scenarios SD and DD may occur with any mode of explosion. The key ingredient to obtain a SN Ia is the density and temperature in the burning time. This determines the ejected quantity of |$^{56}$|Ni.

After the seminal work from Nomoto (1982) with the classical W7 model, the revised model by Iwamoto et al. (1999) is the most frequently adopted in chemical evolution studies. The W7 model adopts a pure deflagration schema in a DD scenario and was conceived to reproduce the observational data, so that there are inherent physical limitations and additional 1D modelling limitations. The yields were calculated for solar metallicity and a 1/10 of the solar value. Iwamoto et al. (1999) also studied a slow delayed detonation with DDT and a channel CD.

Later, Seitenzahl et al. (2013a) presented results for a suite of 3D, high-resolution hydrodynamical simulations of delayed-detonation explosions (DDT) giving four data sets for |${\rm [Fe/H]} = -2$|⁠, −1, −0.3-, and 0. In this model, the detonation is delayed with three phases: (1) there is a deflagration with burning in a high-density core; (2) the energy propagates with subsonic velocity and the star expands, so the density is now lower; and (3) a detonation occurs after a time. Seitenzahl et al. (2013b) found that the yields show a great dependence on metallicity in both CC SN and SNe Ia, finding that the 50  per cent of SNe Ia would have explode as near-MCh WDs to reproduce the observations related to [Mn/Fe].

In recent years, there has been a growing number of studies that have computed models of explosions, providing the yields of elements generated by them. We have selected some of these works to utilize in our research (Leung & Nomoto 2018; Mori et al. 2018; Gronow et al. 2021). These studies employ various thermodynamic conditions, such as central densities, velocities, turbulent velocity fluctuations, and different strengths of the deflagration phase. Additionally, these models also allow for explosions with masses below the sub-MCh. Also, there are two sets from Bravo et al. (2019) giving yields able to reproduce the observations of SN Ia. In turn, Leung & Nomoto (2020) have revised the W7 model for an MCh deflagration explosion with improvements on the nuclear reaction network in comparison with those from Iwamoto et al. (1999), also included in our set of yields.

The yields of different elements ejected used in this work correspond to 12 different sets obtained for different explosion mechanisms and given in Table 3, where we define for each set a number (column 1) and an acronym or name of the yield set (column 2) that we will use later for identifying each one. The third column displays the type of explosion, the fourth column the chosen model or table, the fifth column the source reference, and the sixth column the line style to plot the models in the next figures. In some cases, the authors give more than one table and we have selected one or more sets following their recommendations. All of them are plotted in Fig. 5. In panel (a), we represent the ejected mass of each element as a function of the number of the yield set (column 1 of Table 3); in panel (b), we represent the abundance compared with the solar one from Lodders (2019) as a function of the mass number A of each nuclei. Each set is drawn with a different colour as labelled. The yield for Fe, the most important element in this case, is essentially the same for all of them. For O, except GR20211, Si, and S, all tables give results within an order of magnitude. There are, however, important variations in the yields for |$^{24}$|Mg and |$^{40}$|Ca among the different sets, where models GR20211, LN20182, and BR20192 show the lowest yields values.

4 RESULTS

We have computed 180 models for MWG with 15 different possibilities of DTD functions, representing different scenarios for the binary systems and their evolution, as shown in Table 2, combined with twelve different set of yields for different explosion mechanisms for the SNe Ia, as shown in Table 3. As it is said in the Introduction, not all combinations of DTD-SN yields are compatible with our knowledge about SN Ia scenarios and channels. Thus, yields computed for DDT models with MCh explosions are compatible with DTD computed for SD scenarios, but the DDT explosion is also possible in the DD scenario. In turn, yields for sub-MCh are calculated for DDet channels in an SD scenario. The merging of two WDs as violent mergers may also occur at sub-MCh. However, we have not used any set of yields valid for this case, and this way, yields sets for sub-MCh are here only valid for the DDet explosions in SD scenarios. The combinations sub-MCh yields |$+$| DD DTDs that we have among the 180 computed models are not realistic in this work and these combinations must be interpreted with caution.

Our computations are based on the same model from MOL19, with the same total mass radial distribution, collapse time-scales, and star formation efficiency for MWG. We have, however, updated the stellar yields to the those from Cristallo et al. (2011, 2015) and Limongi & Chieffi (2018), and we have also modified the stellar mean lifetimes. For that reason, we have checked in Appendix  B (Supporting Information), that models give, for the quantities not related with Fe or |$\alpha$|-element abundances, the same results than in previous works, reproducing equally well the MWG data: SFR and the time evolution of CNO abundances for the solar region, and the radial distributions of SFR, mass of gas and stars, and CNO abundances in the Galactic disc, all of them compiled in Mollá et al. (2015, MOL15).

4.1 Supernova rates in the MWG

Once our models are well calibrated, we proceed to analyse the results related to the SNRs. Panel (a) of Fig. 6 shows the evolution of both the CC SN and SN Ia rate as a function of time. Since the same IMF is used for all models, the CC SNR is basically the same for all of them, as shown by the (overlapping) lines located in the upper side of the panel (a). In contrast, the SN Ia rate evolves differently in each model, according to the used DTD. Thus, SN Ia appear at a different time for each DTD (see delay times in Table 2), although it is difficult to see it at the scale of the age of the Universe. The first to appear is the STROLG1 (lavender line), followed by GRSDSCH, CASTRIL, and MAOZ017, while the last ones are RLPPCC00, STROLG3, and STROLG4. In panel (b), it is shown how the ratio between both types of SN changes among models. The lowest value at the present time corresponds to STROLG4 (0.78) and the highest one is for GRSDSCH (1.09). In the same panel, we also represent the ratio |$N_{\rm CC}/N_{\rm Ia}$| from Maoz & Graur (2017), as a green square, and the ratio around |$\sim 5\text{---}7$| (depending on the IMF) advocated by Greggio (2005) for an SSP as a blue triangle.

From an observational point of view, this ratio is dependent on the selected galaxy sample, varying with morphological type and |$B-V$| colour (Li et al., 2011b). Being affected by the star formation history, we do find that it changes according to the radial region, although we only represent here the solar region. Following Mannucci et al. (2005), Li et al. (2011a), and Maoz et al. (2014), the observed ratio between CC and Ia SNe may be as small as 1 (or even lower than 1) or as large as |$\sim 200$|⁠. This range is illustrated by a red dot with error bars. Thus, all used DTDs give equally acceptable results for the present time within the range defined by the error bars.

4.2 The time evolution of metallicity and alpha-over-iron relative abundance on the solar vicinity

In this subsection the results related with the time evolution of the abundance for Iron [Fe/H] and |$\alpha$|-elements [X/Fe] for the 148 models considered as realistic are presented. In Fig. 7, we display the age–metallicity relation for the solar region, the metallicity being represented by [Fe/H]. Some combinations reproduce the data better than others, however, it is very difficult to choose the best model, given the high dispersion of data.

The next step is to analyse the evolution of different |$\alpha$|-elements relative abundances [X/Fe], which incorporate both CC and Ia SNe ejecta. In Fig. 8, we present the relative abundances [X/Fe] for elements O, Mg, Si, and Ca as a function of the stellar age in Gyr. In all plots, each set of a different colours shows models for a given DTD (with the same colour meaning as in previous figures). As before, the different yields from SN Ia have a smaller effect over the plots compared with the DTD and the results show a dispersion around each colour/DTD. From the solar age (dashed vertical line) to the present age, all models follow a similar behaviour with differences only in the absolute abundances for O, Mg, and Si, while Ca shows quite important differences in this last evolution times for some particular SN Ia yield sets. For Ca, LEU2020 gives results well down of the data points from the solar birth time until the present time for all DTDs.

Since the elements produced by massive stars appear first in the ISM, the relative abundances [X/Fe] are high compared to the solar values for early times/old stellar ages. Later, the abundances of the iron-peak elements, produced by SNe Ia, start to increase in the ISM and the relative abundances decrease. From an evolutionary point of view, the relative abundances [X/Fe] begin to decrease just when the first SN Ia explode, and, obviously the time when this occurs is very dependent on the DTD. Following these plots, models using CASTRIL and STROLG1 are the closest ones to the data for ages larger than 10 Gyr (early times). This is expected since these DTDs produce SN Ia very promptly once binary systems are created, as indicated by the short delay times |$\Delta t = 40$| Myr in the first case, and by the parameter |$\xi=10$| Myr in the second one. The STROLG4 models – light blue lines – DTD show a bump for all elements at early times (old ages) that does not exist in the data trends doing this DTD invalid for our purposes.

Fig. 9 displays the same relative abundances [X/Fe] but now as a function of [Fe/H], following the same structure as before. To use [Fe/H] as a proxy of time has the advantage of eliminating the stellar age uncertainties, since only abundances are involved. Kobayashi, Leung & Nomoto (2020, see their figs 17 and 18), considering an MCh case and a sub-MCh case – associated to SD and DD scenarios, respectively – claim that in the plot of [O/Fe] versus [Fe/H] the difference between SD-DTDs and DD-DTDs would appear as a flatter line in the lower metallicities with a strong decline for around solar metallicities in the first case, while for the second case the line would decline smoothly from the lowest metallicities. We see, however, that not only these two behaviours appear for our set of models with different DTDs, and different slopes are observed. The flatter line is the one from STROLG4 – light blue lines –, while the steeper ones correspond to STROLG1, STROLG2, and STROLG3, – lavender, blue, and cyan lines –, all of them corresponding to an SD mechanism. The plateau also appears for STROLG5, GRWDD04 and RLPPC00 with a DD scenario. These results are in agreement with the ones obtained by Vincenzo, Matteucci & Spitoni (2017) and, more recently, by Palicio et al. (2023). Moreover, Kobayashi et al. (2020), by computing models with some different proportions of both scenarios, obtained that a combination with a 25 per cent sub-MCh SNe Ia and a 75 per cent MCh SNe Ia reproduces better the observational constraints that only one scenario.

Not all models are broadly compatible with the observational data. In fact, it seems difficult to find – at naked eye – a model able to reproduce simultaneously the data in all panels (see the next sections). Thus, the best models for Ca seem to be the ones in the lower part of plots (STROLG1 models, lavender lines), while O, Mg, and Si data are better fitted with STROLG5 (coral lines), and other DTDs corresponding to DD scenarios as GRWDD04 (orange lines), or RLPPC00 models (magenta lines), but also by STROLG2 and STROLG3 (blue and cyan lines), or GRESDCH (violet lines), which corresponds to an SD scenario. For Ca, the two sets of yields that in Fig. 8 show values above and below the observational data at lower ages, also show in Fig. 9 the same disagreement with observations.

We see in Figs 7 and 9 that the models do not reach [Fe/H] > 0.2 dex as observed in the stellar data. Since the models show the enrichment only for the the solar region, there are two possibilities to explain stars with metallicity higher than 0.2 dex: (i) they were born at other radial regions and then migrate to the solar region; and (ii) the errors in the metallicities and the stellar ages for these stars are higher than the estimated ones. A more detailed investigation will be given in Section 4.4 to support the first possibility.

4.3 Best models selection

To select the best model capable of reproducing all our data simultaneously, we performed a maximum-likelihood estimation (see e.g. Lehmann & Casella, 1998). We use the data described in Section 2 and compare them with the corresponding results in each model. To perform the calculation, we need to define the probability of an observation given a model. We first find the closest model point to the observational point of interest. The distance is taken in the space of all measurements available for a given observation, in our case the abundances and the age of each star. Assuming that we have uncorrelated observations with normally distributed errors in the measurements, the probability of an observation for a given model is then given by:

where |$M_{\mathrm{ nmod}}$| denotes a given model in the range |$\mathrm{ nmod} \in [1,148]$|⁠, m is a model point in |$M_{\mathrm{ nmod}}$|⁠, and |$X_{\mathrm{ obs}}$| represents each observational point i that has j measurements (abundances and age) with |$\sigma _{i,j}$| their respective errors.

We calculate for each defined model |$M_{\mathrm{ nmod}}$| the likelihood using the following expression:

Numerically, it is convenient to work with the log-likelihood and that is the procedure we adopt. With the calculated likelihood for each model on the grid, we can obtain the one that best explains the available data by finding the model that gives the maximum value. To obtain the models inside the 1|$\sigma$| confidence level, as we are assuming a flat prior for all of them, we added the values of the highest likelihood models that correspond to 39.3 per cent of the volume resulting from the cumulative distribution function. These results are summarized in Fig. 10, where we show the probability distribution obtained based on the likelihood of each model in the range |$\mathrm{ nmod} \in [1,148]$|⁠, with the models inside the |$1\sigma$| confidence interval marked by empty black circles and the best model location marked by the red solid circle. Our analysis has shown that only at the 1|$\sigma$| confidence level are there relevant differences between the models, such that this confidence interval will be considered in subsequent analyses.

In Table 4, we give the best 52 models that at 1|$\sigma$| confidence level are able to reproduce simultaneously the observational data. In this table, column 1 displays the DTD name, column 2 the type of scenario, column 3 the SN-Ia yields set, and column 4 the explosion type. Column 5 displays the corresponding normalized probability associated with the maximum-likelihood procedure. We include in Table 4, column 6, the flag Q that gives the quality of the different Yield|$+$|DTD scenario compatibility, with six categories:

(A) The DTD and the yield are in good agreement from a theoretical perspective, corresponding to the classic SD-MCh-DDT scenario.

(B) These model combinations are also in good agreement from the theoretical perspective, although they correspond to the SD-sub-Mch-DDet scenario.

(C) Models where the yield was calculated by a parametrization associated to the SD-DDT scenario, like high central densities for the WD, while the DTD model came from a DD scenario.

(D) This case refers to models whose DTD were derived from observational analysis. For these models, the yields correspond to different calculations of the same classic SD-MCh-DDT scenario, which may contradict the DD time dependence (⁠|$t^{-1}$|⁠) of the DTD models.

(E) Models whose yields are for a double-detonation mechanism at sub-MCh, but their associated DTD do not include the sub-MCh case as a possible SNIa progenitor. Although the yield and the DTD are both consistent with an SD scenario, in some combinations, the DTD models assume a hydrogen donor while the yields needs a helium donor star.

(F) Models where the yield corresponds to the classic Mch-DDT case, but the DTD assumes the SD scenario with a sub-Mch explosion mechanism.

Categories A and B do not have any incompatibility between the yield and DTD, while the remaining present some warnings, as described. It is interesting to point out that the four best models in Table 4 correspond to the combination DTD function and SN Ia yield that is in agreement with the classic SD and Mch scenario. Within the 52 models which may reproduce all the data plots simultaneously at the 1|$\sigma$| confidence level, the results point to a higher predominance of lower N|$_{\mathrm{ yields}, \mathrm{ SN-Ia}}$|⁠, indicating that models with yields that include MCh or Near MCh (classical scenario) correspond to 75 per cent and are most probable to fit the data. Analysing the DTDs distribution within the 52 most probable models, we see that 40 per cent are in the DD scenario, while the remaining 60  per cent correspond to the SD scenario. We also noticed that in Fig. 10 there are three regions that concentrate the highest probabilities, approximately located in the intervals:

• N|$_{\rm yields, SN-Ia} \in [2,5]$| and N|$_{\rm DTDs} \in [2,4]$|⁠.

• N|$_{\rm yields, SN-Ia} \in [2,5]$| and N|$_{\rm DTDs} \in [9,13]$|⁠.

• N|$_{\rm yields, SN-Ia} \in [10,12]$| and N|$_{\rm DTDs} \in [10,12]$|⁠.

The results of our multidimensional analysis independently confirm previous results in the literature, as those from Kobayashi et al. (2020), where by the analysis of [O/Fe] versus [Fe/H] and considering only the MCh case and a sub-MCh cases, it is found that a combination with a 25 per cent sub-MCh SNe Ia and a 75 per cent MCh SNe Ia reasonably reproduces the observational constraints.

At the 1|$\sigma$| confidence level, we can discard the SN Ia yields BR20191, LN20182, and IW1999 (MCh explosions) and the DTD function STROLG4 (SD), that for the whole set of combinations of SN Ia yields and DTDs functions proposed in this work, none has satisfactory results for them.

The results of the 52 best models are displayed in Figs 11, 12, and 13, for [Fe/H] versus age, [X/Fe] versus age, and [X/Fe] versus [Fe/H], respectively. We also plotted in these figures the best model selected as the highest probability model, given by the combination N|$_{\mathrm{ yields}, \mathrm{ SN-Ia}}= 5$| (LN20181, Near MCh scenario) and N|$_{\mathrm{ DTD}}=11$| (STROLG1, SD scenario). Considering the age evolution (Figs 11 and 12), the best models are a compromise between fitting the data at large ages values and also the high density data at |$\sim 3$| Gyr. For Ca, a fraction of the best models fits the higher ages data and the other fraction the high density data at |$\sim 3$| Gyr. Considering the abundances ratios as a function of [Fe/H] (Fig. 13), most of the best models pass through the high density data at [Fe/H] |$\sim 0.0$| dex, with the exception of some models for Mg and Ca. By analysing these figures, some part of the data dispersion in the [X/Fe] ratios may be explained by the mixing of different SN Ia production channels and yields. This will be addressed in Section 4.4.

Ultimately, the combination of SN Ia yields and DTD functions in Table 4 can be used as reference for future chemical evolution models that intend to reproduce the chemical abundances in the Galactic disc.

4.4 Results for other radial regions and abundances dispersions

As obtained in Section 4.3, 52 best models are able to reproduce all observational data sets simultaneously, considering the 1|$\sigma$| confidence level. It is, therefore, possible that some channels (DTDs or explosion mechanisms) appear simultaneously considering all SN Ia events in a given region. In that case, they would also create some data dispersion. As we noted previously, the relative abundance [X/Fe] data for the |$\alpha$|-elements show a large dispersion that any chemical evolution model may not reproduce with only one region. In order to seek for a possible cause of the high dispersion in the observational data, we look at the dispersion that the 52 best models produce when other radial regions are included in the plot.

The age–metallicity relation is represented in Fig. 14 (top panel). There is a clear dependence on the Galactocentric distance, the inner regions showing higher metallicities than the outer ones for the same age. This results from the negative radial gradient of metallicity presented at the Galactic disc and this dependency is stronger in the past than now. In Fig. 14 (middle panel), we represent the evolution of [O/Fe] as a function of the stellar age for the 52 best models and including different radial regions of the disc, as in the top panel. In the same figure, the bottom panel shows the relation [O/Fe] versus [Fe/H]. The binned data are shown by the black dots with error bars and correspond to the mean and the standard deviation of the data, respectively, weighted by the respective errors. The fact that different radial regions show variable [Fe/H] and [O/Fe] is due to the SFR changing on radius in the disc, with higher rates in the early times for the inner regions, while in the outer regions the SFR starts with low values increasing with time.

Since different regions of the disc do show different patterns in the plane [Fe/H] versus age and in the plane [X/Fe], either versus age or versus [Fe/H], the observed dispersion would be an indicator of mixing of stars originally born in different radial regions than the solar neighbourhood. Additionally, the results point to an inner radial migration towards the solar region. The different DTDs or SN Ia yields sets may also contribute to the dispersion. These results agree with Palla et al. (2022) findings, who by using [Mg/Fe] abundances from AMBRE-HARPS, and comparing them with detailed chemical evolution models, conclude that none set of yields may reproduce the data, even taken into account the stellar migration. We note that Buder et al. (2019) also find evidence for radial migration in the GALAH data at the solar vicinity, where stars older than 8 Gyr have angular momenta lower than the Sun, which implies that they are in eccentric orbits and are originated in the inner disc.

5 CONCLUSIONS

In this paper, we have examined the role of SNe Ia, focusing primarily on the different DTDs determined by the scenarios of binary systems, as well as the yields of elements created by different explosion mechanisms. The findings of this study can be summarized as follows:

• We have analysed carefully a large set of observations related with the |$\alpha$|-element over iron abundances obtained from some modern surveys and other data of the literature. The detailed description of these data can be seen at Appendix  A.

• An updated version of the MulChem chemical evolution model applied to the MWG was used, where we have now taken the stellar yields from Cristallo et al. (2011, 2015) (the fruity set), for low- and intermediate-mass stars, and from Limongi & Chieffi (2018) for massive stars with several rotation stellar velocities. The calibration of the MWG model by using these stellar yields is given in Appendix  B.

• To these prescriptions, we have added 15 possible DTDs for both DD and SD scenarios, and 12 sets of SNe Ia yields from recent works with different mechanisms of explosion. The combination of both inputs provide 180 different models, 148 of them considered realistic enough, although some of them are better than others, considering the current knowledge of SNe Ia production. The models results have been compared with the observational data presented in the first part of work.

• For this comparison, a multidimensional maximum-likelihood procedure was implemented, where the abundances of O, Mg, Si, Ca, and Fe and, additionally, the stellar age are fitted simultaneously. From the analysis, 52 best models are able to reproduce the data sets simultaneously at the 1|$\sigma$| confidence level, which represent a |$\sim 35$| per cent of the realistic models (a |$\sim 29~{{\ \rm per\ cent}}$| of the total number of computed models). From these best models, the most probable DTD is the one called STROLG1 from Strolger et al. (2020), defined by the parameters |$\xi=10$| Myr, |$\omega =600$|⁠, and |$\alpha =220$|⁠, with a very prompt delay. The best set of SN Ia yields seems to be LN20181 (MCh) from Leung & Nomoto (2018). The combination STROLG1 (SD) |$+$| LN20181 could indicate that a single degenerate with a time of delay until the explosion from the merger (SD with MED) could be a real channel for the SN Ia production.

• Within the 52 best models able to reproduce all observational data sets simultaneously, considering the 1|$\sigma$| confidence level, the SN Ia yields that include MCh or Near MCh (classical scenario) correspond to 75 per cent and are most probable to fit the data. For the DTDs within the 52 most probable models, 60 per cent correspond to the SD scenario, while the remaining 40 per cent are in the DD scenario. It is, therefore, possible that some of these channels (DTDs or explosion mechanisms), taking into consideration the 52 best models, appear simultaneously considering all SN Ia events in a given region.

• The dispersion given by chemical evolution models is high for different radial regions of the disc. This might suggest that the observed dispersion is related to uncertainties in distance determinations, or with the stellar migration in the disc. In our models, differences between radial regions are due to the SFR variations. The data dispersion may be also explained by a mix of DTDs (different scenarios to create SN Ia binary systems), or by a mix of explosions channels. The observed dispersion is probably due to both effects: radial mixing and several simultaneous SN Ia scenarios/channels.

ACKNOWLEDGEMENTS

We are very grateful to the reviewer whose comments and suggestions have helped to improve the manuscript.

This work is result of the grants PGC-2018-0913741-B-C22, MDM-2015-0500, MDM-2017-0737 (Unidad de Excelencia María de Maeztu CAB), and PID2019-107408GB-C41, PID2019-107408GB-C42, PID2022-136598NB-C33, funded by MCIN/AEI/10.13039/501100011033. OC acknowledges funding support from FAPEMIG grant APQ-00915-18. LG acknowledges financial support from AGAUR, CSIC, MCIN, and AEI 10.13039/501100011033 under projects PID2020-115253GA-I00, PIE 20215AT016, CEX2020-001058-M, and 2021-SGR-01270.

This work has made use of the computing facilities available at the Laboratory of Computational Astrophysics of the Universidade Federal de Itajubá (LAC-UNIFEI). The LAC-UNIFEI is maintained with grants from CAPES, CNPq and FAPEMIG.

This work made use of the Third Data Release of the GALAH Survey (Buder et al. 2021). The GALAH Survey is based on data acquired through the Australian Astronomical Observatory, under programs: A/2013B/13 (The GALAH pilot survey); A/2014A/25, A/2015A/19, A2017A/18 (The GALAH survey phase 1); A2018A/18 (Open clusters with HERMES); A2019A/1 (Hierarchical star formation in Ori OB1); A2019A/15 (The GALAH survey phase 2); A/2015B/19, A/2016A/22, A/2016B/10, A/2017B/16, A/2018B/15 (The HERMES-TESS program); and A/2015A/3, A/2015B/1, A/2015B/19, A/2016A/22, A/2016B/12, A/2017A/14 (The HERMES K2-follow-up program). We acknowledge the traditional owners of the land on which the AAT stands, the Gamilaraay people, and pay our respects to elders past and present. This paper includes data that has been provided by AAO Data Central (datacentral.org.au).

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

DATA AVAILABILITY

• Appendix  A (Supporting Information) gives the whole description of the data and the normalization of different data sets.

• The whole set of 180 models are available in electronic format upon request to the authors.

• The DTD functions and the SN Ia yields are included in the package starmatrix, an Open Source python code (Bazán & Mollá, 2022) that compute the new ejected elements by each SSP (see https://github.com/xuanxu/starmatrix).

• The calibration of the Galactic chemical evolution model is given in Appendix  B (Supporting Information).

Footnotes

REFERENCES

APPENDIX A: NORMALIZED DATA

The description of the data sets used in this work are provided in Table 1. The main point is that these abundances are given as [X/Y] |$=\log {(X/Y)}-\log {(X/Y)_{\odot }}$|⁠, and obtained by using different solar abundances, as shown in Table A1. Therefore, we have re-normalized them to the same solar abundances scale from Lodders (2019).

After this, we have still seen differences in their scale, since we consider that all abundances [X/Fe] must be zero when the iron abundance or metallicity [Fe/H] = 0 and when the stellar age |$\tau = 4.6$| Gyr, that is, the data must be around the solar values in the solar metallicity and age. We have, therefore, analysed each cloud of data for each set and element, to search for its centre and then to shift it to locate it in the solar value.

In order to do this, we have selected the points around the solar age (3–6 Gyr), Galactocentric radius (7–9 kpc) or metallicity, in each case. We have plotted the corresponding distribution [X/Fe] and [Fe/H] by searching for the mean value by fitting a Gaussian and obtained the mean and standard deviation from the fit. As this averaged abundances should be zero, we apply shifts values |$\Delta$| [X/Fe]|$_{\rm shift}$| and |$\Delta$| [Fe/H]|$_{\rm shift}$| that correspond to the averaged values obtained from the fits. We show an example of the method in Fig. A1 with data from GALAH for [O/Fe]. In this case, we applied a shift of [O/Fe]|$\sim -0.39$| dex and <[Fe/H]>|$\sim +0.29$| dex and the result is displayed in the figure. In all cases we checked that this procedure lead to the same results as calculating the mean and the standard deviation of the sample.

We have repeated this process for [Fe/H] and all relative abundances [X/Fe]. The resulting distributions are shown in Fig. A2 for [Fe/H] as a function of age, in Figs A3–A6, for [X/Fe] as a function of the stellar age, and in Figs A7–A9, as a function of metallicity [Fe/H]. In these figures, the each dataset is represented by different symbols: green squares (Chen et al. 2000), grey stars (HARPS-GTO), and density (GALAH). For reference, in all of the figures we plot as dashed lines the Sun position. We have also plotted all sets together for each distribution in the last panel of each figure and these combined data sets are used in the statistical analysis given in Section 4.

A1 The age–metallicity relation

In this section, we show the data relative to the time evolution of the Iron abundance [Fe/H], usually taken as the measurement of the global metallicity, for the solar region. Data are taken from different surveys as given in Chen et al. (2000, CHEN), Delgado Mena et al. (2015, 2017, HARPS), and Buder et al. (2021, GALAH). Fig. A2 shows four panels with data from CHEN, GALAH, and HARPS and the bottom right panel gives the corresponding figure from all survey data sets together.

A2 The evolution of the α-element abundances over iron along the stellar age and [Fe/H]

In this subsection, we have obtained the relative abundances [X/Fe] for O, Mg, Si, and Ca as a function of the stellar age, obtained for different methods by several authors. All the data used in this work are complete in the sense that each figure has the same stars. Figs A3–A6 display the abundances ratios [O/Fe], [Mg/Fe], [Si/Fe], and [Ca/Fe] versus age in Gyr, respectively. It is clear from these figures the importance of using different databases for this study, since CHEN and HARPS add stars with different ages than GALAH, increasing the densities in the plots at higher ages.

Finally, we show the available data for the |$\alpha$|-element abundances over iron as a function of the iron abundance or metallicity [Fe/H], presented in Figs A7–A10, for the abundances ratios [O/Fe], [Mg/Fe], [Si/Fe], and [Ca/Fe] as a function of [Fe/H], respectively.

From all these figures it is possible to conform that all the databases used in this work are compatible with each other and also that they are complementary.

APPENDIX B: CALIBRATION OF MODEL FOR THE SOLAR REGION AND GALACTIC DISC

B1 The solar vicinity

In Mollá et al. (2015, hereinafter MOL15), a compilation of data from different authors (see references in MOL15) for MWG was presented. A binning process resulted in a set of binned data for each type of observation. This was described in detail in the Appendix from MOL15, where a table of references of these resulting binned distributions were given. The observational data basically included the solar neighbourhood time evolution of metallicity, SFR, and CNO abundances, in addition to the present state of the disc, including the radial distributions of surface densities for stars, gas, SFR, and CNO abundances.

In this work, we adopted the data from references described in MOL15 for the solar vicinity time evolution of SFR. The data are binned for each Gyr with |$\Delta t=\pm 0.5\, { \rm Gyr}$|⁠. The resulting SFR values were normalized to the most recent values of the SFR in the solar region at the present time, which was obtained as SFR|$_{\odot } \sim 0.266$| M|$_{\odot }$| yr|$^{-1}$|⁠, at a Galactocentric distance |$R = 8$| kpc.

The SFR and the C, N, and O time evolution for the solar vicinity are plotted in Figs B1 and B2. In these figures, all the 180 models give, as expected, similar results and are in agreement with the observational data obtained in MOL15. It is evident that these quantities do not dependent on the SN Ia rates nor on the SN Ia yields – much lower than the stellar yields from massive stars/CC SNs or from intermediate stars – shown in Fig. 5 of the main paper. Since we have updated the stellar yields for low- and intermediate-mass stars, as well as massive stars, it is important to check that these new models still reproduce well the solar region observational data. We have checked that this is true by comparing the results for the 180 models – which appear as only one – with the data. In fact, we demonstrate that the star formation history (as the C, N, and O abundances along the evolutionary time) shows good agreement compared with observational data, and effectively there is no change in the results depending on the DTD or on the SN Ia yield set.

B2 The Galaxy disc

The data for radial distributions of surface densities in the MWG disc for H i, molecular gas, stellar profile, and the SFR surface densities were taken from table A3 in MOL15. The solar stellar surface density is between 33 and 64 |$\mathrm{ M}_{\odot }$| pc|$^{-2}$| from different authors, as, for example, McMillan (2011), Bovy & Rix (2013), and Burch & Cowsik (2013). These values depend mostly on the scale length of the disc, |$R_{\rm d}$|⁠, which is in the range [2.1–5.4] kpc. Therefore, we assume a certain uncertainty of |$\sim 0.20$| dex in the absolute values of this radial distribution. The SFR binned radial distribution was also normalized to the SFR|$_{\odot }$|⁠. For what refers to the C, N, and O abundances radial distributions, we have used the data of table A4 from MOL15. In all cases, the error bars include the statistical error that corresponds to the dispersion of the original data coming from different authors, plus the systematical uncertainty due to the instrumental setup and to the techniques to estimate the final data from observations.

Next figures refer to the radial distribution of different quantities throughout the Galactic disc. We show the radial distributions of C, N, and O in Fig. B3. Fig. B4 displays the surface density radial distributions of the SFR, the stellar profile, and both gas phases, diffuse and molecular. In all figures, the corresponding data are those compiled in MOL15 with their error bars, which indicate the dispersion of the original data from different authors. In both figures, the results of all models are identical and the lines overlap each other. In all cases, the models still present satisfactory results in agreement with the observational data, taking into account the data spread and uncertainties.

We did not expect large differences between models for C, N, and O, which basically come from stars of different masses, although there are small quantities coming from SN Ia. As explained above, we have included Figs B2 and B1 as Figs B3 and B4, to demonstrate that using different sets for stellar yields than MOL15 and MOL17 give essentially the same evolution for the solar region and the Galactic disc, that is, the new model is well calibrated compared with the data for these regions.