11institutetext:  SISSA, Via Bonomea 265, I-34136 Trieste, Italy
11email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22institutetext: CNR-IFN Padova, Via Trasea 7, I-35131
Padova, Italy
22email<EMAIL_ADDRESS>33institutetext: Dipartimento di
Fisica e Astronomia, Università degli studi di Padova, Vicolo Osservatorio 3,
Padova, Italy
33email<EMAIL_ADDRESS>44institutetext: INAF-OAPD, Vicolo
dell’osservatorio 5, Padova, Italy 55institutetext: INAF-OATs, Via G. B.
Tiepolo 11, I-34143 Trieste, Italy
55email<EMAIL_ADDRESS>66institutetext: IFPU - Institute for fundamental
physics of the Universe, Via Beirut 2, 34014 Trieste, Italy 77institutetext:
INFN, Sezione di Padova, Via Marzolo 8, I–35131, Padova, Italy
88institutetext: Center for Interdisciplinary Exploration & Research in
Astrophysics (CIERA) and Department of Physics & Astronomy, Northwestern
University, Evanston, IL 60208, USA 99institutetext: INFN-Sezione di Trieste,
via Valerio 2, 34127 Trieste, Italy

# On the effects of the Initial Mass Function on Galactic chemical enrichment

S. Goswami 11 A. Slemer 223344 P. Marigo 33 A. Bressan 11 L. Silva 5566 M.
Spera 11 3 3 7 7 8 8 L. Boco 115566 V. Grisoni 1155 L. Pantoni 11 A. Lapi
11556699

(Accepted by A&A)

###### Abstract

Context.  There is mounting evidence that the stellar initial mass function
(IMF) could extend much beyond the canonical $\mbox{$M_{\rm i}$}\sim
100\,\mbox{${M}_{\odot}$}$ limit, but the impact of such hypothesis on the
chemical enrichment of galaxies still remains to be clarified.

Aims. We aim to address this question by analysing the observed abundances of
thin- and thick-disc stars in the Milky Way with chemical evolution models
that account for the contribution of very massive stars dying as pair
instability supernovae.

Methods. We built new sets of chemical yields from massive and very massive
stars up to $\mbox{$M_{\rm i}$}\sim 350\,\mbox{${M}_{\odot}$}$, by combining
the wind ejecta extracted from our hydrostatic stellar evolution models with
explosion ejecta from the literature. Using a simple chemical evolution code
we analyse the effects of adopting different yield tables by comparing
predictions against observations of stars in the solar vicinity.

Results. After several tests, we focus on the [O/Fe] ratio which best
separates the chemical patterns of the two Milky Way components. We find that
with a standard IMF, truncated at $\mbox{$M_{\rm i}$}\sim
100\,\mbox{${M}_{\odot}$}$, we can reproduce various observational constraints
for thin-disc stars, but the same IMF fails to account for the [O/Fe] ratios
of thick-disc stars. The best results are obtained by extending the IMF up to
$\mbox{$M_{\rm i}$}=350\,\mbox{${M}_{\odot}$}$ and including the chemical
ejecta of very massive stars, in the form of winds and pair instability
supernova explosions.

Conclusions. Our study indicates that PISN could have played a significant
role in shaping the chemical evolution of the Milky Way thick disc. By
including their chemical yields it is easier to reproduce not only the level
of the -enhancement but also the observed slope of thick-disc stars in the
[O/Fe] vs. [Fe/H] diagram. The bottom line is that the contribution of very
massive stars to the chemical enrichment of galaxies is potentially quite
important and should not be neglected in chemical evolution models.

###### Key Words.:

stars: abundances – stars: massive – Galaxy: abundances – Galaxy: disc –
Galaxy: solar neighborhood – Galaxy: evolution

##  1 Introduction

The chemical evolution of the Milky Way is one of the most important
astrophysical topics because it provides direct hints on the more general
question of how galaxies formed and evolved. It is also an anchor for galaxy
models because of the possibility to study their properties through the
analysis of individual stars. Indeed this field of research is continuously
growing providing on one side an increasing amount of observational abundance
data that contribute to enlarge and sharpen the whole picture (Gilmore et al.
2012; Bensby et al. 2014; Majewski et al. 2017; de Laverny et al. 2013) and,
on the other, a growing number of interpretative tools that go from simple
chemical evolutionary models to more complex chemo-hydrodynamical models (e.g.
Valentini et al. 2019). A key ingredient to interface model predictions with
observations are the stellar chemical yields, describing the contribution of
stars of different types to the metal enrichment of the galaxies. Other
physical processes play of course an important role in the chemical evolution
of the galaxies, like the functional forms adopted to describe the stellar
birthrate function, the gas inflows and outflows, the mixing of newly ejected
elements with the surrounding medium, the relative displacement of individual
stars from their original positions, etc. However stellar yields keep a key
role because, being the result of the evolution of stars, they may provide a
tight link between stellar and galactic time-scales. The contribution of
individual stars to the metal enrichment has been the subject of many studies
in the past (e.g. Matteucci 2014). The chemical enrichment from stars takes
place when elements newly produced by nuclear reactions in the deep stellar
interiors are ejected into the interstellar medium, via stellar winds or
supernova explosions.

Low- and intermediate-mass stars, $\mbox{$M_{\rm i}$}\simeq
0.8\,\mbox{${M}_{\odot}$}-6\,\mbox{${M}_{\odot}$}$, never reach the carbon
ignition temperature and enrich the interstellar medium (ISM) mainly during
the Red Giant Branch (RGB) and Asymptotic Giant Branch (AGB) phases (Marigo
2001; Cristallo et al. 2009, 2011, 2015; Karakas 2010; Karakas & Lugaro 2016;
Karakas et al. 2018; Ventura et al. 2013, 2017, 2018; Pignatari et al. 2016;
Slemer et al. 2017; Ritter et al. 2018, and references therein).

Stars in a narrow mass interval, $\mbox{$M_{\rm i}$}\simeq
6\,\mbox{${M}_{\odot}$}-10\,\mbox{${M}_{\odot}$}$, are able to burn carbon in
their non- or mildly degenerate cores and experience the so-called Super-AGB
phase (Ritossa et al. 1996; García-Berro et al. 1997; Iben et al. 1997).
Depending on the efficiency of stellar winds and the growth of the core mass,
the final fate of Super-AGB stars splits in two channels that lead either to
the formation of an O-Ne-Mg white dwarf ($7\lesssim\mbox{$M_{\rm
i}$}/\mbox{${M}_{\odot}$}\lesssim$ 9) or to an electron capture supernova
($9\lesssim\mbox{$M_{\rm i}$}/\mbox{${M}_{\odot}$}\lesssim$ 10; Hurley et al.
2000; Siess 2007; Poelarends et al. 2008).

Massive stars ($\mbox{$M_{\rm i}$}\simeq 10-120\,\mbox{${M}_{\odot}$}$)
experience more advanced nuclear burnings (Ne, O, Si) up to the formation of
an iron core, that eventually implodes producing either a successful core-
collapse supernova (CCSN) explosion, or the direct collapse into a black hole
as a failed SN (Woosley & Weaver 1995; Fryer 1999; Chieffi & Limongi 2004;
Limongi & Chieffi 2006; Nomoto et al. 2006; Fryer et al. 2006; Fryer &
Kalogera 2001; Heger & Woosley 2002; Heger et al. 2003a; Fryer et al. 2012;
Janka 2012; Ugliano et al. 2012; Ertl et al. 2016; Pignatari et al. 2016;
Ritter et al. 2018; Limongi & Chieffi 2018). On rare occasions that depend on
the characteristics of the progenitor and the details of the explosion, stars
in this mass range may give rise to powerful hypernovae (Izzo et al. 2019).

Very massive objects (VMO), $100\lesssim\mbox{$M_{\rm
i}$}/\mbox{${M}_{\odot}$}\lesssim 300\,\mbox{${M}_{\odot}$}$, enter the pair
instability regime during central oxygen burning, which may trigger their
final thermonuclear explosion (Heger & Woosley 2002; Umeda & Nomoto 2002;
Kozyreva et al. 2014b; Woosley & Heger 2015; Woosley 2017).

The chemical ejecta contributed by all stars of different initial masses and
evolutionary stages are key ingredients of galactic chemical evolution
studies. In most cases, the adopted stellar IMF is truncated at around
$\mbox{$M_{\rm i}$}\simeq 100\,\mbox{${M}_{\odot}$}$, which implies that no
chemical contribution from VMO is taken into account. This common assumption
is due to the fact that very few evolutionary models exist for stars with
$\mbox{$M_{\rm i}$}>100\,\mbox{${M}_{\odot}$}$. With the exception of zero-
metallicity stars (e.g., Haemmerlé et al. 2018; Yoon et al. 2012; Ohkubo et
al. 2009, 2006; Lawlor et al. 2008; Marigo et al. 2003, 2001), the lack of
systematic evolutionary studies of VMO has effectively limited the exploration
area of chemical evolution models, which focussed on the role of population
III stars only (e.g., Cherchneff & Dwek 2010, 2009; Rollinde et al. 2009;
Ballero et al. 2006; Matteucci & Pipino 2005; Ricotti & Ostriker 2004). In
fact, the occurrence of VMO and their final fates were thought to apply only
to primordial population III stellar populations (Bond et al. 1982, 1984;
Heger & Woosley 2002; Nomoto et al. 2013).

In recent years our understanding of the evolution of massive and very massive
stars has dramatically changed as a result of important discoveries. We can
now rely on studies focussed on young super star clusters (Evans et al. 2010;
Walborn et al. 2014; Schneider et al. 2018; Crowther 2019; Crowther et al.
2016) and on the identification of massive stellar black holes hosted in
binary systems (Abbott et al. 2016, 2020; Spera & Mapelli 2017; Spera et al.
2015). All this evidence points to an IMF that extends up to VMO, either as a
genuine top-heavy IMF, or in virtue of an early efficient merging in binary
systems (Senchyna et al. 2020). Further support for a top-heavy IMF in certain
environmental conditions comes also from recent finding that the low observed
16O/18O isotopic ratios in starburst galaxies can be reproduced only by models
that assume an excess of massive stars (Romano et al. 2017). However, while
there is strong evidence supporting the existence of stars with initial mass
up to $\mbox{$M_{\rm i}$}\simeq 350\,\mbox{${M}_{\odot}$}$, other evidence
exists in favour of the opposite scenario, i.e. that the IMF is bottom-heavy,
like in the case of observations of the gravity-sensitive narrow-band
integrated indices in local ellipticals (van Dokkum & Conroy 2012; van Dokkum
et al. 2017). These latter observations are quite challenging from the
chemical evolution point of view because massive elliptical galaxies are among
the most metal-rich stellar systems known and it is difficult to explain their
chemical pattern with such a bottom-heavy IMF (Bressan et al. 1994; Matteucci
1994; Thomas et al. 2005). In parallel, it has also been suggested that the
IMF may be a self regulating process leading to a variation of its form in
different environmental conditions. Early models based on an IMF varying with
time and with galactic environment (Padoan et al. 1997) could better explain
many of the observed properties of elliptical and low surface brightness
galaxies, such as $\alpha$-enhancement, downsizing, fundamental plane etc.
(Chiosi et al. 1998; Jimenez et al. 1998). Since then, more evidence of a
dependence of the IMF on the environment has been accumulated, leading to the
definition of galaxy wide IMF (gwIMF) and time integrated gwIMF (IGIMF) that
may deviate from a canonical IMF, depending on metallicity and star formation
activity (Kroupa & Weidner 2003; Weidner & Kroupa 2005; Kroupa 2008; Recchi et
al. 2009; Marks et al. 2012; Jeřábková et al. 2018; Hosek et al. 2019).

Motivated by these studies our group has recently carried out the first
systematic analysis of stellar models for massive stars and VMO, extending up
to $\mbox{$M_{\rm i}$}=350\,\mbox{${M}_{\odot}$}$, for a large grid of initial
metallicity, from $\mbox{$Z_{\rm i}$}=0.0001$ to $\mbox{$Z_{\rm i}$}=0.06$
(Tang et al. 2014; Chen et al. 2015; Costa et al. 2020). These models,
computed with the PARSEC code, are based on updated input physics and, above
all, rely on the most recent advances in the theory of stellar winds, which
critically affect the evolution of massive stars and VMO.

In this work we use the same PARSEC models of massive stars and VMO to derive
their stellar wind ejecta, which are then incorporated in a chemical evolution
model to analyse the impact of varying the slope of the IMF and its upper mass
limit. Additionally, we test different sets of chemical ejecta and check their
ability in reproducing the chemical characteristics of the thin and thick-disc
stars.

For this purpose we first collected from the literature various sets of
ejecta, that cover the contributions from low-, intermediate-mass and massive
stars, the latter being usually provided up to $\mbox{$M_{\rm i}$}\simeq
100\,\mbox{${M}_{\odot}$}$. Then, we constructed new tables of chemical
ejecta, extending the initial mass range to include the contribution of VMO,
up to $\mbox{$M_{\rm i}$}=350\,\mbox{${M}_{\odot}$}$. The new set of chemical
yields of massive stars and VMO contains the stellar wind ejecta, based on the
PARSEC models, suitably tailored to explosive yields for CCSN, pulsational
pair instability SN (PPISN), and of pair instability SN (PISN) taken from the
recent literature.

The structure of the paper is as follows. Basic stellar classification and
relevant mass limits are recalled in Sect. 2. In Sect. 3 we introduce our new
set of chemical ejecta for massive and very massive stars. After an outline of
our PARSEC code, we describe the method adopted to combine PARSEC stellar
evolution models of massive stars with explosive models available in the
literature. Then we discuss the resulting ejecta due to both stellar winds and
explosions (CCSN, PPISN, and PISN). The full content of the ejecta tables,
made available to the user, is detailed in Appendix A. Section 4 summarizes
the main characteristics of the AGB yields obtained with our COLIBRI code.
Section 5 introduces and compares various sets of chemical ejecta of AGB and
massive stars taken from the recent literature. Section 6 briefly describes
the observational sample of thin- and thick-disc stars in the solar vicinity
and the adopted chemical evolution code used for the interpretation of their
abundances. In Sect. 7 we analyse the predictions of chemical evolution models
calculated adopting different sets of chemical ejecta. Their performance is
tested through the use of various diagnostics, with particular focus on the
observed [O/Fe] vs. [Fe/H] diagram populated by thin- and thick-disc stars.
Finally, Sect. 8 recaps our study and its main conclusions.

##  2 Stellar classes and mass limits

The final fate of stars depends primarily on their initial mass, $M_{\rm i}$,
and metallicity, $Z_{\rm i}$. To characterise the chemical contributions of
stars it is convenient to group them in classes as a function of $M_{\rm i}$,
according to evolutionary paths and final fates. Let us introduce a few
relevant limiting masses that define each stellar family. Mass limits and
other relevant quantities used throughout the paper are also defined in Table
1.

It should be noticed that the mass ranges specified below should not be
considered as strict, but rather approximate limits, since they significantly
depend on the efficiency of processes like convective mixing and stellar winds
and, especially for massive stars, also on the initial chemical composition.

We define with $M_{\rm AGB}\simeq 6\,\mbox{${M}_{\odot}$}$ the maximum initial
mass for a star to build a highly electron-degenerate C-O core after the end
of the He-burning phase. This class comprises low- and intermediate-mass
stars, which then proceed through the AGB phase leaving a C-O WD as compact
remnant.

Stars with $\mbox{$M_{\rm i}$}>M_{\rm AGB}$ are able to burn carbon in mildly
or non-degenerate conditions. Those stars that build an electron-degenerate
O-Ne-Mg core are predicted to enter the Super-AGB phase, undergoing recurrent
He-shell flashes and powerful mass loss, similarly to the canonical AGB phase.
If stellar winds are able to strip off the entire H-rich envelope while the
core mass is still lower than $\simeq 1.38\,\mbox{${M}_{\odot}$}$, then the
evolution will end as an O-Ne-Mg WD (Nomoto 1984; Iben et al. 1997). We denote
with $M_{\rm SAGB}\simeq 7\,\mbox{${M}_{\odot}$}$ the upper mass limit of this
class of stars (Herwig 2005; Siess 2006, 2007; Doherty et al. 2014).

Stars with $\mbox{$M_{\rm i}$}>M_{\rm SAGB}$ and having an electron-degenerate
O-Ne-Mg core that is able to grow in mass up to the critical value of $\simeq
1.38\mbox{${M}_{\odot}$}$, are expected to explode as electron capture
supernovae (ECSN; Nomoto 1984; Poelarends et al. 2008; Leung et al. 2020).

Let us denote with $M_{\rm mas}$ the minimum initial mass for a star to avoid
electron degeneracy in the core after carbon burning. We note that following
this definition the progenitors of electron capture SN cover the range $M_{\rm
SAGB}<\mbox{$M_{\rm i}$}<M_{\rm mas}$.

Stars with $M_{\rm mas}\leq\mbox{$M_{\rm i}$}<M_{\rm VMO}\simeq
100\,\mbox{${M}_{\odot}$}$ are able to proceed through all hydrostatic nuclear
stages up to Si-burning, with the formation of a Fe core that eventually
undergoes a dynamical collapse triggered by electron-captures and
photodisintegrations (Woosley & Weaver 1995; Thielemann et al. 2011).

Very massive stars with $\mbox{$M_{\rm i}$}\geq M_{\rm VMO}$ may experience
electron-positron pair creation instabilities before and during oxygen
burning, with a final fate that is mainly controlled by the mass of the helium
core, ${M}_{\rm He}$ (Heger & Woosley 2002; Heger et al. 2003a; Nomoto et al.
2013; Kozyreva et al. 2014b; Woosley 2016, 2017; Leung et al. 2019), resulting
in a successful/failed CCSN or thermonuclear explosion. For further details,
we refer to Sect. 3.6.

##  3 Chemical ejecta of massive and very massive stars using PARSEC models

In this work our reference set of evolutionary tracks for massive and very
massive stars is taken from the large database of Padova and TRieste Stellar
Evolution Code
(PARSEC)111https://people.sissa.it/~sbressan/parsec.html,222http://stev.oapd.inaf.it/cgi-
bin/cmd. The PARSEC code is extensively described elsewhere (Bressan et al.
2012; Costa et al. 2019a, b) and here we will provide only a synthetic
description of the relevant input physics.

###  3.1 Evolutionary models

The PARSEC database includes stellar models with initial masses from
$8\,\mbox{${M}_{\odot}$}$ to $350\,\mbox{${M}_{\odot}$}$ and metallicity
values $\mbox{$Z_{\rm
i}$}=0.0001,0.0005,0.001,0.004,0.006,0.008,0.017,0.02,0.03$ (Chen et al. 2015;
Tang et al. 2014). The adopted reference solar abundances are taken from
Caffau et al. (2011), with a present-day solar metallicity
$Z_{\odot}=0.01524$. Note that the latter value does not correspond to the
initial metallicity of the Sun, which instead is predicted to be $Z_{\rm
protoSUN}=0.017$ (see Sect. 7.1). For all metallicities the initial chemical
composition of the models is assumed to be scaled-solar. The isotopes included
in the code are: $\mathrm{H}$, $\mathrm{D}$, $\mathrm{{}^{3}He}$,
$\mathrm{{}^{4}He}$, $\mathrm{{}^{7}Li}$, $\mathrm{{}^{7}Be}$,
$\mathrm{{}^{12,13}C}$,$\mathrm{{}^{14,15}N}$, $\mathrm{{}^{16,17,18}O}$,
$\mathrm{{}^{19}F}$, $\mathrm{{}^{20,21,22}Ne}$,$\mathrm{{}^{23}Na}$
$\mathrm{{}^{24,25,26}Mg}$, $\mathrm{{}^{26,27}Al}$, $\mathrm{{}^{28}Si}$.
Opacity tables are from Opacity Project At Livermore
(OPAL)333http://opalopacity.llnl.gov/ team (Iglesias & Rogers 1996, and
references therein) for $4.2\leq\log(T/\mathrm{K})\leq 8.7$, and from ÆSOPUS
tool444http://stev.oapd.inaf.it/aesopus (Marigo & Aringer 2009), for
$3.2\leq\log(T/\mathrm{K})\leq 4.1$. Conductive opacities are included
following Itoh et al. (2008). Neutrino energy losses by electron neutrinos are
taken from Munakata et al. (1985), Itoh & Kohyama (1983) and Haft et al.
(1994). The equation of state is from
freeeos555http://freeeos.sourceforge.net/ code version 2.2.1 by Alan W. Irwin.

The mass loss prescriptions employed in PARSEC are the law of de Jager et al.
(1988) for red super giants (RSG; $T_{\rm eff}\leq 12\,000$ K), the Vink et
al. (2000) relations for blue super giants (BSG; $T_{\rm eff}>12\,000$ K), and
Gräfener (2008) and Vink et al. (2011) during the transition phase from O-type
to Luminous Blue Variables (LBV) and RSG, and finally to Wolf Rayet (WR)
stars.

Figure 1: Pre-SN mass $M_{\rm fin}$ as function of initial $M_{\rm i}$ for different values of the initial metallicity $Z_{\rm i}$, as indicated. Stellar tracks are taken from PARSEC $\mathrm{VI.1}$ models. Table 1: Description of the main quantities used in this work. name | definition
---|---
WD |  White dwarf
ECSN |  Electron capture supernova
CCSN |  Core collapse supernova
PPISN |  Pulsation pair instability supernova
PISN |  Pair instability supernova
DBH |  Direct collapse to black hole
$Z_{\rm i}$ |  Initial metallicity.
$M_{\rm i}$ |  Mass of the star at the zero-age main sequence
$M_{\rm fin}$ |  Mass of the star at the beginning of central carbon burning (almost equivalent to the pre-SN mass)
$M_{\rm rem}$ |  Mass of the remnant
$M_{\rm cut}$ |  Mass-cut, in a pre-supernova model, enclosing the entire mass that will collapse and form the compact remnant
$M_{\rm He}$ |  Mass of the He-core at the beginning of central carbon burning
$M_{\rm CO}$ |  Mass of CO-core at the beginning of central carbon burning
$M_{\rm AGB}$ |  Maximum mass for a star to experience the AGB phase and leave a C-O WD
$M_{\rm SAGB}$ |  Maximum mass for a star to evolve through the Super-AGB and leave an O-Ne-Mg WD
$M_{\rm mas}$ |  Minimum mass for a star to experience all hydrostatic nuclear burnings up to the Si-burning stage, with the formation of an iron core which eventually collapses, leading to either a successful CCSN or a failed SN
$M_{\rm VMO}$ |  Mass boundary between massive stars and VMO

###  3.2 Calculation of the chemical ejecta

Since with the adopted version of the PARSEC code, models of massive and very
massive stars are evolved until central carbon exhaustion, we combine our
evolutionary tracks with extant explosive models covering a range of initial
masses that corresponds to different final fates (CCSN, failed SN, PPISN and
PISN). Following the work by Slemer (2016), for each stellar model of given
initial mass $M_{\rm i}$ we first compute the amount of ejected mass of the
element j due to the stellar winds, $E_{j}^{\rm w}(M_{\rm i})$, and then the
contributions of the associated supernova channels, $E_{j}^{\rm sn}(M_{\rm
i})$, as detailed in Sect. 3.5. The total ejecta $E_{j}(M_{\rm i})$ are given
by

$E_{j}(M_{\rm i})=E_{j}^{\rm w}(M_{\rm i})+E_{j}^{\rm sn}(M_{\rm i})$ (1)

The complete tables of wind and explosion ejecta for massive and very massive
stars ($8\leq\mbox{$M_{\rm i}$}/\mbox{${M}_{\odot}$}\leq 350$) and four values
of the initial metallicity ($0.0001\leq\mbox{$Z_{\rm i}$}\leq 0.02$) are
described in Appendix A. They are available online
666http://gofile.me/55RoU/qnITvlSQU.

We consider the most important chemical species and their isotopes from H to
Zn. The species explicitly included in PARSEC nuclear networks are all the
isotopes from $\mathrm{{}^{1}H}$ to $\mathrm{{}^{28}Si}$. Heavier elements are
present in the initial chemical composition, according to the adopted scaled-
solar mixture (see Sect. 3), and are not affected by the nuclear reactions and
mixing events during the hydrostatic H- and He-burning phases.

In the following we detail how we calculate the wind and explosion ejecta,
$E_{j}^{\rm w}(M_{\rm i})$ and $E_{j}^{\rm sn}(M_{\rm i})$, that appear in
Eq.(1).

###  3.3 Wind ejecta

The wind ejecta of a species $j$ contributed by a star of initial mass $M_{\rm
i}$ is computed with the following equation:

$E^{\rm w}_{j}(M_{\rm i})=\int_{0}^{\tau_{\rm C}}\dot{M}(M_{\rm
i},t)X_{j}^{\rm s}(t)dt\,$ (2)

where the integral is performed over the stellar lifetime, from the zero age
main sequence (ZAMS) up to the stage of carbon ignition, $\tau_{\rm C}$. For a
given $M_{\rm i}$ the quantities $\dot{M}(M_{\rm i},t)$ and $X_{j}^{\rm s}(t)$
denote, respectively, the mass-loss rate and surface abundance (in mass
fraction) of the species j, at the current time t.

The total amount of mass lost by a star during its hydrostatic evolution,
$\mbox{$M_{\rm i}$}-\mbox{${M}_{\rm fin}$}$, can be appreciated from Fig. 1,
which shows the pre-SN mass ($M_{\rm fin}$) as a function of $M_{\rm i}$, for
a few selected values of the initial metallicity.

Figure 12, in the appendix, illustrates the fractional wind ejecta, $E^{\rm
w}_{j}(M_{\rm i})/M_{\rm i}$, of the main chemical species considered in the
PARSEC models, as a function of $M_{\rm i}$ and $Z_{\rm i}$. For $M_{\rm
i}\leq 100\,\mbox{${M}_{\odot}$}$ the wind ejecta generally increase with
initial mass and metallicity, which is explained by the strengthening of
stellar winds at higher luminosities and larger abundances of metals. This
applies to H, He, N, Ne, Na, Mg, Al, and Si.

The metallicity trends for C and O reverse in the case of VMO. Compared to the
predictions for $\mbox{$Z_{\rm i}$}=0.02$, at low metallicities and high
initial masses, $\mbox{$Z_{\rm i}$}\leq 0.006$ and $\mbox{$M_{\rm
i}$}>100-200\,\mbox{${M}_{\odot}$}$, the wind ejecta may be larger by even one
order of magnitude for C and up to two orders of magnitude for O. This result
is explained when considering the stage at which stars of different $M_{\rm
i}$ and $Z_{\rm i}$ enter the WC and WO phases, which are characterized by
powerful winds enriched in C and O.

At $\mbox{$Z_{\rm i}$}=0.02$ all VMO experience high mass loss before entering
in the WC regime, which is attained close to the end of the He-burning phase.
We note that these models are not expected to go through the WO regime. As a
consequence, their ejecta are characterized by low amounts of primary C and O.
Conversely, at lower metallicities, $\mbox{$Z_{\rm i}$}=0.0001$ and $0.004$,
due to the relatively weak stellar winds during the early evolutionary stages,
VMO reach the WC and WO regimes with a much larger mass, hence producing
higher ejecta of C and O.

###  3.4 Explosion ejecta of electron capture supernovae

To account for the ECSN channel we take advantage of the recent revision of
the PARSEC code (Costa et al. 2020), who extended the sequence of hydrostatic
nuclear burnings up to oxygen. In this way we can check which models develop a
degenerate O-Ne-Mg core after the carbon burning phase. As mentioned in Sect.
4, we did not follow the Super-AGB phase and the corresponding yields are
taken from Ritter et al. (2018) using the models with $\mbox{$M_{\rm
i}$}=6\,\mbox{${M}_{\odot}$}$ and $7\,\mbox{${M}_{\odot}$}$.

As to the ECSN channel we proceed as follows. Given the severe uncertainties
that affect the definition of the mass range for the occurrence of ECSN (e.g.,
Doherty et al. 2017; Poelarends et al. 2008) and the modest chemical
contribution expected from the explosive nucleosynthesis (e.g., Wanajo et al.
2009), we adopt a simple approach. For each $Z_{\rm i}$ we look over the mass
range $8\leq\mbox{$M_{\rm i}$}/\mbox{${M}_{\odot}$}\leq 10$, and assign the
ECSN channel to the PARSEC models that, after the carbon-burning phase,
develops a degenerate core with mass close to the critical value of
$1.38\,\mbox{${M}_{\odot}$}$. In the metallicity range under consideration
($0.0001\leq\mbox{$Z_{\rm i}$}\leq 0.02$), this condition is met by PARSEC
models with $8\leq\mbox{$M_{\rm i}$}/\mbox{${M}_{\odot}$}\leq 9$.

The ECSN explosion ejecta are taken from the work of Wanajo et al. (2009,
their table 2) using the FP3 model as suggested by the authors. The
nucleosynthesis results derive from a neutrino-driven explosion of a
collapsing O–Ne–Mg core of mass $=1.38\,\mbox{${M}_{\odot}$}$, with a stellar
progenitor of $\mbox{$M_{\rm i}$}=8.8\,\mbox{${M}_{\odot}$}$ (Nomoto 1984).
According to this model, the total mass ejected by the explosion is quite low,
$\simeq 1.39\times 10^{-2}\,\mbox{${M}_{\odot}$}$. This fact, together with
the neutron-richness of the ejecta, lead to a modest production of radioactive
56Ni and hence of stable 56Fe ($\approx 0.002-0.004\,\mbox{${M}_{\odot}$}$).
In addition, the ECSN yields are characterized by a minor production of
-elements (e.g., O and Mg), and an appreciable synthesis of heavier species
like 64Zn and some light p-nuclei (e.g., 74Se, 78Kr, 84Sr, and 92Mo). Assuming
no fall back during the explosion, the ECSN event is expected to produce a
neutron star with mass $\mbox{${M}_{\rm rem}$}=1.366\,\mbox{${M}_{\odot}$}$.
Finally, we add the ejecta of the layers above ${M}_{\rm rem}$ with the
chemical composition predicted by the corresponding PARSEC model.

The specific ECSN model adopted for each $Z_{\rm i}$ serves as a bridge
between AGB and massive stars, to avoid a coarse mass-interpolation of the
ejecta in the transition region $\mbox{$M$}_{\rm AGB}<\mbox{$M_{\rm
i}$}/\mbox{${M}_{\odot}$}<M_{\rm mas}$ (see Table 1). We note that, with a
canonical IMF extending up to 350 ${M}_{\odot}$, the weight of the
$\mbox{$M$}_{\rm AGB}-M_{\rm mas}$ range should be not that large and only a
few elements would be affected, that are not the focus of the present paper. A
more careful consideration of this mass interval is postponed to a future
study.

###  3.5 Explosion ejecta of core collapse supernovae

Figure 2: ${M}_{\rm He}$(left panel) and ${M}_{\rm CO}$(right panel) as a
function of $M_{\rm i}$, for different values of $Z_{\rm i}$. Lines show the
data extracted from PARSEC stellar evolution models, while crosses represent
the models of Chieffi & Limongi (2004).

Stars in this class have $M_{\rm mas}<\mbox{$M_{\rm i}$}<\mbox{${M}_{\rm
VMO}$}$. The upper limit, ${M}_{\rm VMO}$, corresponds to a star that reaches
$M_{\rm He}\sim 32\mbox{${M}_{\odot}$}$ after central He burning and enters
the pair-instability regime during O-burning (Woosley 2017), thus avoiding the
standard evolutionary path to the silicon burning stage (see Sect. 3.6). At
solar composition, this mass limit is typically $\mbox{${M}_{\rm VMO}$}\sim
100\,\mbox{${M}_{\odot}$}$ but it is expected to vary with metallicity, as it
is affected by mass loss during the early evolutionary phases.

Of particular relevance for this mass range is the determination of the
explodability of a model, i.e. the conditions that lead to a successful SN or
to a failed SN. In recent years there have been many attempts to explore the
dependence of the outcome of the supernova collapse on the input physics, with
the final goal to possibly determine a relation between the explodability and
the main stellar parameters, in particular the initial or pre-supernova mass
of the star (Fryer 1999; Fryer et al. 2006; Fryer & Kalogera 2001; Heger &
Woosley 2002; Heger et al. 2003a; O’Connor & Ott 2011; Fryer et al. 2012;
Janka 2012; Ugliano et al. 2012; Ertl et al. 2016).

We briefly recall that Fryer et al. (2012) provided simple relations of the
explodability with the final C-O core mass (${M}_{\rm CO}$), O’Connor & Ott
(2011) introduced the compactness criterion as a threshold for the
explodability and Ertl et al. (2016) introduced a two-parameter explodability
criterion. Fryer et al’s models depend only on the C-O core mass after carbon
burning and on the pre-supernova mass of the star. Instead, the other models
predict a non-monotonic behaviour of the explodability with the core mass,
with the existence of islands of explodability intermixed with islands of
failures.

A related issue concerns the material that may fall back onto the surface of
the proto-neutron star after the explosion, leading eventually to the
formation of a black hole and possibly also to a failed SN (Fryer et al.
2012).

It is clear that the present theoretical scenario is heterogeneous and there
is no unanimous consensus of different authors on the explodability of a
massive star following the collapse of its Fe core. This is mainly due to the
fact that when the models are near the critical conditions for explosion they
become critically sensitive to slight variations in the input micro-physics
and numerical treatments (Burrows et al. 2018).

All these facts make it hard to unambiguously set a threshold mass between
successful and failed explosions. However, since indications exist that a
reasonable limit could be in the range
$25\,\mbox{${M}_{\odot}$}\lesssim\mbox{$M_{\rm i}$}\lesssim
30\,\mbox{${M}_{\odot}$}$ (e.g. O’Connor & Ott 2011; Sukhbold et al. 2016), we
assume that massive stars with $\mbox{$M_{\rm i}$}\geq
30\,\mbox{${M}_{\odot}$}$ will fail to explode.

The other important parameter needed to obtain the ejecta is the remnant mass,
${M}_{\rm rem}$. To derive this quantity one could use the observed relation
between the ejecta of $\mathrm{{}^{56}Ni}$ and the pre-supernova mass of CCSN
(Umeda & Nomoto 2008; Utrobin & Chugai 2009; Utrobin et al. 2010). We note
that also this relation is affected by some uncertainty, in particular on the
determination of the pre-supernova mass ${M}_{\rm fin}$. For this reason, we
prefer to use, as the calibrating value for the models, the estimated value of
the $\mathrm{{}^{56}Ni}$ mass ejected by SN1987A, $\mathrm{{}^{56}Ni}\sim
0.07\mbox{${M}_{\odot}$}$ (Nomoto et al. 2013; Prantzos et al. 2018).

Given the explodability criterion and the ejected mass of
$\mathrm{{}^{56}Ni}$, we adopt suitable explosion models to derive the
corresponding ejecta. For this purpose we use the CCSN models by Limongi &
Chieffi (2003) and Chieffi & Limongi (2004) (hereafter CL04) because they
tabulate the explosion isotopes as a function of the internal mass coordinate.

Each stellar model of the PARSEC grid is characterized by four known
parameters, namely: $M_{\rm i}$, $Z_{\rm i}$, ${M}_{\rm fin}$and ${M}_{\rm
CO}$. We use the mass of the C-O core, ${M}_{\rm CO}$, to match the PARSEC
models to CL04 ones for $\mbox{$Z_{\rm i}$}=0.0001,0.001,0.006,0.02$. These
are the only values of $Z_{\rm i}$ in common between CL04 and PARSEC. Once
identified the CL04 explosion model that corresponds to a given ${M}_{\rm
CO}$, it is straightforward to integrate from the external layers inward until
the desired ejecta of $\mathrm{{}^{56}Ni}$ is reached. The corresponding mass
coordinate of the inner layer provides the mass cut, ${M}_{\rm cut}$, and
hence the explosion ejecta.

This scheme needs to be made a bit more articulated to take into account that
for the same $M_{\rm i}$ the PARSEC and CL04 models do not predict exactly the
same ${M}_{\rm CO}$. Simple interpolations are therefore applied.

We proceed as follows. For each PARSEC model of metallicity $Z_{\rm i}$ we
identify in the corresponding CL04 grid the two explosive models that bracket
the mass of the core, $M_{\rm CO1}<\mbox{${M}_{\rm CO}$}<M_{\rm CO2}$, with
pre-explosive masses $M_{\rm fin1}$ and $M_{\rm fin2}$, respectively. Using
the $\mathrm{{}^{56}Ni}$ criterion we derive the corresponding mass cuts,
$M_{\rm cut}(M_{\rm CO1})$ and $M_{\rm cut}(M_{\rm CO2})$, and the explosive
ejecta, integrating from $M_{\rm cut}(M_{\rm CO1,2})$ to $M_{\rm fin1,2}$.
Finally, we use ${M}_{\rm CO}$ of the PARSEC model as interpolating variable
to obtain $M_{\rm cut}(\mbox{$M_{\rm i}$},\mbox{$Z_{\rm i}$},\mbox{${M}_{\rm
CO}$})$ and the explosion ejecta $E^{\rm sn}_{j}(\mbox{$M_{\rm
i}$},\mbox{$Z_{\rm i}$},\mbox{${M}_{\rm CO}$})$ for all chemical species under
consideration.

To estimate the mass of the remnant, ${M}_{\rm rem}$, we assume that in
successful CCSN the efficiency of fall-back is negligible, as shown by recent
hydrodynamical simulations (Ertl et al. 2016). It follows that
$\mbox{${M}_{\rm rem}$}=\mbox{${M}_{\rm cut}$}$ for successful CCSN and
$\mbox{${M}_{\rm rem}$}=\mbox{${M}_{\rm fin}$}$ for failed SN. As to the
nature of the compact remnant, we assign a neutron star for $\mbox{${M}_{\rm
rem}$}<2.9\,\mbox{${M}_{\odot}$}$, or a black hole otherwise (Tews & Schwenk
2020; Kalogera & Baym 1996).

Before closing this section, a few remarks are worth. The first applies to the
matching parameter ${M}_{\rm CO}$. Figure 2 compares the PARSEC values of
${M}_{\rm He}$ and ${M}_{\rm CO}$ with those derived from CL04 models, as a
function of $M_{\rm i}$. We note that the values of ${M}_{\rm He}$ and
${M}_{\rm CO}$ of our PARSEC models are slightly larger than predicted by
CL04. This is due to the fact that in PARSEC we adopt a slightly more
efficient core overshooting. An implication of this difference will be
discussed later (Sect. 5).

The second is that we assume that differences in the internal configurations
between PARSEC models and pre-explosive models with the same $M_{\rm CO}$, do
not impact the final nucleosynthesis outcome. Differences are in fact expected
because the two evolutionary models adopt different mass-loss rates, mixing
schemes and distribution of heavy elements.

The latter difference is considered when computing the yields of newly
produced elements by properly accounting for the initial composition. As to
the differences in input physics, we note the following. The stellar models of
the CL04 grid were computed at constant mass, while our PARSEC tracks include
mass loss by stellar winds for $\mbox{$M_{\rm i}$}\geq
14\,\mbox{${M}_{\odot}$}$. However, this difference should not affect our
results because, besides the fact that we match the models using ${M}_{\rm
CO}$(which somewhat alleviates the problem of the different mixing scheme),
mass loss is not so important for the progenitors of successful CCSN with
$\mbox{$M_{\rm i}$}\leq 25\,\mbox{${M}_{\odot}$}$, and especially for
$\mbox{$Z_{\rm i}$}\leq 0.006$. Powerful stellar winds affect the pre-
supernova evolution of more massive stars, but in this case the matching with
CCSN explosive models is not required, as these stars fail to explode and only
their wind ejecta are considered. We conclude this section noting that, for
example, differences between the new yields of massive stars obtained in this
work and those of Limongi & Chieffi (2018, LC18) for zero rotational velocity
are generally less than those between LC18 models with different rotational
velocities, as discussed in Sect. 5.

###  3.6 Explosion ejecta of pulsational pair instability and pair instability
supernovae

Very massive stars that develop a final helium core mass in the range between
$\sim 32\,\mbox{${M}_{\odot}$}$ and $\sim 64\,\mbox{${M}_{\odot}$}$ are
expected to enter the domain of pulsational pair-instability supernovae
(PPISN), before ending their life with a successful or failed core-collapse
supernova (Woosley et al. 2002; Chen et al. 2014; Yoshida et al. 2016; Woosley
2017). During the pair-instability phase, several strong pulses may eject a
significant fraction of the star’s residual envelope and, possibly, a small
fraction of the core mass. In contrast, the thermonuclear ignition of oxygen
in stars with helium core masses between $\mathrm{\sim
64\mbox{${M}_{\odot}$}}$ and $\sim 135\,\mbox{${M}_{\odot}$}$ leads to a pair-
instability supernova (PISN), assimilated to a single strong pulse that
disrupts the entire star, leaving no remnant behind (Heger & Woosley 2002;
Heger et al. 2003b).

PISN have been usually associated to the first, extremely metal-poor stellar
generations (e.g. Karlsson et al. 2008). However, recent stellar evolution
models suggest that PISN could occur also for stars with initial metallicity
$\mbox{$Z_{\rm i}$}\approx\mbox{${Z}_{\odot}$}/3$, which implies that they are
potentially observable even in the local universe (Yusof et al. 2013; Kozyreva
et al. 2014b). For these reasons, PPISN and PISN may play a key role to
understand the chemical evolution of the Galaxy.

While the physical mechanisms behind PPISN and PISN are quite well understood,
severe uncertainties affect the range of helium (or, equivalently, carbon-
oxygen) core masses that drive stars to enter the pair-instability regime
(e.g. Leung et al. 2019; Farmer et al. 2019; Marchant & Moriya 2020; Renzo et
al. 2020; Costa et al. 2020). In this work we adopt the indications from
Woosley (2017), who suggests $32\lesssim\mbox{${M}_{\rm
He}$}/\mbox{${M}_{\odot}$}\lesssim 64$ for PPISN and
$64\lesssim\mbox{${M}_{\rm He}$}/\mbox{${M}_{\odot}$}\lesssim 135$ for PISN.

We model PPISN as a super-wind phase that ejects the surface layers, without
any appreciable synthesis of new elements, until the star collapses to a BH.
For each PARSEC model with a given helium core mass, the corresponding
${M}_{\rm rem}$ is obtained by interpolation in ${M}_{\rm He}$, between the
values tabulated by Woosley (2017). Then, the total PPISN ejecta are estimated
by integrating in mass the PARSEC structures from ${M}_{\rm rem}$ to ${M}_{\rm
fin}$. Finally, we add the PARSEC wind ejecta.

For PISN, we calculate the explosion ejecta from the zero-metallicity pure-
helium stellar models provided by Heger & Woosley (2002), as a function of the
mass $M_{\mathrm{He-star}}$. Similarly to the case of PPISN, the helium core
mass, ${M}_{\rm He}$, of our PARSEC tracks is taken as interpolating variable
to perform the match with the explosion models and derive the ejecta. Finally,
we add the PARSEC contributions of all layers from ${M}_{\rm He}$ to ${M}_{\rm
fin}$ and the wind ejecta.

We use the models of Heger & Woosley (2002), computed at $\mbox{$Z_{\rm
i}$}=0$, to derive the PISN ejecta for VMO with $\mbox{$Z_{\rm i}$}>0$. This
assumption is reasonable since the total ejected mass of metals777According to
standard terminology, metals refer to chemical species heavier than helium.,
$M_{\mathrm{metals}}$, comprises most of the PISN ejecta, with a fractional
contribution, $M_{\mathrm{metals}}/M_{\mathrm{He-star}}$, that is generally
larger than $97\%$ for all tabulated models.

It is worth noticing that stars with $\mbox{$M_{\rm i}$}\geq\
90\,\mbox{${M}_{\odot}$}$ enter the PPISN and PISN regimes already at
$\mbox{$Z_{\rm i}$}=0.006$ (see Fig. 3). Our models agree with earlier
theoretical findings (e.g., Kozyreva et al. 2014b) and support the hypothesis
that some superluminous supernovae recently observed at metallicity
$\approx\mbox{${Z}_{\odot}$}/3$, may be explained through the pair-instability
mechanism, provided the IMF extends to VMOs (Woosley et al. 2007; Kozyreva et
al. 2014a). In this respect, we note that, according to IGIMF analysis by
Jeřábková et al. (2018), galaxies with metallicity $\mbox{\rm[{\rm Fe}/{\rm
H}]}<0$ dex and $\mathrm{SFR>1}$ $\mbox{${M}_{\odot}$}/{\rm yr}$ are
characterized by a top-heavy and bottom-light IGIMF, as compared to the
canonical one. Using the IGIMF grid provided by Jeřábková et al. (2018) and
assuming, for example, a metallicity [Fe/H] = -1 dex and a SFR $\sim$2
Msun/yr, the slope in the high mass tail would be x $\sim$ -2.1. Extrapolating
this slope up to 200 ${M}_{\odot}$ instead of 150 ${M}_{\odot}$, the value
adopted in the grid, this will allow about 66700 stars with 10
${M}_{\odot}$$\leq$ $M_{\rm i}$ $\leq$120 ${M}_{\odot}$and 2000 stars with 120
${M}_{\odot}$$\leq$ $M_{\rm i}$ $\leq$200 ${M}_{\odot}$, to be born every 3
Myr (the lifetime of the most massive stars).

Figure 3: Final fate of massive and very massive stars as a function of
$M_{\rm i}$ and $Z_{\rm i}$. Green dot is a successful SN, from ECSN if the
background is light green or CCSN if it is yellow; red dot is a BH from a
failed CCSN; red dot in a black box is a BH from PPISN; yellow star is a
thermonuclear explosion from PISN and black dot is a DBH.

###  3.7 Ejecta of very massive stars that directly collapse to black holes

If a star is massive enough to build a helium core with $\mbox{${M}_{\rm
He}$}>135\,\mbox{${M}_{\odot}$}$, no material will be able to avoid the direct
collapse into a black hole (DBH), induced by the pair creation instability.
Under these conditions no explosive ejecta are produced (Fryer & Kalogera
2001; Heger & Woosley 2002; Nomoto et al. 2013), and the only chemical
contribution comes from wind ejecta. With the adopted mass-loss rates in
PARSEC, these objects appear only at a low metallicity, $\mbox{$Z_{\rm
i}$}=0.0001$, and initial masses $\mbox{$M_{\rm
i}$}>200\,\mbox{${M}_{\odot}$}$. Conversely, at larger metallicities stars
with $\mbox{$M_{\rm i}$}>200\,\mbox{${M}_{\odot}$}$ avoid the DBH channel,
since mass loss is efficient enough to drive their He-core masses into the
regimes of PISN or failed CCSN.

##  4 Ejecta of AGB stars

We have complemented the ejecta of massive stars with those of AGB stars
computed with the COLIBRI code (Marigo et al. 2013), in the mass range
$0.7\lesssim\mbox{$M_{\rm i}$}/\mbox{${M}_{\odot}$}\lesssim 6$ and for the
same values of the PARSEC metallicities, already mentioned in Sect. 3.1. These
models follow the whole thermally pulsing phase, TP-AGB, up to the ejection of
the entire envelope by stellar winds. The initial conditions are taken from
the PARSEC grid of stellar models at the first thermal pulse or at an earlier
stage on the Early-AGB. COLIBRI and PARSEC share the same input physics (e.g.,
opacity, equation of state, nuclear reaction rates, mixing-length parameter)
and the numerical treatment to solve the structures of the atmosphere and the
convective envelope. For these reasons, the PARSEC+COLIBRI combination
provides a dense, homogeneous and complete grid of models for low- and
intermediate-mass stars (roughly $\simeq 70$ values of $M_{\rm i}$ for each
metallicity value).

In COLIBRI models, the parameters describing the main processes that affect
the TP-AGB phase, such as the mass-loss rates and the efficiency of the third
dredge-up, have been thoroughly calibrated with observations of AGB stars in
the Galaxy, Magellanic Clouds, and low-metallicity nearby galaxies (Girardi et
al. 2010; Rosenfield et al. 2014, 2016; Marigo et al. 2017; Lebzelter et al.
2018; Pastorelli et al. 2019, 2020; Marigo et al. 2020). The COLIBRI yields
account for the chemical changes due to the first, second, third dredge-up
episodes and hot-bottom burning in the most massive AGB stars ($\mbox{$M_{\rm
i}$}\gtrsim 3-4\,\mbox{${M}_{\odot}$}$), and include the same chemical species
as in PARSEC, from 1H to 28Si (see Sect. 3.1).

Finally, super-AGB stars are not treated explicitly here and their ejecta are
taken from Ritter et al. (2018), for stars with $\mbox{$M_{\rm
i}$}=6,7\,\mbox{${M}_{\odot}$}$ and $\mbox{$Z_{\rm
i}$}=0.0001,0.001,0.006,0.02$ 888https://github.com/NuGrid/NuPyCEE. The
chemical composition of the ejecta is the result of third dredge-up episodes
and hot-bottom burning. An overshoot scheme is applied to the borders of
convective regions, including the bottom of the pulse-driven convection zone.
As a consequence, the intershell composition is enriched with primary 16O
($\approx 15\%$) in Ritter et al. (2018) computations, much more than in
standard models without overshoot (${}^{16}\mathrm{O}\approx 1-2\%$), like in
Karakas (2010).

Table 2: Sets of chemical ejecta adopted in the chemical evolution models label | AGB stars | massive stars | rotation | PPISN/PISN/DBH | colour/line
---|---|---|---|---|---
MTW | M20 | TW | No | Yes | yellow
KTW | K10 | TW | No | No | blue (dashed)
Rr | R18 | R18r | No | No | green (continous)
Rd | R18 | R18d | No | No | green (dashed)
MLr | M20 | LC18 | Yes | No | cyan

##  5 Chemical ejecta from other authors

Here we present various combinations of chemical ejecta taken from the
literature and compare the main trends as a function of $M_{\rm i}$ and
$Z_{\rm i}$. They are summarized in Table 2. The different sets of ejecta are
then incorporated in our chemical evolution model of the Milky Way (see Sect.
6.2).

As to the yields from AGB stars, we consider three sets, namely: M20 (from the
COLIBRI code, Sect. 4), K10 (Karakas 2010) and R18 (Ritter et al. 2018). K10
provides the ejecta of AGB stars in the mass range $1\lesssim\mbox{$M_{\rm
i}$}/\mbox{${M}_{\odot}$}\lesssim 6$ for four metallicities ($\mbox{$Z_{\rm
i}$}=0.0001,0.004,0.008,0.02$). To obtain the yields at $\mbox{$Z_{\rm
i}$}=0.001$ and $\mbox{$Z_{\rm i}$}=0.006$, we interpolate in metallicity
between their original tables. R18 provide the ejecta of AGB and Super-AGB
stars in the mass range $1\leq\mbox{$M_{\rm i}$}/\mbox{${M}_{\odot}$}\leq 7$
for five metallicities ($\mbox{$Z_{\rm i}$}=0.0001,0.001,0.006,0.01,0.02$).

As to the yields of massive stars, including both wind and explosion
contributions, we consider three sets, namely: R18 (Ritter et al. 2018), L18
(Limongi & Chieffi 2018) with and without rotation, and TW that refers to the
new ejecta from this work (see Sect. 3). Nomoto et al. (2013) also published
yields for massive stars, which however do not include the wind contributions,
and therefore we did not consider them in our analysis.

R18 computed the ejecta of massive stars in the mass range
$12\leq\mbox{$M_{\rm i}$}/\mbox{${M}_{\odot}$}\leq 25$, for the same initial
metallicities as their AGB models, i.e. $\mbox{$Z_{\rm
i}$}=0.0001,0.001,0.006,0.01,0.02$, and for two models of explosion
conditions, rapid (Rr) or delayed (Rd), respectively (Fryer et al. 2006).

LC18 calculated the ejecta in the mass range $13\leq\mbox{$M_{\rm
i}$}/\mbox{${M}_{\odot}$}\leq 120$ for three different rotational velocities
($V_{\rm rot}=0,150,300$ km/s), and four metallicities ($\mbox{\rm[{\rm
Fe}/{\rm H}]}=0,-1,-2,-3$ dex). Here we use the version of their ejecta for
$\mbox{$Z_{\rm i}$}=0.0001,0.004,0.008,0.02$ publicly available on Github
NuPyCEE repository999https://github.com/NuGrid/NuPyCEE. Both R18 and LC18 sets
have the noticeable property that wind and explosion ejecta derive from
homogeneous stellar evolution models.

The distinguishing feature of our TW ejecta is that they range in mass beyond
the classical limit of $\mbox{$M_{\rm i}$}\simeq 100\,\mbox{${M}_{\odot}$}$,
extending up to $\mbox{$M_{\rm i}$}=350\,\mbox{${M}_{\odot}$}$, hence opening
the possibility to investigate the chemical role of VMO in terms of stellar
winds, PPISN and PISN explosions, and DBH channel.

To make a meaningful comparison among the different sets of ejecta we opt to
use a dimensionless quantity, defined as the ratio between the newly produced
yield of a given species j and the stellar initial mass, $P_{j}(\mbox{$M_{\rm
i}$})$:

$P_{j}(\mbox{$M_{\rm i}$})=\left[E_{j}(\mbox{$M_{\rm i}$})-(\mbox{$M_{\rm
i}$}-\mbox{${M}_{\rm rem}$})\,X_{j,0}\right]/\mbox{$M_{\rm i}$}$ (3)

where the total ejecta $E_{j}(\mbox{$M_{\rm i}$})$ is defined by Eq.(1), and
$X_{j,0}$ is the initial stellar abundance (in mass fraction) of the element
j.

The different sets are compared in Figs. 13 \- 17 as a function of $M_{\rm i}$
and a few values of the initial metallicity, $\mbox{$Z_{\rm
i}$}=0.0001,0.001,0.006,0.02$ respectively. We note that, for comparison
purposes only, the LC18 yields for $\mbox{$Z_{\rm i}$}=0.006$ are obtained
through a metallicity interpolation. The various panels show the ejecta of
4He, 12C, 14N, 16O, 20Ne, 24Mg, 28Si, S, Ar, Ca, Ti and Fe.

The AGB ejecta of 4He, 12C, 14N exhibit significant differences among
different sets. At $\mbox{$Z_{\rm i}$}=0.0001$ the K10 ejecta are much larger
than those of R18 and M20. In general, the production of 12C and 14N predicted
by K10 is much higher than M20, reaching a factor of ten for 14N at
$\mbox{$Z_{\rm i}$}=0.006$. At increasing metallicity the differences become
less pronounced. At $\mbox{$Z_{\rm i}$}=0.02$ the trend reverses with K10
predicting the lowest yields, but for the most massive AGB stars with hot-
bottom burning. These discrepancies are mainly the result of different input
physics (e.g., molecular opacities, mixing length parameter), mass-loss
prescriptions, as well as differences in the efficiency of the third dredge-
up.

In the domain of massive stars, our TW ejecta agree fairly well with non-
rotating LC18 at both $\mbox{$Z_{\rm i}$}=0.0001$ and $\mbox{$Z_{\rm
i}$}=0.006$, while the comparison slightly worsens at $\mbox{$Z_{\rm
i}$}=0.02$, likely because the effect of mass loss becomes important at higher
metallicity. We note that our TW set produces slightly larger fractions of
16O, 20Ne and 24Mg than non-rotating LC18 at any $Z_{\rm i}$.

At $\mbox{$Z_{\rm i}$}=0.0001$, the rotating LC18 models yield a much larger
fraction of 14N and, to a much less extent, 16O and 12C, compared to the non-
rotating set. This trend remains at increasing metallicity ($\mbox{$Z_{\rm
i}$}=0.006$ and $\mbox{$Z_{\rm i}$}=0.02$) but the over-production of 14N
appears less pronounced. Conversely, species such as 20Ne and 24Mg are
produced less by stars with rotation.

The comparison between TW and R18 shows that, at $\mbox{$Z_{\rm i}$}=0.0001$,
there is a fairly good agreement for 14N, 16O, 28Si, S, Ar and Ca. At the same
time, our TW ejecta produce less 12C, more 20Ne and 24Mg than R18. At
$\mbox{$Z_{\rm i}$}=0.006$ the TW predictions for 12C, which agree well with
non-rotating LC18, are about twice the R18 ejecta. We note that at this
metallicity R18 presents a notable Fe production, higher by more than a factor
of three compared to TW and LC18. A less pronounced, but still large Fe yield
is predicted by R18 also at lower metallicity ($\mbox{$Z_{\rm i}$}=0.0001$).
Possible consequences of such Fe over-production will be discussed later in
the sections devoted to chemical evolution models.

We note that none of the stellar models in the LC18 grid reaches ${M}_{\rm
He}$ high enough to enter the pair-instability regime. Recalling that stars
with $\mbox{$M_{\rm i}$}>30\,\mbox{${M}_{\odot}$}$ fail to explode and
collapse to a BH, it follows that in the mass range, $30\leq\mbox{$M_{\rm
i}$}/\mbox{${M}_{\odot}$}\leq 120$, LC18 ejecta are only due to stellar winds
and become null for $\mbox{$M_{\rm i}$}>120\,\mbox{${M}_{\odot}$}$. Since the
maximum mass in the R18 grid is $M_{\rm i}$2̄5 ${M}_{\odot}$, beyond this
limit all R18 ejecta are zero. It follows that we can analyse the total ejecta
of VMO, $\mbox{$M_{\rm i}$}\gtrsim 100\,\mbox{${M}_{\odot}$}$, by only
referring to our MTW set.

At $\mbox{$Z_{\rm i}$}=0.02$, the most important contributions from VMO
correspond to 4He, 12C and 14N. For all the other elements there is no
significant production because, due to the large mass-loss rates at this
metallicity, our models do not enter the PISN channel. The large production of
14N is due to the large convective cores of the high mass stars, that mix the
CNO products into the external regions, subsequently exposed to the ISM by
mass-loss. We note that, for $M_{\rm i}$$\leq 30\ \mbox{${M}_{\odot}$}$, the
14N production in our models is only slightly larger than that of LC18 models
with zero rotational velocity and so significantly less than that of their
models with high rotational velocities. However at $M_{\rm i}$$\gtrsim
30\mbox{${M}_{\odot}$}$, our models have a 14N wind production that is almost
identical to that of LC18 models with high rotational velocities. Since the
latter is about 60 % larger than that of their zero rotation models at the
same $M_{\rm i}$ , we identify the reason for this difference in the larger
mass-loss rates and consequent more rapid pealing and faster core mass
decrease in the LC18 models. At $\mbox{$Z_{\rm i}$}=0.006$ VMO are expected to
eject appreciable amounts of newly produced 4He and 12C, while the yield of
14N decreases considerably. At the same time, other species provide notable
contributions, such as 16O, 24Mg, 28Si, S, Ar, Ca, Ti and Fe. The yields of
these nuclides increase at higher $M_{\rm i}$. This is a clear effect of the
occurrence of PPISN and PISN which is favoured at lower metallicities. At
$\mbox{$Z_{\rm i}$}=0.001$ and $\mbox{$Z_{\rm i}$}=0.0001$ the ejecta of all
models with $\mbox{$M_{\rm i}$}>100\,\mbox{${M}_{\odot}$}$ have typical
signatures of these explosive events. We note that at these metallicities the
yields of 56Fe may reach extremely high values, up to
$20\,\mbox{${M}_{\odot}$}$ and $40\,\mbox{${M}_{\odot}$}$, respectively.

,

,

Figure 4: From top-left to bottom-right: evolution of SFR, gas fraction, SNII
and SNIa rates, and total gas metallicity. Different lines correspond to the
best chemical evolution models obtained with different sets of yields, as
indicated. Vertical bars at an age of 13 Gyr show their present-day estimated
values. In the bottom-right panel the solid triangle at the age of 8.4 Gyr
marks the protosolar metallicity, $Z_{\rm protoSUN}=0.017$, resulting from the
PARSEC calibration. The parameters of the chemical evolution models are listed
in Table 3.

Figure 5: Observed [Fe/H] distributions of thin-disc stars (solid brown) and
thick-disc stars (dotted yellow). Predicted thin-dick MDF, obtained with the
MTW set (solid blue), are superimposed. All histograms are self-normalized.
Vertical lines mark the median values of the distributions.

##  6 Chemical evolution of the Milky Way: thin and thick discs

We aim to investigate how the adoption of different chemical ejecta affect our
interpretation of the observed abundances in Milky Way (MW) thin- and thick-
disc stars.

###  6.1 The observed abundance data

The large amount of data collected over the years for stars in the solar
vicinity led to the definition of different Galactic components, namely: the
thin and thick discs, the halo and the -enhanced metal-rich population
(Allende Prieto et al. 2008; Gilmore et al. 2012; Zucker et al. 2012; de
Laverny et al. 2013). The populations of the MW disc can be distinguished in
various ways. For example, by adopting kinematic parameters, Jurić et al.
(2008) pointed out that the stellar number density distribution of the MW
could be well reproduced with two components with different scale heights
above the Galactic plane: the thick disc with a scale height of
$\mathrm{\simeq 900\,pc}$, and the thin disc with a scale height
$\mathrm{\simeq 300\,pc}$.

At the same time, chemical abundances reveal the existence of clearly separate
sequences of -elements as a function of [Fe/H], with thick-disc stars
generally belonging to a high [/Fe] (-enhanced) sequence, while thin-disc
stars exhibiting a lower [/Fe] ratio at the same [Fe/H] (e.g. Plevne et al.
2020; Grisoni et al. 2017; Kawata & Chiappini 2016; Bekki & Tsujimoto 2011;
Feltzing et al. 2003; Prochaska et al. 2000).

The classification based on kinematical properties and the one based on
abundance measurements provide somewhat different results so that it is not
clear which is the best way to group these stars (Boeche et al. 2013). In this
respect, it has often been pointed out that chemical evolution leaves a
persistent imprint that hardly changes, while kinematic properties are more
likely to vary as they may be affected by dynamical interactions (Schönrich &
Aumer 2017; Vera-Ciro et al. 2016).

In recent years, ages of individual stars have been measured with sufficient
accuracy to be used as robust population indicators, like in the case of star
clusters (Fuhrmann 2011). We emphasize that age cannot be taken as a proxy for
metallicity, rather it is a complementary independent parameter that concurs
to define the full population box, i.e. the distribution of stars in age and
abundances which, together with spatial and kinematic parameters, gives the
information necessary to reconstruct the star formation history in a galaxy.

In this work, we use the homogeneous set of data of disc stars provided by
Bensby et al. (2014), who conducted a high-resolution and high signal-to-noise
spectroscopic analysis of 714 F and G dwarf and subgiant stars in the solar
neighbourhood. This study is particularly suited for our purpose because,
based on the analysis of the their kinematical properties by Casagrande et al.
(2011), each star in the sample is assigned a relative membership probability,
$TD/D$ defined as the ratio between the thick-disc and thin-disc probability,
and $TD/H$ defined as the ratio between the thick-disc and halo probability.

Bensby et al. (2014) classified stars with $TD/D>2$ – having the probability
of belonging to the thick disc of least twice that of belonging to the thin
disc –, as potential thick-disc stars, while those with $TD/D<0.5$ as
potential thin-disc stars. Some of the thick disc stars where then assigned to
the halo population following the $TD/H<0.5$ kinematical criterion.

Adopting the same kinematical criteria, we count 387 thin disc stars, 203
thick disc stars and 36 halo stars. We discard 88 stars with $0.5<TD/D<2$. The
existence of at least two distinct disc sequences is clearly visible in the
abundance patterns that define the so-called -enhancement, as illustrated in
the [O/Fe], [Mg/Fe], [Si/Fe] and [Ca/Fe] vs. [Fe/H] diagrams of Fig. 6.
Interestingly, the [O/Fe] vs. [Fe/H] diagram clearly shows not only that the
two disc populations draw separate sequences, but also that the slopes of two
branches are different at increasing [Fe/H]. Interpreting the different slopes
with chemical evolution models is quite challenging, as already noted by
Kubryk et al. (2015).

###  6.2 Chemical evolution models

To analyze the evolution of the thin and thick disc populations we use the
chemical evolution model CHE-EVO (Silva et al. 1998), a one-zone open chemical
evolution code, that follows the time evolution of the gas abundances of the
elements, including infall of primordial gas. It has been used in several
contexts to provide the input star formation and metallicity histories to
interpret the spectro-photometric evolution of both normal and starburst
galaxies (e.g. Vega et al. 2008; Silva et al. 2011; Fontanot et al. 2009; Lo
Faro et al. 2013, 2015; Hunt et al. 2019). The equation describing the
evolution of the gas masses reads:

$\dot{M}_{{\rm g},j}=\dot{M}_{{\rm g},j}^{\rm SF}+\dot{M}_{{\rm g},j}^{\rm
FB}+\dot{M}_{{\rm g},j}^{\rm Inf}$ (4)

where, for the specie j, $\dot{M}_{{\rm g},j}^{\rm SF}$ represents the rate of
gas consumption by star formation, $\dot{M}_{{\rm g},j}^{\rm FB}$ is the rate
of gas return to ISM by dying stars and $\dot{M}_{{\rm g},j}^{\rm Inf}$ refers
to the infall rate of pristine material. For the star formation rate (SFR) we
adopt a modified Kennicutt-Schmidt law (Schmidt 1959; Kennicutt 1998):

$\psi(t)=\nu\,M_{\rm g}(t)^{k}$ (5)

where $\nu$ is the efficiency of star formation, $M_{\rm g}$ is the mass of
the gas and $k$ is the exponent of the star formation law. The quantity
$\dot{M}_{{\rm g},j}^{\rm FB}$ is calculated from the yields tables described
in Sect. 5, integrating the contributions of dying stars at any time-step. The
gas infall is assumed exponential, with an e-folding time scale $\tau_{\rm
inf}$ and a chemical composition equal to the primordial one (e.g. Grisoni et
al. 2017).

Type Ia supernovae are also taken into account according to the single
degenerate scenario and computed following the standard formalism first
introduced by Matteucci & Greggio (1986). The contribution of these sources to
the chemical enrichment is regulated by the parameter $A_{\rm SNIa}$ which
sets the fraction of the number of binary systems with total mass in the 3
${M}_{\odot}$-16 ${M}_{\odot}$range, effectively contributing to the SNIa
rate. Our adopted ejecta for SNIa are taken from Iwamoto et al. (1999). As to
AGB, massive, and very massive stars, we adopt the ejecta described in Sects.
3-5, and reported in Table 2. In addition to the evolution of elemental gas
masses, metallicity and gas fraction, the code provides also the evolution of
SNII and SNIa rates, and total mass in stars ($M_{\rm s}$).

A fundamental assumption of chemical evolution models is the IMF:

$\phi(\mbox{$M_{\rm i}$})=\frac{dn}{d\log(\mbox{$M_{\rm
i}$})}\propto\mbox{$M_{\rm i}$}^{-x}.$ (6)

We use a Kroupa-like three-slope power law IMF with $x=0.3$ for
$0.1\leq\mbox{$M_{\rm i}$}/\mbox{${M}_{\odot}$}\leq\,0.5$, $x=1.3$ for
$0.5\leq\mbox{$M_{\rm i}$}/\mbox{${M}_{\odot}$}\leq\,1$, while we vary the
slope for $\mbox{$M_{\rm i}$}>1\,\mbox{${M}_{\odot}$}$, as well as the upper
mass limit of the IMF, ${M_{\rm UP}}$, to search for best fitting models with
different ejecta combinations. We consider IMF slopes between $x=1.7$ (Kroupa
et al. 1993) and $x=1.3$ (Kroupa 2001; Chabrier 2003).

###  6.3 Previous Analyses of the MW thin and thick discs

There have been in the past and in the more recent literature many attempts to
explain the different chemical evolutionary paths of different MW components,
in particular of the thin and thick discs. The outcome of these studies is
that the observed different chemical evolutionary paths are related to
differences in the main physical processes that drive galaxy evolution, among
which the most significant are the gas accretion time-scale and the star
formation efficiency and, possibly, radial migration (Larson 1972; Lynden-Bell
1975; Pagel & Edmunds 1981; Matteucci & Greggio 1986; Matteucci & Brocato
1990; Ferrini et al. 1994; Prantzos & Aubert 1995; Chiappini et al. 1997;
Portinari & Chiosi 1999; Chiappini et al. 2001; Bekki & Tsujimoto 2011; Micali
et al. 2013; Sahijpal 2014; Snaith et al. 2014; Grisoni et al. 2017; Grand et
al. 2018; Spitoni et al. 2021). Good agreement between observations and
theoretical predictions for the Galaxy is obtained by models assuming that the
disc formed by infall of gas (Chiosi 1980; Matteucci & Francois 1989;
Chiappini et al. 1997). The formation of the different components is
associated to distinct sequential main episodes of gas accretion (infall
phases) that first rapidly accumulates in the central regions and then, more
slowly, in the more external ones, according to the so-called inside-out
scenario (Chiappini et al. 2001). In particular, the three-infall model,
devised by Micali et al. (2013), is able to reproduce the abundance patterns
of the MW halo, thick and thin disc at once. In this model, the halo forms in
a first gas infall episode of short timescale (0.2 Gyr) and mild star
formation efficiency, $\nu$=2 Gyr-1, lasting for about 0.4 Gyr. It is
immediately followed by the thick disc formation, characterized by a somewhat
longer infall timescale (1.2 Gyr), a larger duration of about 2 Gyr and a
higher star formation efficiency, $\nu$=10 Gyr-1. Finally, the star formation
continues in the thin disc with a longer infall timescale (6 Gyr in the solar
vicinity) and is still continuing nowadays, with star formation efficiency
$\nu$=1 Gyr-1. The [O/Fe] vs. [Fe/H] path is thus continuous across the
regions populated by halo, thick and thin disc stars. While in Micali et al.
(2013) the chemical enrichment is continuous across the three different infall
stages, Grisoni et al. (2017) used also an alternative scheme where the thin
and thick disc components evolve separately, in a parallel approach (see also
Chiappini (2009)). In the parallel approach, the disc populations are assumed
to form in parallel but to proceed at different rates. The gas infall
exponentially decreases with a timescale that is 0.1 Gyr and 7 Gyr, for the
thick and thin disc, respectively. This alternative approach better reproduces
the presence of the metal-rich $\alpha$-enhanced stars in the [Mg/Fe] vs.
[Fe/H] diagram obtained with the recent AMBRE data (Mikolaitis et al. 2017).

In our analysis, we will assume that the thin and thick disc populations
evolve separately as in the parallel model approach adopted by Grisoni et al.
(2017). This is clearly an oversimplification, because these stellar
components occupy the same volume in the solar neighbourhood, but nevertheless
it will allow at least to check our models against individually well-separated
populations. Alternatively, we could have considered the two populations
together and tried to obtain a model that recovers a sort of average path, as
done several times in the past. However, the evidence that the two populations
are different is so strong that reproducing their averaged properties is even
less meaningful.

##  7 Analysis of stellar abundances in the solar vicinity

###  7.1 Model constraints for the thin disc

For each combination of ejecta in Table 2, we build a large library of
chemical evolution models by varying the following parameters: the star
formation rate efficiency ($0.2\leq\nu\leq 2.0$) the normalization factor for
the SNIa ($0.02\leq A_{\rm Ia}\leq 0.1$), and the infall timescale
($0.1\leq\tau_{\rm inf}\leq 10$). For simplicity we set $k=1$ for all models.
With this choice for parameter ranges, the chemical evolution models are able
to bracket a few basic observational constraints for the thin disc, namely:

  * •

The current SFR of the MW. To a large degree it corresponds to that the thin
disc, and is estimated as $\rm{SFR}=0.65-3.0\mbox{${M}_{\odot}$}\,{\rm
yr}^{-1}$ (Robitaille & Whitney 2010).

  * •

The current gas fraction. It is assumed to be $M_{\rm g}/(M_{\rm g}+M_{\rm
s})\sim 0.2$ (Kubryk et al. 2015).

  * •

The current SNII rate. We set $R_{\rm SNII}=2\pm 1$ SNII events per century
(Prantzos et al. 2011).

  * •

The current SNIa rate. We set $R_{\rm SNIa}=0.4\pm 0.2$ SNIa events per
century (Prantzos et al. 2018).

  * •

The protosolar metallicity $Z_{\rm protoSUN}$. Though there are in literature
contrasting opinions concerning the significance of the differences between
Sun’s abundances and those of solar twin stars (e.g. Bensby et al. 2014;
Botelho et al. 2020), we take the protosolar metallicity as representative of
that of other disc stars with similar [Fe/H], and thus as a constraint for
disc chemical evolution models. The protosolar metallicity represents the bulk
metallicity of the molecular cloud out of which the Sun was born. It does not
coincide with the current photospheric solar metallicity, $Z_{\odot}\simeq
0.0134-0.0152$, from Asplund et al. (2009) and Caffau et al. (2011)
respectively, as a result of chemical sedimentation effects over a time of
about 4.6 Gyr (the present Sun’s age). Here we adopt $Z_{\rm protoSUN}=0.017$,
which is the initial metallicity of the PARSEC $1\,\mbox{${M}_{\odot}$}$ model
that best reproduces the currently observed Sun’s properties when using the
Caffau et al. (2011) solar mixture (Bressan et al. 2012). Assuming an age of
$13\,\mathrm{Gyr}$ for the formation of the Galaxy (e.g. Savino et al. 2020),
it follows that the Galactic age at the birth of the proto-Sun is $t_{\rm
protoSUN}=8.4\,\mathrm{Gyr}$. At this epoch, the metallicity in the solar
vicinity is $Z_{\rm protoSUN}$.

In addition to the constraints just mentioned above, we also require the
models to reproduce the observed ([Fe/H]) metallicity distribution function
(MDF) of thin-disc stars derived from the Bensby et al. (2014) data.

In total, we build a database of about 1300 models for each combination of
yields. Then, after identifying the group of models that satisfy all the
constraints within the uncertainty bars, for each set of yields we find the
best model that fits the observed path in the [O/Fe] vs. [Fe/H] plane (Bensby
et al. 2014), with the aid of a $\chi$-square minimization technique. Finally,
we perform a fine tuning and group the best models so as to minimize the
change of the parameters obtained with different yield sets, without worsening
the fits. This allows us to better isolate and analyze the effects produced by
changing either the yields or the IMF.

The evolution of the SFR, the gas mass fraction, the SNII and SNIa rates and
the gas metallicity ($Z$) of the selected models are shown in Fig. 4. The
predicted thin-disc MDF of our reference model, MTW, is compared with the
observed one in Fig. 5. In the same figure we also plot the observed MDF of
the thick-disc MDF, for sake of comparison. The vertical bars in the figure
mark the location of the median values of the distributions. Let us now
analyse the results obtained with the various yield combinations. The adopted
parameters of the selected models are summarized in Table 3.

  * •

MTW yields. The MTW model, that adopt the yields determined in this work, is
our reference model. It is shown in dark yellow in Fig. 4. This model
reproduces all the constraints fairly well, but for the gas fraction that
approaches the observed lower limit. The predicted MDF of the MTW model is
very similar to the observed one, with a median value in very close agreement
to that of the observed one. The thin disc MDF is significantly different from
that of the thick disc. The IMF used to reproduce the observed constraints has
a slope of $x=1.5$ in the high-mass range (Kroupa et al. 1993), and
$\mbox{${M_{\rm UP}}$}=120\,\mbox{${M}_{\odot}$}$.

  * •

KTW yields. The KTW model is shown in blue dashed in Fig. 4. The IMF used in
this case is the same as in the MTW model. Likewise, this model reproduces all
the constraints satisfactorily, including the the observed MDF, which is not
shown here for sake of clarity. We note that the metallicity of the KTW model
increases with time a bit faster than the MTW model. At ages of 2 Gyr and 8
Gyr the metallicity $Z$ in the KTW model is about 4% and 6% higher,
respectively, than predicted with MTW. This small difference is due mainly to
the larger AGB yields of 12C and 14N in K10 than in M20.

  * •

Rr and Rd yields. The Rr and Rd models are shown in green in Fig. 4. Both
provide similarly good fits as the other sets, but the model parameters are
quite different from the other cases, as shown in Table 3. First, the mass
limit is set to $\mbox{${M_{\rm UP}}$}=25\,\mbox{${M}_{\odot}$}$, in agreement
with their yield tables. Second, since the Fe production by CCSN in some of
these models is significantly higher than predicted by other authors, in order
to reproduce the observed [Fe/H] distribution and $Z_{\rm protoSUN}$ at the
Galaxy age, we need to decrease the SNIa efficiency factor to ASNIa=0.025 and
to adopt a flatter IMF in the high-mass regime.

  * •

MLr yields. The MLr model is shown in cyan in Fig. 4. Also this model
reproduces well the above observational constraints without the need to change
the chemical evolution parameters.

In general, we conclude that with all yield combinations it is possible to
find chemical evolution models able to nicely reproduce the observed
constraints for thin-disc component in the solar vicinity. As to the
metallicity distribution traced by [Fe/H], models generally tend to slightly
under-populate the central bins at lower metallicities, at the same time
extending a bit beyond the maximum observed [Fe/H]. Being of opposite sign,
the two deviations roughly compensate each other, allowing a good reproduction
of the observed median value. In principle, the agreement could be improved
even more through a fine tuning of the chemical evolution parameters. However,
we note the discrepancies between the predicted and observed thin-disc MDFs
are much lower than the intrinsic differences between the thin- and thick-disc
observed distributions. Therefore, we may reasonably consider our results as
fair models for the thin disc and proceed with the analysis of the predicted
abundance ratios.

Table 3: Parameters of the best chemical evolution models for each ejecta
combination, described in Table 2.

Ejecta | Chemical parameter | IMF
---|---|---
SET | ${\nu}$ | ${k}$ | $\tau_{\rm inf}$ | $A_{\rm SNIa}$ | ${M_{\rm UP}}$ | ${x}$
MTW | 0.8 | 1.0 | 6.0 | 0.04 | 120 | 1.5
KTW | 0.8 | 1.0 | 6.0 | 0.04 | 120 | 1.5
Rr | 0.4 | 1.0 | 3.0 | 0.025 | 25 | 1.3
Rd | 0.4 | 1.0 | 3.0 | 0.025 | 25 | 1.3
MLr | 0.8 | 1.0 | 6.0 | 0.04 | 120 | 1.5

  * *

The efficiency of star formation, , is in units of Gyr-1,
${\tau_{\rm inf}}$ is expressed in units of Gyr, and ${M_{\rm UP}}$ is given
in ${M}_{\odot}$.

,

,

Figure 6: Comparison of predicted abundances, derived with our best chemical
evolution models (see Fig. 4) with observations of thin-disc stars (blue
points) and thick-disc stars (magenta points) and halo stars (cyan triangles)
in the solar neighbourhood from Bensby et al. (2014). Plotted in the figure
are also the sample of metal-poor (black dots) and very metal-poor (gray dots)
stars from Nissen et al. (2002) and (Cayrel et al. 2004), respectively. Note
that the MTW and KTW models overlap almost completely.

###  7.2 Chemical evolution of the thin disc

For each set of chemical yields, we compare the elemental abundances predicted
by the corresponding best model with observations of thin-disc, thick-disc and
halo low-metallicity stars (Fig. 6). The four diagrams show the trends of a
few selected species, O, Mg, Si and Ca, as a function of [Fe/H]. We recall
that all the data are normalized according to the scaled-solar abundances
given by Caffau et al. (2011), which sets the reference solar mixture in
PARSEC.

  * •

MTW yields. This model is able to reproduce fairly well the observed trend in
the [O/Fe] diagram, with a modest over-production at the highest [Fe/H]. Such
small discrepancy should be attributed to the chemical yields from massive
stars. In fact, in the COLIBRI models used here (M20 yields), TP-AGB stars
produce a negligible amount of primary oxygen, as the chemical composition of
the intershell is the standard one and it contains no more than $1\%-2\%$ of
16O (see, e.g. Herwig 2000, models without overshooting).

Taking the MTW model as representative of the main thin-disc branch, we see
that at decreasing [Fe/H] the difference in [O/Fe] between thin and thick-disc
stars increases. At $\mathrm{\mbox{\rm[{\rm Fe}/{\rm H}]}\leq-1}$ the thin-
disc branch disappears, and the MTW model approaches the lower boundary of the
halo population.

Following the probability $TD/D$-criterion (Bensby et al. 2014), we note that
a second, less populated, thin-disc branch overlaps with the sequence of
thick-disc stars in the same abundance diagram. Clearly, The MTW model is not
able to reproduce this secondary branch.

Concerning the [Mg/Fe] ratio, the MTW model runs through the lower border of
the data, showing a well-known difficulty likely related with the 24Mg yields
(Romano et al. 2010). The [Si/Fe] ratio is fairly well recovered by this
model, while the [Ca/Fe] is under-produced.

  * •

KTW yields. This model behaves similarly to the MTW model for all the four
abundance ratios. We recall that KTW and MTW models only differ for the AGB
yields in use. Small differences appear at the lowest [Fe/H] values where the
increase of KTW metallicity with time takes place somewhat faster than
predicted by the MTW model (see Fig. 4 and the discussion in Sect 7.1). Since
MTW and KTW models share the same chemical evolution parameters and massive
star yields, differences in the trends of chemical species should be likely
ascribed to reaching somewhat different metallicities at the same evolutionary
time.

  * •

Rr and Rd yields. While successfully reproducing basic constraints of the MW
thin disc, these models fail to recover the evolution of the selected
abundance ratios (see Fig. 6). Even adopting a low SNIa efficiency parameter,
models exhibit a substantial deficit in 16O, 24Mg, 28Si and Ca relative to Fe.
In the attempt to solve the discrepancy we explored a wide range of chemical
evolution parameters, but we were unable to find better models than those
shown in Figs. 4 and 6. The results worsen for the rapid case and, since the
ratios [/Fe] run much flatter than observed, we argue that the issue may be
linked to the iron yields of CCSN. In fact, we find that at $\mbox{$Z_{\rm
i}$}<0.02$, iron production by CCSNs predicted by R18 is significantly higher
than in L18 explosive models (used in MTW, KTW and MLr), by a factor from two
to four.

  * •

MLr yields. The MLr model reproduces fairly well the [O/Fe] data of the thin
disc. The abundance of Mg is clearly under-predicted as already found by
Prantzos et al. (2018), while predictions for Si and Ca are able to populate
the regions of both thin- and thick-disc components. We note the MLr model
shows a general tendency to produce abundance ratios running with slightly
steeper slopes than observed.

From all the tests carried out with different combinations of chemical yields
we can draw a few conclusions. First, the chemical species discussed here are
marginally dependent on the chemical yields of AGB stars, and therefore cannot
be considered useful diagnostics for testing the goodness of low- and
intermediate-mass evolutionary models. This is not surprising since we do not
expect that AGB stars synthesize Fe and Ca, while they may be contributors of
Mg isotopes, which are present both in the dredged material and involved in
the Mg-Al cycle when hot-bottom burning is active in AGB stars with $M_{\rm
i}$¿ 3-4 ${M}_{\odot}$(Slemer et al. 2017; Marigo et al. 2013; Ventura &
D’Antona 2009). As discussed above, the production of some primary oxygen by
AGB stars depends on the inclusion of convective overshoot at the boundaries
of the pulse-driven convective zone (Herwig 2000), which applies to R18, but
not to M20 and K10 yields. In the context of this work, the role of AGB stars
as oxygen producers is not critical irrespective of the selected yield set.
This reinforces the conclusion that we need to consider other more suitable
elements, such as carbon and nitrogen, to compare and check different sets of
AGB yields. An in-depth analysis of AGB yields is postponed to a dedicated
future work.

It follows that the abundance trends investigated in this work are critically
dependent on the chemical yields from massive stars and VMO. Therefore, the
reader should keep in mind that, even when not explicitly stated, the
discussion that follows mainly deals with the effects produced by chemical
yields of stars with $\mathrm{M_{i}>8\mbox{${M}_{\odot}$}}$.

Once the chemical evolution models are calibrated on a few basic observables
of thin-disc stars, it is possible to reproduce fairly well the enrichment
paths of [O/Fe] and [Si/Fe] with most of the yield sets. As to the [Mg/Fe]
ratio, we meet the long-lasting problem of underproduction (e.g. Prantzos et
al. 2018), which appears somewhat less pronounced with the MTW yields. As to
the [Ca/Fe] the observations are better reproduced by massive star models with
rotation. Using the R18 yields, none of the elemental ratios is well
reproduced.The discrepancy is likely due to the high Fe yields in some of the
Ritter et al. (2018) explosion models, which makes it hard to recover the
observed [/Fe] ratios at a given [Fe/H].

Figure 7: Observed thick-disc MDF (Bensby et al. 2014, yellow thick solid
line) compared with the predictions of model TD1 (Table 4). The thin-disc MDF
is also plotted for comparison. Vertical lines are the median values of the
corresponding distributions.

###  7.3 Chemical evolution of the thick disc

We turn now our attention to the chemical evolution of thick-disc stars. The
diversity between thin (blue) and thick (magenta) disc stars is clearly
visible in the [O/Fe] vs. [Fe/H] diagram of Fig. 6. The models adopted for the
thin-disc stars are clearly not able to account for the [O/Fe] enrichment
history of thick-disc stars. Furthermore, stellar age determinations indicate
that thick-disc stars are, on average, older than thin-disc stars.

The origin of this marked difference may be linked to different star formation
histories, similarly to other components of the MW, such as the Galactic
Bulge, whose stars are even more metal-rich than thick-disc stars (Bensby et
al. 2017). Some chemo-dynamical models suggest that thick-disc stars may have
originated in the inner regions of the MW, characterized by an earlier and
faster enrichment so that they had the time to migrate into the solar vicinity
and become kinematically hotter. Under this hypothesis, these stars would be
naturally characterized by chemical compositions and, in particular, levels of
-enhancement different from the native stellar populations of the solar
neighbourhood (e.g. Aumer et al. 2017; Schönrich & Aumer 2017).

More recently the analysis of Gaia data (Gaia Collaboration et al. 2018),
provided evidence that the peculiar chemical composition of the thick disc
could have originated in an early merger of the MW with a satellite galaxy,
the Gaia Enceladus Sausage galaxy (Belokurov et al. 2018; Helmi et al. 2018;
Haywood et al. 2018).

Contrary to the thin disc case, where we can constrain the chemical evolution
models on a variety of observational data such as the current SFR, gas
fraction etc., for the thick disc the most robust constraint is the MDF. The
model that best reproduces the observed thick disc MDF, keeping the IMF
parameters identical to our model MTW for the thin disc, ${M_{\rm UP}}$ = 120
and $x$ = 1.5 and the same chemical evolution parameters $k$= 1 and
ASNIa=0.04, is obtained using $\tau_{inf}$=0.5 Gyr and $\nu$=1.4Gyr-1. The
resulting theoretical thick disc MDF (blue histogram) is compared with the
observed one (yellow histogram) in Fig. 7. The parameters of this model, named
TD1 are listed in Table 4. In order to fit the high [Fe/H] tail of the
observed MDF we need to include a galactic wind at an epoch of tGW=2.5 Gyr, to
expel the residual gas (about 10% of the total mass) and stop the SF. While
this assumption is introduced to cope with the simplicity of our chemical
evolution code, we note that the adopted tGW corresponds to the epoch of the
end of the second infall episode in the Micali et al. (2013) model. We note
that models with shorter $\tau_{inf}$ show a clear star number excess below
[Fe/H] $\sim$ -1, with respect to the observed one while, models with a larger
$\tau_{inf}$ show an excess in the high metallicity tail of the MDF. The
[O/Fe] vs. [Fe/H] abundance patterns obtained with model TD1, are shown in
Figure 8. This model is able to fit the observed region occupied by thick disc
stars but the slope is different from the observed one. In the [O/Fe] vs.
[Fe/H] diagram, it underestimates the high [O/Fe] values at low [Fe/H] and, it
overestimates the low [O/Fe] values at high [Fe/H]. This model is, however,
able to reproduce the region occupied by the most metal-poor stars, likely
belonging to the halo population, that will be discussed later.

Figure 8: Comparison of the observed [O/Fe] vs. [Fe/H] ratios of the thick
disc (magenta dots) with those predicted by model TD1 discussed in the text
(see Table4).

####  7.3.1 Models with stellar rotation

An improvement of the fit of the model TD1 in the [O/Fe] vs. [Fe/H] diagram
could be obtained adopting yields that account for rotational mixing in
massive stars, as recently suggested by Romano et al. (2019, 2020). In support
of this indication, we recall that, only yields including stellar rotation,
e.g. model MLr, are actually able to produce slopes similar to the ones
observed for the thick-disc stars (except for Mg; see Fig. 6). For this
purpose we calculated new models based on the chemical evolution parameters of
model TD1, but using the Limongi & Chieffi (2018) yields for different initial
rotational velocities, $V_{\rm rot}=0$, $V_{\rm rot}=150$, $V_{\rm
rot}=300\,\mathrm{km\,s^{-1}}$ and the averaged rotation yields, already
discussed in model MLr. The [O/Fe] vs. [Fe/H] diagram obtained with the new
models is shown in Fig. 9. Using yields from stars with higher rotational
velocity we obtain higher [O/Fe] ratios. This simply reflects the fact that
massive stars produce larger amount of oxygen at increasing $V_{\rm rot}$ (see
Figs. 13-16), while iron, being mainly contributed by SNIa, barely changes.
Figure 9 shows that enhancing the fraction of stars with high rotational
velocity may help explaining the observed [O/Fe], especially at metallicities
$\mbox{\rm[{\rm Fe}/{\rm H}]}\lesssim-1$. Looking at Fig. 9, we may see that
the non rotating model computed with Limongi & Chieffi (2018) yields fits the
lower envelope of the thin disc stars, of some thick disc outliers at
$\mbox{\rm[{\rm Fe}/{\rm H}]}\simeq-2$ and of the few halo stars with
$\mbox{\rm[{\rm Fe}/{\rm H}]}<-2$. By increasing the rotational velocity the
models shifts toward higher [Fe/H] values. In particular using the model
computed adopting the averaged yields obtained from the three different
rotational velocity values (cyan line, Prantzos et al. (2018)), our model fits
the border between thin and thick disc stars. Then, by further increasing the
rotational velocity, all the thick disc stars can be fitted. These models are
also able to reproduce the observed [O/Fe] values of the halo stars.

Figure 9: [O/Fe] vs. [Fe/H] abundances predicted by massive star models without and with rotation. The blue, orange and green lines refer to Limongi & Chieffi (2018) yields for initial rotational velocities $V_{\rm rot}$=0, $V_{\rm rot}$=150 and $V_{\rm rot}$=300 km/s, respectively. The cyan line represents the MLr model with rotation-averaged yields.  Table 4: IMF parameters of the chemical evolution models for the thick disc.$n_{*}$ denotes the corresponding number of stars born in the selected mass range per 106${M}_{\odot}$ of gas converted into stars. model | $\mbox{${M_{\rm UP}}$}\ [\mbox{${M}_{\odot}$}]$ | $x$ | $\frac{n_{*}}{1E6\mbox{${M}_{\odot}$}}$ | $\frac{n_{*}}{1E6\mbox{${M}_{\odot}$}}$
---|---|---|---|---
|  |  | (10-120${M}_{\odot}$) | ($120-\mbox{${M_{\rm UP}}$}$)
TD1 | 120 | 1.5 | 5448 | 0
TD2 | 120 | 1.7 | 3465 | 0
TD3 | 200 | 1.7 | 3449 | 30
TD4 | 200 | 1.4 | 6626 | 108
TD5 | 350 | 1.2 | 9226 | 356
all models: $\nu$=1.4Gyr-1; $k=1$; $\tau_{inf}$=0.5 Gyr; ASNIa=0.04

Figure 10: Same as Figure 7 but for the thick disc model TD4.

Figure 11: Comparison of the observed [O/Fe] vs. [Fe/H] (left panel) and
[Mg/Fe] vs. [Fe/H] (right panel) ratios of the thick disc stars, with the
predictions of models with different IMF parameters as listed in Table4. We
see that increasing ${M_{\rm UP}}$, so as to include the contributions of VMO
(winds, PPISN and PISN), thick-disc stars can be reproduced pretty well.
Moreover, changing both the IMF slope and ${M_{\rm UP}}$, different
populations can be recovered.

####  7.3.2 Effects due to IMF variations and VMO

Another possible explanation for the different [O/Fe] evolution of thin- and
thick-disc stars could be linked to different IMFs of the two populations. As
already discussed in Sect. 3.6, the IGIMF predicts that, at low metallicity
and high star formation rates, the IMF may become steeper at low masses and
flatter at high masses than the canonical one (Marks et al. 2012; Jeřábková et
al. 2018). This suggests that, in the early phases of the thick disc
evolution, the IMF could be well populated up to ${M_{\rm UP}}$. To
investigate this aspect, we make use of our new MTW set of yields, which
includes also the chemical contributions of VMO (winds, PPISN/PISN explosions
and DBH). This gives us the possibility to investigate the effects of changing
both the IMF exponent $x$ (for $\mbox{$M_{\rm i}$}>1\,\mbox{${M}_{\odot}$}$),
and the upper mass limit ${M_{\rm UP}}$, that can be pushed up to
$350\,\mbox{${M}_{\odot}$}$.

Starting from the chemical evolution parameters calibrated to reproduce the
MDF shown in Figure 7, we vary the IMF parameters ($x$ and ${M_{\rm UP}}$),
shifting the upper mass limit into the range ($120\leq\mbox{${M_{\rm
UP}}$}/\mbox{${M}_{\odot}$}\leq 350$). The adopted values, listed in Table 4,
are meant to test slopes typically found in the literature ($1.3\leq x\leq
1.7$ for $\mbox{$M_{\rm i}$}>1\,\mbox{${M}_{\odot}$}$; see Sect. 6.2), and to
obtain models that are able to bracket the observed thick disc data shown in
Figure 8. The theoretical thick-disc MDFs somewhat change with the IMF
variation (Fig. 10), however, the impact is not dramatic within the explored
ranges of the IMF parameters, and the corresponding models actually bracket
the median [Fe/H] of the observed distribution.

The abundance patterns of thick disc stars predicted by these new models are
compared with observations in the [O/Fe] and [Mg/Fe] vs. [Fe/H] diagrams shown
in Figure 11. We note that, in the [O/Fe] vs. [Fe/H] diagram, model TD2, with
x=1.7, runs almost parallel, below model TD1, eventually being able to
reproduce the branch with low [O/Fe] ratios at [Fe/H] =0. Models TD1 and TD2
are also able to reproduce the lower envelope of observed halo stars with
[Fe/H] $\geq$ -3. However, while below this metallicity, the data show a
plateau or even a decreasing branch, the [O/Fe] values of the above two models
continue to increase. We recall that the majority of these halo stars are sub-
giant or giant stars from Cayrel et al. (2004) and their [O/Fe] values have
been corrected for 3D effects, as suggested in Nissen et al. (2002).

The behaviour of the models change significantly if they include the ejecta
from VMO. In Table 4, we report also the expected number of stars per 106
${M}_{\odot}$ of gas converted into stars, that are born with $10\
\mbox{${M}_{\odot}$}\ \leq\mbox{$M_{\rm i}$}\ \leq\ 120$ ${M}_{\odot}$ and
with $120\ \mbox{${M}_{\odot}$}\ \leq\mbox{$M_{\rm i}$}\leq$ ${M_{\rm UP}}$,
resulting from the corresponding IMF parameters. It is easy to see from Table
4, that the fractional contribution of VMOs to the total number of massive
stars progressively increase as we move from model TD3 to model TD4 and TD5.
We first note that even a small fraction of stars born as VMO is enough to
change the initial evolution of the [O/Fe] ratios (model TD3). This is due to
the fact that low metallicity VMOs are able to produce enough Fe to decrease
the [O/Fe] ratio. Then, as the fractional number of VMOs increases, the models
are able to reach lower and lower [O/Fe] ratios. Models TD3, TD4 and TD5
reproduce well the distribution of halo stars with [Fe/H] $\leq$ -2.

When the metallicity rises above $Z_{\rm i}$ $\sim$ 0.0001 there is a jump in
the 16O production with respect to that of Fe. Indeed, while 16O is copiously
produced in almost all PISN models, for Fe this is true only for the more
massive ones, with Ej(Fe) $\sim$ 20 ${M}_{\odot}$ for ${M}_{\rm He}$ $\sim$
115 ${M}_{\odot}$ and Ej(Fe) $\sim$ 0.5 ${M}_{\odot}$ for ${M}_{\rm He}$
$\sim$ 85 ${M}_{\odot}$(Heger & Woosley 2002). This alternative behaviour of
16O and Fe production by PISNe at decreasing ${M}_{\rm He}$, gives rise to a
rapid jump in the [O/Fe] ratio at [Fe/H] $\sim$ -2. After the jump, the [O/Fe]
ratio reaches a peak well within the observed thick disc data, and then
decreases following the observed slope of the data. At [Fe/H] $\geq$ -2,
models TD3 and TD4 bracket the majority of the thick disc, data leaving
outside only a few objects with high [O/Fe] ratios. Model TD5 can be
considered as an extreme case, that is able to reproduce the lower and upper
envelopes of the [O/Fe] ratios, at the lower and higher metallicities,
respectively.

A similar evolution can be seen also in the case of the [Mg/Fe] vs. [Fe/H]
ratios, as shown in the right panel of Figure 11. We note that models without
VMOs, TD1 and TD2, are not able to reproduce the [Mg/Fe] ratios observed at
very low metallicity [Fe/H] $\leq$ -2. On the contrary, models TD3, TD4 and
TD5 run across the observed metal-poor halo stars, well reproducing the
growing [Mg/Fe] ratio at increasing [Fe/H]. Thus, in the very early stages of
the thick disc chemical evolution, the models with VMOs are able to reproduce
the observed $\alpha$-poor (Mg) stars, for instance, those with
$\mbox{\rm[{\rm Mg}/{\rm Fe}]}\sim 0.0$ at $\mbox{\rm[{\rm Fe}/{\rm
H}]}\lesssim-1$. For [Fe/H] $>$ -2, models TD3, TD4 and TD5 jump inside the
thick disc [Mg/Fe] observations, providing an overall good fit to the
evolution of halo and thick disc stars.

In summary, while no single TD model is able to reproduce the [O/Fe] and
[Mg/Fe] observations of the halo and thick disc stars at once, a combination
of TD models with varying IMF, in line with the predictions of the IGIMF
already discussed, is able to reproduce the location and the dispersion
observed among halo stars and the thick disc chemical evolution pattern.

Finally, we note that, the inclusion of PISNe may also help in solving the
issue of Mg underprediction present in most of the published yields tables
(Timmes et al. 1995; Portinari et al. 1998; Prantzos et al. 2018).

In conclusion, our analysis shows that IMF variations, in terms of slope but,
more importantly, of ${M_{\rm UP}}$, may significantly affect the predictions
of chemical evolution models. In turn, this may impact on our understanding of
the observed MW populations.

##  8 Summary and concluding remarks

In this work we investigate the impact of the ejecta of massive and very
massive stars on the predictions of chemical evolution models. We constructed
different sets of chemical yields that include both winds and explosion
contributions. To this aim, we collected various explosion yields available in
the literature and combine them with the wind ejecta computed with our code
PARSEC or other authors. A novelty of this work is that we investigate the
effect of increasing the IMF upper mass limit up to $\mbox{${M_{\rm
UP}}$}=350\,\mbox{${M}_{\odot}$}$, which is well beyond the standard value of
$\mbox{${M_{\rm UP}}$}\sim 100\,\mbox{${M}_{\odot}$}$, typically used in most
studies. This allows us to explore the impact of VMO, which are expected to
eject significant amounts of newly produced elements through both powerful
stellar winds and PPISN/PISN explosions. For completeness, we also include the
chemical yields of AGB stars computed with our COLIBRI (Marigo et al. 2013)
and from other studies. We note that the impact of AGB stars on the chemical
species considered in this work (Fe, O, Mg, Si, Ca) is minor. An analysis
aimed to highlight the role of the chemical yields from low- and intermediate-
mass stars is postponed to a dedicated future work.

The different sets of chemical yields are then incorporated in our chemical
evolution code CHE-EVO (Silva et al. 1998) to analyse the chemical evolution
of thin- and thick-disc stars in the solar vicinity, for which we can rely on
accurate abundance measurements and kinematical classification (Bensby et al.
2014).

For each set of chemical yields, we ran grids of chemical evolution models to
single out the parameters that best fit the main observational constraints of
the MW thin disc, namely: the present-day star formation rate and gas
fraction, the rates of CCSN and SNIa, the metallicity of the Sun at its birth
epoch $\simeq 4.6$ Gyr ago, and the observed MDF. For all sets of chemical
yields, we were able to find suitable combinations of input parameters that
match all the constraints reasonably well (see Table 2). In the best-fit
models the Schmidt star formation law, has a typical efficiency in the range
$0.4<\nu<0.8$ and an exponent $k=1$. The gas infall time-scale generally
varies in the range $3<\tau/\mathrm{Gyr}<6$.

Once the chemical evolution model is calibrated for each set of yields, we
move to test them against the observations. We focus on the trends of the
abundance ratios for a few -elements ([O/Fe], [Mg/Fe], [Si/Fe], [Ca/Fe]) as a
function of the metallicity, traced by [Fe/H].

A general agreement is found for [O/Fe] and [Si/Fe], while the predictions for
[Mg/Fe] and [Ca/Fe] run below the observed data. The best results are obtained
when including the yields from rotating massive stars (Limongi & Chieffi 2018)
and those based on our PARSEC models for the hydrostatic phases. With the
Ritter et al. (2018) yields, all abundance ratios are significantly below the
observed values, despite the fact that the corresponding calibrated chemical
evolution model reproduces the thin-disc constraints. After carefully
examining the problem, we conclude that a possible cause is the large iron
production of these sets of yields.

Most best-fit models are able to recover reasonably well the slope of the bulk
of the thin-disc stars. With our MTW set of yields, the calibrated IMF
exponent is $x=1.5$, which is intermediate between the results of Kroupa et
al. (1993, $x=1.7$) and the top-heavier IMFs of Kroupa (2001, $x=1.3$) and
Chabrier (2003, $x=1.3$).

In this study we did not find evidence that the IMF exponent is degenerate
with other chemical evolution parameters, e.g. $\nu$ or $k$, but we cannot
exclude that more complex chemical evolution models may provide a different
indication. Recently Valentini et al. (2019) pointed out that hydrodynamical
models support either of the Kroupa IMF slopes, depending on the set of
observational data adopted for comparison. In this respect, we note that a
steeper slope of the IMF is more suitable to reproduce thin-disc stars, as
suggested by various other investigations (Matteucci & Francois 1989;
Matteucci 2001; Grisoni et al. 2018; Matteucci et al. 2019, 2020).

Thick-disc stars exhibit a [O/Fe] vs. [Fe/H] pattern different from that of
the thin-disc stars, characterised by a steeper slope and a larger degree of
-enhancement, with higher [O/Fe] ratios. This feature is much less evident in
the [Mg/Fe] diagram. In this respect, we note that there are stars, likely
belonging to the thin disc according to the kinematical classification, that
has higher [O/Fe] values and fall in the region populated by the thick-disc
stars. Our thin-disc models do not reproduce such stars. Possible explanations
of this anomalous behaviour are that, either they have formed in different
environmental conditions compatible with a top-heavier IMF (e.g. in the IGIMF
scenario for a low metallicity and/or high SFR) or, rather, the kinematical
parameters of these stars varied in such a way that they are now classified as
thin-disc members.

None of the chemical evolution models calibrated on the thin disc is able to
reproduce the [O/Fe] trend of the thick disc. To reproduce the higher values
of the $\alpha$\- enhancement seen in thick disc stars one need to change the
chemical evolution parameters in particular reduce the infall time-scale to
$\tau_{inf}$=0.5 Gyr and increase the star formation efficiency to $\nu$=1.4
Gyr-1, in agreement with the literature (Grisoni et al. 2019, 2020a, 2020b).
The MDF of thick disc stars is well reproduced with these parameters but, in
order to better fit its fast decreasing tail at the higher metallicity with
our simple chemical evolution model, we need to impose a galactic wind at an
age of 2.5 Gyr. Tough the MDF is well reproduced, the tun on the [O/Fe] vs.
[Fe/H] diagram is not satisfactory. The model enters the region populated by
thick disc stars but with a slope that is flatter than that of the trend of
the data. In particular, the model cannot reproduce the high $\alpha$\-
enhanced thick disc stars at low metallicities. The model cannot reproduce
also the pattern of halo stars, which have been included here to check its
early chemical evolution.

Since it has been recently suggested that a higher $\alpha$-enhancement could
be obtained by considering yields from fast rotating stars (Limongi & Chieffi
2018; Romano et al. 2019), we have calculated new models changing only the
yields and using those that include effects of rotation (Limongi & Chieffi
2018; Prantzos et al. 2018). Our analysis shows that we can well reproduce the
[O/Fe] vs. [Fe/H] diagram of thick disc by increasing only the fraction of
rotating stars. The LC18 models without rotation run along the lower envelope
of thick disc while those with intermediate and high rotational velocity
reproduce its upper envelope as shown in Figure 9. This in fact may give rise
to a degeneracy between chemical evolution parameters ($\tau_{inf}$ and $\nu$)
and stellar rotational velocities. The low metallicity halo stars in this
diagram cannot be reproduced by our model.

An alternative way to reproduce the observations of the thick disc is offered
by the recent growing evidence for a significant dependence of the IMF on the
environment (Marks et al. 2012; Jeřábková et al. 2018). In the IGIMF theory
the IMF may become top heavier at decreasing metallicity and increasing star
formation rate. If this is true, the different MDFs of the thin and thick
discs could be not only the result but also one of the reasons of the
different chemical evolution patterns of the two MW components. In particular
in the early phases of the thick disc evolution, the IMF could be well
populated up to an ${M_{\rm UP}}$ that is able to switch on the chemical
contribution from PISNe. Motivated by the above reasons we have used our MTW
yields to analyzed the effects of varying both the slope and ${M_{\rm UP}}$.
Using the same chemical evolution parameters of model TD1, we computed other
models to check the effects of such IMF variations. Among the calculated new
models we have selected four representative cases that are able to match the
properties of thick disc stars. One with a steeper IMF (x=1.7) and canonical
${M_{\rm UP}}$ =120 which defines the lower boundary of the thick disc in the
[O/Fe] vs. [Fe/H] diagram, other two with ${M_{\rm UP}}$ =200 and x=1.7 and
1.4, respectively and a final one with ${M_{\rm UP}}$ =350 and x=1.2.

The latter three models, characterized by a ${M_{\rm UP}}$ that falls within
the VMOs, are able to bracket the thick disc observations in both the [O/Fe]
vs. [Fe/H]and [Mg/Fe] vs. [Fe/H] diagrams. At the same time, while with
canonical ${M_{\rm UP}}$ we are not able to reproduce the properties of metal-
poor halo stars in neither of the above diagrams, this is possible with the
models that include the contribution of PISNe. Halo stars, in the [O/Fe] and
[Mg/Fe] vs. [Fe/H] diagrams, show a constant or even decreasing trend of the
$\alpha$-enhancement at decreasing [Fe/H]. This trend can be reproduced if the
models include enough Fe contribution coming from PISNe at lower metallicities
([Fe/H] $\leq$-2). While there is a dispute concerning the reality of this
decreasing branch in the [O/Fe] diagram (Nissen et al. 2002; Cayrel et al.
2004; Micali et al. 2013) this constant or even decreasing trend is seen also
in the [Mg/Fe] diagram. In any case recent analysis of O and Fe measurements
in disc and halo stars by Amarsi et al. (2019) indicate that the [O/Fe] ratio
reaches a constant plateau of [O/Fe] $\sim$0.6 below [Fe/H] $\sim$-1 and it is
thus not compatible with the constant growth, at decreasing [Fe/H], indicated
by our models with canonical ${M_{\rm UP}}$. It is worth to note here that the
models that include the effects of rotation are facing the same problem.
Instead, an increasing fraction of VMOs, eventually in line with the
prediction of the IGIMF theory, will produce models that progressively show
lower [$\alpha$/Fe] ratios, thus naturally explaining the halo stars
dispersion and, at the same time, matching the high [$\alpha$/Fe] ratios of
thick disc stars.

In conclusion, our study provides a clear indication that PISN may have played
a significant role in shaping the chemical evolution path of halo and thick-
disc stars. Especially in view of the possible flattening of the IMF at low
metallicity and high SFR, their contribution should not be neglected in
chemical evolution models. In the light of these results, we may reasonably
expect that similar data, which actually exist for nearby extremely metal-poor
star bursting galaxies (Kojima et al. 2020), might host the chemical signature
of very massive stars, and that would witness their key role in the early
phases of galaxy evolution.

The complete TW sets of chemical ejecta for massive and very massive stars,
described in Appendix A, are made publicly available through the URL: xxxxx.

###### Acknowledgements.

We thank F. Matteucci, D. Romano and the anonymous referee for useful
suggestions that helped to improve the manuscript. PM and AS acknowledge the
support from the ERC Consolidator Grant funding scheme (project STARKEY, grant
agreement n. 615604). AB and LS acknowledge PRIN MIUR 2017 - 20173ML3WW_001.
MS acknowledges funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie-Skłodowska-Curie grant agreement No.
794393. This research made use of Astropy, a community-developed core Python
package for Astronomy (Astropy Collaboration et al. 2013, 2018) and
matplotlib, a Python library for publication quality graphics (Hunter 2007).
This work has been partially supported by PRIN MIUR 2017 prot. 20173ML3WW 002,
‘Opening the ALMA window on the cosmic evolution of gas, stars and
supermassive black holes’. AL acknowledges the MIUR grant ‘Finanziamento
annuale individuale attivita base di ricerca’ and the EU H2020-MSCA-ITN-2019
Project 860744 ‘BiD4BEST: Big Data applications for Black hole Evolution
STudies’. VG acknowledges financial support at SISSA from the European Social
Fund operational Programme 2014/2020 of the autonomous region Friuli Venezia
Giulia.

## References

  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2020) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, Phys. Rev. Lett., 125, 101102
  * Allende Prieto et al. (2008) Allende Prieto, C., Majewski, S. R., Schiavon, R., et al. 2008, Astronomische Nachrichten, 329, 1018
  * Amarsi et al. (2019) Amarsi, A. M., Nissen, P. E., & Skúladóttir, Á. 2019, A&A, 630, A104
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  * Aumer et al. (2017) Aumer, M., Binney, J., & Schönrich, R. 2017, MNRAS, 470, 3685
  * Ballero et al. (2006) Ballero, S. K., Matteucci, F., & Chiappini, C. 2006, New A, 11, 306
  * Bekki & Tsujimoto (2011) Bekki, K. & Tsujimoto, T. 2011, ApJ, 738, 4
  * Belokurov et al. (2018) Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611
  * Bensby et al. (2017) Bensby, T., Feltzing, S., Gould, A., et al. 2017, A&A, 605, A89
  * Bensby et al. (2014) Bensby, T., Feltzing, S., & Oey, M. S. 2014, A&A, 562, A71
  * Boeche et al. (2013) Boeche, C., Chiappini, C., Minchev, I., et al. 2013, A&A, 553, A19
  * Bond et al. (1982) Bond, J. R., Arnett, W. D., & Carr, B. J. 1982, in NATO Advanced Science Institutes (ASI) Series C, Vol. 90, NATO Advanced Science Institutes (ASI) Series C, ed. M. J. Rees & R. J. Stoneham, 303–311
  * Bond et al. (1984) Bond, J. R., Arnett, W. D., & Carr, B. J. 1984, ApJ, 280, 825
  * Botelho et al. (2020) Botelho, R. B., Milone, A. d. C., Meléndez, J., et al. 2020, MNRAS, 499, 2196
  * Bressan et al. (1994) Bressan, A., Chiosi, C., & Fagotto, F. 1994, ApJS, 94, 63
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127
  * Burrows et al. (2018) Burrows, A., Vartanyan, D., Dolence, J. C., Skinner, M. A., & Radice, D. 2018, Space Sci. Rev., 214, 33
  * Caffau et al. (2011) Caffau, E., Ludwig, H.-G., Steffen, M., Freytag, B., & Bonifacio, P. 2011, Sol. Phys., 268, 255
  * Casagrande et al. (2011) Casagrande, L., Schönrich, R., Asplund, M., et al. 2011, A&A, 530, A138
  * Cayrel et al. (2004) Cayrel, R., Depagne, E., Spite, M., et al. 2004, A&A, 416, 1117
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Chen et al. (2014) Chen, K.-J., Woosley, S., Heger, A., Almgren, A., & Whalen, D. J. 2014, ApJ, 792, 28
  * Chen et al. (2015) Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068
  * Cherchneff & Dwek (2009) Cherchneff, I. & Dwek, E. 2009, ApJ, 703, 642
  * Cherchneff & Dwek (2010) Cherchneff, I. & Dwek, E. 2010, ApJ, 713, 1
  * Chiappini (2009) Chiappini, C. 2009, in IAU Symposium, Vol. 254, The Galaxy Disk in Cosmological Context, ed. J. Andersen, Nordströara, B. m, & J. Bland -Hawthorn, 191–196
  * Chiappini et al. (1997) Chiappini, C., Matteucci, F., & Gratton, R. 1997, ApJ, 477, 765
  * Chiappini et al. (2001) Chiappini, C., Matteucci, F., & Romano, D. 2001, ApJ, 554, 1044
  * Chieffi & Limongi (2004) Chieffi, A. & Limongi, M. 2004, ApJ, 608, 405
  * Chiosi (1980) Chiosi, C. 1980, A&A, 83, 206
  * Chiosi et al. (1998) Chiosi, C., Bressan, A., Portinari, L., & Tantalo, R. 1998, A&A, 339, 355
  * Costa et al. (2020) Costa, G., Bressan, A., Mapelli, M., et al. 2020, arXiv e-prints, arXiv:2010.02242
  * Costa et al. (2019a) Costa, G., Girardi, L., Bressan, A., et al. 2019a, A&A, 631, A128
  * Costa et al. (2019b) Costa, G., Girardi, L., Bressan, A., et al. 2019b, MNRAS, 485, 4641
  * Cristallo et al. (2011) Cristallo, S., Piersanti, L., Straniero, O., et al. 2011, ApJS, 197, 17
  * Cristallo et al. (2009) Cristallo, S., Straniero, O., Gallino, R., et al. 2009, ApJ, 696, 797
  * Cristallo et al. (2015) Cristallo, S., Straniero, O., Piersanti, L., & Gobrecht, D. 2015, ApJS, 219, 40
  * Crowther (2019) Crowther, P. A. 2019, Galaxies, 7, 88
  * Crowther et al. (2016) Crowther, P. A., Caballero-Nieves, S. M., Bostroem, K. A., et al. 2016, MNRAS, 458, 624
  * de Jager et al. (1988) de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259
  * de Laverny et al. (2013) de Laverny, P., Recio-Blanco, A., Worley, C. C., et al. 2013, The Messenger, 153, 18
  * Doherty et al. (2014) Doherty, C. L., Gil-Pons, P., Lau, H. H. B., Lattanzio, J. C., & Siess, L. 2014, MNRAS, 437, 195
  * Doherty et al. (2017) Doherty, C. L., Gil-Pons, P., Siess, L., & Lattanzio, J. C. 2017, PASA, 34, e056
  * Ertl et al. (2016) Ertl, T., Janka, H.-T., Woosley, S. E., Sukhbold, T., & Ugliano, M. 2016, ApJ, 818, 124
  * Evans et al. (2010) Evans, C. J., Walborn, N. R., Crowther, P. A., et al. 2010, ApJ, 715, L74
  * Farmer et al. (2019) Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, ApJ, 887, 53
  * Feltzing et al. (2003) Feltzing, S., Bensby, T., & Lundström, I. 2003, A&A, 397, L1
  * Ferrini et al. (1994) Ferrini, F., Molla, M., Pardi, M. C., & Diaz, A. I. 1994, ApJ, 427, 745
  * Fontanot et al. (2009) Fontanot, F., Somerville, R. S., Silva, L., Monaco, P., & Skibba, R. 2009, MNRAS, 392, 553
  * Fryer (1999) Fryer, C. L. 1999, ApJ, 522, 413
  * Fryer et al. (2012) Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91
  * Fryer & Kalogera (2001) Fryer, C. L. & Kalogera, V. 2001, ApJ, 554, 548
  * Fryer et al. (2006) Fryer, C. L., Young, P. A., & Hungerford, A. L. 2006, ApJ, 650, 1028
  * Fuhrmann (2011) Fuhrmann, K. 2011, MNRAS, 414, 2893
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1
  * García-Berro et al. (1997) García-Berro, E., Ritossa, C., & Iben, Icko, J. 1997, ApJ, 485, 765
  * Gilmore et al. (2012) Gilmore, G., Randich, S., Asplund, M., et al. 2012, The Messenger, 147, 25
  * Girardi et al. (2010) Girardi, L., Williams, B. F., Gilbert, K. M., et al. 2010, ApJ, 724, 1030
  * Gräfener (2008) Gräfener, G. 2008, in Clumping in Hot-Star Winds, ed. W.-R. Hamann, A. Feldmeier, & L. M. Oskinova, 103
  * Grand et al. (2018) Grand, R. J. J., Bustamante, S., Gómez, F. A., et al. 2018, MNRAS, 474, 3629
  * Grisoni et al. (2020a) Grisoni, V., Cescutti, G., Matteucci, F., et al. 2020a, MNRAS, 492, 2828
  * Grisoni et al. (2019) Grisoni, V., Matteucci, F., Romano, D., & Fu, X. 2019, MNRAS, 489, 3539
  * Grisoni et al. (2020b) Grisoni, V., Romano, D., Spitoni, E., et al. 2020b, MNRAS, 498, 1252
  * Grisoni et al. (2018) Grisoni, V., Spitoni, E., & Matteucci, F. 2018, MNRAS, 481, 2570
  * Grisoni et al. (2017) Grisoni, V., Spitoni, E., Matteucci, F., et al. 2017, MNRAS, 472, 3637
  * Haemmerlé et al. (2018) Haemmerlé, L., Woods, T. E., Klessen, R. S., Heger, A., & Whalen, D. J. 2018, MNRAS, 474, 2757
  * Haft et al. (1994) Haft, M., Raffelt, G., & Weiss, A. 1994, ApJ, 425, 222
  * Haywood et al. (2018) Haywood, M., Di Matteo, P., Lehnert, M. D., et al. 2018, ApJ, 863, 113
  * Heger et al. (2003a) Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003a, ApJ, 591, 288
  * Heger et al. (2003b) Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003b, ApJ, 591, 288
  * Heger & Woosley (2002) Heger, A. & Woosley, S. E. 2002, ApJ, 567, 532
  * Helmi et al. (2018) Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85
  * Herwig (2000) Herwig, F. 2000, A&A, 360, 952
  * Herwig (2005) Herwig, F. 2005, ARA&A, 43, 435
  * Hosek et al. (2019) Hosek, Matthew W., J., Lu, J. R., Anderson, J., et al. 2019, ApJ, 870, 44
  * Hunt et al. (2019) Hunt, L. K., De Looze, I., Boquien, M., et al. 2019, A&A, 621, A51
  * Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543
  * Iben et al. (1997) Iben, Icko, J., Ritossa, C., & García-Berro, E. 1997, ApJ, 489, 772
  * Iglesias & Rogers (1996) Iglesias, C. A. & Rogers, F. J. 1996, ApJ, 464, 943
  * Itoh & Kohyama (1983) Itoh, N. & Kohyama, Y. 1983, ApJ, 275, 858
  * Itoh et al. (2008) Itoh, N., Uchida, S., Sakamoto, Y., Kohyama, Y., & Nozawa, S. 2008, ApJ, 677, 495
  * Iwamoto et al. (1999) Iwamoto, K., Brachwitz, F., Nomoto, K., et al. 1999, ApJS, 125, 439
  * Izzo et al. (2019) Izzo, L., de Ugarte Postigo, A., Maeda, K., et al. 2019, Nature, 565, 324
  * Janka (2012) Janka, H.-T. 2012, Annual Review of Nuclear and Particle Science, 62, 407
  * Jeřábková et al. (2018) Jeřábková, T., Hasani Zonoozi, A., Kroupa, P., et al. 2018, A&A, 620, A39
  * Jimenez et al. (1998) Jimenez, R., Padoan, P., Matteucci, F., & Heavens, A. F. 1998, MNRAS, 299, 123
  * Jurić et al. (2008) Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864
  * Kalogera & Baym (1996) Kalogera, V. & Baym, G. 1996, ApJ, 470, L61
  * Karakas (2010) Karakas, A. I. 2010, MNRAS, 403, 1413
  * Karakas & Lugaro (2016) Karakas, A. I. & Lugaro, M. 2016, ApJ, 825, 26
  * Karakas et al. (2018) Karakas, A. I., Lugaro, M., Carlos, M., et al. 2018, MNRAS, 477, 421
  * Karlsson et al. (2008) Karlsson, T., Johnson, J. L., & Bromm, V. 2008, The Astrophysical Journal, 679, 6
  * Kawata & Chiappini (2016) Kawata, D. & Chiappini, C. 2016, Astronomische Nachrichten, 337, 976
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, ApJ, 498, 541
  * Kojima et al. (2020) Kojima, T., Ouchi, M., Rauch, M., et al. 2020, arXiv e-prints, arXiv:2006.03831
  * Kozyreva et al. (2014a) Kozyreva, A., Blinnikov, S., Langer, N., & Yoon, S.-C. 2014a, A&A, 565, A70
  * Kozyreva et al. (2014b) Kozyreva, A., Yoon, S.-C., & Langer, N. 2014b, A&A, 566, A146
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
  * Kroupa (2008) Kroupa, P. 2008, in Astronomical Society of the Pacific Conference Series, Vol. 390, Pathways Through an Eclectic Universe, ed. J. H. Knapen, T. J. Mahoney, & A. Vazdekis, 3
  * Kroupa et al. (1993) Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545
  * Kroupa & Weidner (2003) Kroupa, P. & Weidner, C. 2003, ApJ, 598, 1076
  * Kubryk et al. (2015) Kubryk, M., Prantzos, N., & Athanassoula, E. 2015, A&A, 580, A126
  * Larson (1972) Larson, R. B. 1972, Nature, 236, 21
  * Lawlor et al. (2008) Lawlor, T. M., Young, T. R., Johnson, T. A., & MacDonald, J. 2008, MNRAS, 384, 1533
  * Lebzelter et al. (2018) Lebzelter, T., Mowlavi, N., Marigo, P., et al. 2018, A&A, 616, L13
  * Leung et al. (2019) Leung, S.-C., Nomoto, K., & Blinnikov, S. 2019, The Astrophysical Journal, 887, 72
  * Leung et al. (2020) Leung, S.-C., Nomoto, K., & Suzuki, T. 2020, ApJ, 889, 34
  * Limongi & Chieffi (2003) Limongi, M. & Chieffi, A. 2003, ApJ, 592, 404
  * Limongi & Chieffi (2006) Limongi, M. & Chieffi, A. 2006, ApJ, 647, 483
  * Limongi & Chieffi (2018) Limongi, M. & Chieffi, A. 2018, VizieR Online Data Catalog, J/ApJS/237/13
  * Lo Faro et al. (2013) Lo Faro, B., Franceschini, A., Vaccari, M., et al. 2013, ApJ, 762, 108
  * Lo Faro et al. (2015) Lo Faro, B., Silva, L., Franceschini, A., Miller, N., & Efstathiou, A. 2015, MNRAS, 447, 3442
  * Lynden-Bell (1975) Lynden-Bell, D. 1975, Vistas in Astronomy, 19, 299
  * Majewski et al. (2017) Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94
  * Marchant & Moriya (2020) Marchant, P. & Moriya, T. J. 2020, A&A, 640, L18
  * Marigo (2001) Marigo, P. 2001, A&A, 370, 194
  * Marigo & Aringer (2009) Marigo, P. & Aringer, B. 2009, A&A, 508, 1539
  * Marigo et al. (2013) Marigo, P., Bressan, A., Nanni, A., Girardi, L., & Pumo, M. L. 2013, MNRAS, 434, 488
  * Marigo et al. (2003) Marigo, P., Chiosi, C., & Kudritzki, R. P. 2003, A&A, 399, 617
  * Marigo et al. (2020) Marigo, P., Cummings, J. D., Curtis, J. L., et al. 2020, Nature Astronomy [arXiv:2007.04163]
  * Marigo et al. (2017) Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77
  * Marigo et al. (2001) Marigo, P., Girardi, L., Chiosi, C., & Wood, P. R. 2001, A&A, 371, 152
  * Marks et al. (2012) Marks, M., Kroupa, P., Dabringhausen, J., & Pawlowski, M. S. 2012, MNRAS, 422, 2246
  * Matteucci (1994) Matteucci, F. 1994, A&A, 288, 57
  * Matteucci (2001) Matteucci, F. 2001, The chemical evolution of the Galaxy, Vol. 253
  * Matteucci (2014) Matteucci, F. 2014, Saas-Fee Advanced Course, 37, 145
  * Matteucci & Brocato (1990) Matteucci, F. & Brocato, E. 1990, ApJ, 365, 539
  * Matteucci & Francois (1989) Matteucci, F. & Francois, P. 1989, MNRAS, 239, 885
  * Matteucci & Greggio (1986) Matteucci, F. & Greggio, L. 1986, A&A, 154, 279
  * Matteucci et al. (2019) Matteucci, F., Grisoni, V., Spitoni, E., et al. 2019, MNRAS, 487, 5363
  * Matteucci & Pipino (2005) Matteucci, F. & Pipino, A. 2005, MNRAS, 357, 489
  * Matteucci et al. (2020) Matteucci, F., Vasini, A., Grisoni, V., & Schultheis, M. 2020, MNRAS, 494, 5534
  * Micali et al. (2013) Micali, A., Matteucci, F., & Romano, D. 2013, MNRAS, 436, 1648
  * Mikolaitis et al. (2017) Mikolaitis, Š., de Laverny, P., Recio-Blanco, A., et al. 2017, A&A, 600, A22
  * Munakata et al. (1985) Munakata, H., Kohyama, Y., & Itoh, N. 1985, ApJ, 296, 197
  * Nissen et al. (2002) Nissen, P. E., Primas, F., Asplund, M., & Lambert, D. L. 2002, A&A, 390, 235
  * Nomoto (1984) Nomoto, K. 1984, ApJ, 277, 791
  * Nomoto et al. (2013) Nomoto, K., Kobayashi, C., & Tominaga, N. 2013, ARA&A, 51, 457
  * Nomoto et al. (2006) Nomoto, K., Tominaga, N., Umeda, H., Kobayashi, C., & Maeda, K. 2006, Nuclear Physics A, 777, 424
  * O’Connor & Ott (2011) O’Connor, E. & Ott, C. D. 2011, ApJ, 730, 70
  * Ohkubo et al. (2009) Ohkubo, T., Nomoto, K., Umeda, H., Yoshida, N., & Tsuruta, S. 2009, ApJ, 706, 1184
  * Ohkubo et al. (2006) Ohkubo, T., Umeda, H., Maeda, K., et al. 2006, ApJ, 645, 1352
  * Padoan et al. (1997) Padoan, P., Nordlund, A., & Jones, B. J. T. 1997, MNRAS, 288, 145
  * Pagel & Edmunds (1981) Pagel, B. E. J. & Edmunds, M. G. 1981, ARA&A, 19, 77
  * Pastorelli et al. (2020) Pastorelli, G., Marigo, P., Girardi, L., et al. 2020, MNRAS, 498, 3283
  * Pastorelli et al. (2019) Pastorelli, G., Marigo, P., Girardi, L., et al. 2019, MNRAS, 485, 5666
  * Pignatari et al. (2016) Pignatari, M., Herwig, F., Hirschi, R., et al. 2016, ApJS, 225, 24
  * Plevne et al. (2020) Plevne, O., Önal Tas, Ö., Bilir, S., & Seabroke, G. M. 2020, ApJ, 893, 108
  * Poelarends et al. (2008) Poelarends, A. J. T., Herwig, F., Langer, N., & Heger, A. 2008, ApJ, 675, 614
  * Portinari & Chiosi (1999) Portinari, L. & Chiosi, C. 1999, A&A, 350, 827
  * Portinari et al. (1998) Portinari, L., Chiosi, C., & Bressan, A. 1998, A&A, 334, 505
  * Prantzos et al. (2018) Prantzos, N., Abia, C., Limongi, M., Chieffi, A., & Cristallo, S. 2018, MNRAS, 476, 3432
  * Prantzos & Aubert (1995) Prantzos, N. & Aubert, O. 1995, A&A, 302, 69
  * Prantzos et al. (2011) Prantzos, N., Boehm, C., Bykov, A. M., et al. 2011, Rev. Mod. Phys., 83, 1001
  * Prochaska et al. (2000) Prochaska, J. X., Naumov, S. O., Carney, B. W., McWilliam, A., & Wolfe, A. M. 2000, AJ, 120, 2513
  * Recchi et al. (2009) Recchi, S., Calura, F., & Kroupa, P. 2009, A&A, 499, 711
  * Renzo et al. (2020) Renzo, M., Farmer, R. J., Justham, S., et al. 2020, MNRAS, 493, 4333
  * Ricotti & Ostriker (2004) Ricotti, M. & Ostriker, J. P. 2004, MNRAS, 350, 539
  * Ritossa et al. (1996) Ritossa, C., Garcia-Berro, E., & Iben, Icko, J. 1996, ApJ, 460, 489
  * Ritter et al. (2018) Ritter, C., Herwig, F., Jones, S., et al. 2018, MNRAS, 480, 538
  * Robitaille & Whitney (2010) Robitaille, T. P. & Whitney, B. A. 2010, ApJ, 710, L11
  * Rollinde et al. (2009) Rollinde, E., Vangioni, E., Maurin, D., et al. 2009, MNRAS, 398, 1782
  * Romano et al. (2020) Romano, D., Franchini, M., Grisoni, V., et al. 2020, A&A, 639, A37
  * Romano et al. (2010) Romano, D., Karakas, A. I., Tosi, M., & Matteucci, F. 2010, A&A, 522, A32
  * Romano et al. (2019) Romano, D., Matteucci, F., Zhang, Z.-Y., Ivison, R. J., & Ventura, P. 2019, MNRAS, 490, 2838
  * Romano et al. (2017) Romano, D., Matteucci, F., Zhang, Z. Y., Papadopoulos, P. P., & Ivison, R. J. 2017, MNRAS, 470, 401
  * Rosenfield et al. (2014) Rosenfield, P., Marigo, P., Girardi, L., et al. 2014, ApJ, 790, 22
  * Rosenfield et al. (2016) Rosenfield, P., Marigo, P., Girardi, L., et al. 2016, ApJ, 822, 73
  * Sahijpal (2014) Sahijpal, S. 2014, International Journal of Astronomy and Astrophysics, 4, 491
  * Savino et al. (2020) Savino, A., Koch, A., Prudil, Z., Kunder, A., & Smolec, R. 2020, A&A, 641, A96
  * Schmidt (1959) Schmidt, M. 1959, ApJ, 129, 243
  * Schneider et al. (2018) Schneider, F. R. N., Sana, H., Evans, C. J., et al. 2018, Science, 359, 69
  * Schönrich & Aumer (2017) Schönrich, R. & Aumer, M. 2017, MNRAS, 472, 3979
  * Senchyna et al. (2020) Senchyna, P., Stark, D. P., Charlot, S., et al. 2020, arXiv e-prints, arXiv:2008.09780
  * Siess (2006) Siess, L. 2006, A&A, 448, 717
  * Siess (2007) Siess, L. 2007, A&A, 476, 893
  * Silva et al. (1998) Silva, L., Granato, G. L., Bressan, A., & Danese, L. 1998, ApJ, 509, 103
  * Silva et al. (2011) Silva, L., Schurer, A., Granato, G. L., et al. 2011, MNRAS, 410, 2043
  * Slemer (2016) Slemer, A. 2016, Investigating chemical enrichment scenarios fro massive stars, PhD Thesis, Department of Physics and Astronomy, Padova University
  * Slemer et al. (2017) Slemer, A., Marigo, P., Piatti, D., et al. 2017, MNRAS, 465, 4817
  * Snaith et al. (2014) Snaith, O. N., Haywood, M., Di Matteo, P., et al. 2014, ApJ, 781, L31
  * Spera & Mapelli (2017) Spera, M. & Mapelli, M. 2017, MNRAS, 470, 4739
  * Spera et al. (2015) Spera, M., Mapelli, M., & Bressan, A. 2015, MNRAS, 451, 4086
  * Spitoni et al. (2021) Spitoni, E., Verma, K., Silva Aguirre, V., et al. 2021, arXiv e-prints, arXiv:2101.08803
  * Sukhbold et al. (2016) Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H. T. 2016, ApJ, 821, 38
  * Tang et al. (2014) Tang, J., Bressan, A., Rosenfield, P., et al. 2014, MNRAS, 445, 4287
  * Tews & Schwenk (2020) Tews, I. & Schwenk, A. 2020, ApJ, 892, 14
  * Thielemann et al. (2011) Thielemann, F.-K., Hirschi, R., Liebendörfer, M., & Diehl, R. 2011, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 812, Lecture Notes in Physics, Berlin Springer Verlag, ed. R. Diehl, D. H. Hartmann, & N. Prantzos, 153–232
  * Thomas et al. (2005) Thomas, D., Maraston, C., Bender, R., & Mendes de Oliveira, C. 2005, ApJ, 621, 673
  * Timmes et al. (1995) Timmes, F. X., Woosley, S. E., & Weaver, T. A. 1995, ApJS, 98, 617
  * Ugliano et al. (2012) Ugliano, M., Janka, H.-T., Marek, A., & Arcones, A. 2012, ApJ, 757, 69
  * Umeda & Nomoto (2002) Umeda, H. & Nomoto, K. 2002, ApJ, 565, 385
  * Umeda & Nomoto (2008) Umeda, H. & Nomoto, K. 2008, ApJ, 673, 1014
  * Utrobin & Chugai (2009) Utrobin, V. P. & Chugai, N. N. 2009, A&A, 506, 829
  * Utrobin et al. (2010) Utrobin, V. P., Chugai, N. N., & Botticella, M. T. 2010, ApJ, 723, L89
  * Valentini et al. (2019) Valentini, M., Borgani, S., Bressan, A., et al. 2019, MNRAS, 485, 1384
  * van Dokkum et al. (2017) van Dokkum, P., Conroy, C., Villaume, A., Brodie, J., & Romanowsky, A. J. 2017, The Astrophysical Journal, 841, 68
  * van Dokkum & Conroy (2012) van Dokkum, P. G. & Conroy, C. 2012, ApJ, 760, 70
  * Vega et al. (2008) Vega, O., Clemens, M. S., Bressan, A., et al. 2008, A&A, 484, 631
  * Ventura & D’Antona (2009) Ventura, P. & D’Antona, F. 2009, A&A, 499, 835
  * Ventura et al. (2013) Ventura, P., Di Criscienzo, M., Carini, R., & D’Antona, F. 2013, MNRAS, 431, 3642