11institutetext: Chiaki Kobayashi 22institutetext: Centre for Astrophysics
Research, Department of Physics, Astronomy and Mathematics University of
Hertfordshire, College Lane, Hatfield AL10 9AB, UK 22email:
<EMAIL_ADDRESS>33institutetext: Philip Taylor 44institutetext:
Research School of Astronomy and Astrophysics, Australian National University,
Canberra, ACT 2611, Australia

# Chemo-Dynamical Evolution of Galaxies

Chiaki Kobayashi and Philip Taylor Chiaki Kobayashi

###### Abstract

Stars are fossils that retain the history of their host galaxies. Elements
heavier than helium are created inside stars and are ejected when they die.
From the spatial distribution of elements in galaxies, it is therefore
possible to constrain the physical processes during galaxy formation and
evolution. This approach, Galactic archaeology, has been popularly used for
our Milky Way Galaxy with a vast amount of data from Gaia satellite and multi-
object spectrographs to understand the origins of sub-structures of the Milky
Way. Thanks to integral field units, this approach can also be applied to
external galaxies from nearby to distant universe with the James Webb Space
Telescope. In order to interpret these observational data, it is necessary to
compare with theoretical predictions, namely chemodynamical simulations of
galaxies, which include detailed chemical enrichment into hydrodynamical
simulations from cosmological initial conditions. These simulations can
predict the evolution of internal structures (e.g., metallicity radial
gradients) as well as that of scaling relations (e.g., the mass-metallicity
relations). After explaining the formula and assumptions, we will show some
example results, and discuss future prospects.

##  1 Introduction

Figure 1:  Schematic view of chemical enrichment in galaxies.

Explaining the origin of the elements is one of the scientific triumphs
linking nuclear physics with astrophysics. As Fred Hoyle predicted, carbon and
heavier elements (‘metals’ in astrophysics) were not produced during the Big
Bang but instead created inside stars. So-called $\alpha$ elements (O, Ne, Mg,
Si, S, Ar, and Ca) are mainly produced by core-collapse supernovae, while
iron-peak elements (Cr, Mn, Fe, and Ni) are more produced by thermonuclear
explosions, observed as Type Ia supernovae (SNe Ia; Kobayashi et al. 2020a,
hereafter K20). The production depends on the mass of white dwarf (WD)
progenitors, and a large fraction of SNe Ia should come from near-
Chandrasekhar (Ch) mass explosions (see Kobayashi et al. 2020b for
constraining the relative contribution between near-Ch and sub-Ch mass SNe
Ia). Among core-collapse supernovae, hypernovae ($>\kern-11.99998pt\lower
4.73611pt\hbox{$\sim$}~{}10^{52}$ erg) produce a significant amount of Fe as
well as Co and Zn, and a significant fraction of massive stars
($>\kern-11.99998pt\lower 4.73611pt\hbox{$\sim$}~{}20M_{\odot}$) should
explode as hypernovae in order to explain the Galactic chemical evolution
(GCE; Kobayashi et al. 2006, hereafter K06).

Heavier elements are produced by neutron-capture processes. The slow neutron-
capture process (s-process) occurs in asymptotic giant branch (AGB) stars
(e.g., Busso et al. 1999; Herwig 2005; Cristallo et al. 2011; Karakas &
Lattanzio 2014), while the astronomical sites of rapid neutron-capture process
(r-process) have been debated. The possible sites are neutron-star (NS)
mergers (NSMs, Wanajo et al. 2014; Just et al. 2015), magneto-rotational
supernovae (MRSNe, Nishimura et al. 2015; Mösta et al. 2018; Reichert et al.
2021), magneto-rotational hypernovae (MRHNe, Yong et al. 2021), accretion
disks/collapsars (Siegel et al. 2019), and common envelope jets supernovae
(Grichener et al. 2021). Light neutron-capture elements (e.g., Sr) are also
produced by electron-capture supernovae (ECSNe, Wanajo et al. 2013),
$\nu$-driven winds (Arcones et al. 2007; Wanajo 2013), and rotating-massive
stars (Frischknecht et al. 2016; Limongi & Chieffi 2018).

The cycles of chemical enrichment are schematically shown in Figure 1, where
each cycle produces different elements and isotopes with different timescales.
In a galaxy, not only the total amount of metals, i.e. metallicity
$Z\equiv\Sigma_{i>{\rm He}}\,m_{i}/m$, but also elemental abundance ratios,
[X/Fe]$\equiv{\log N_{X}/N_{\rm Fe}/(N_{X,\odot}/N_{{\rm Fe},\odot}})$, evolve
as a function of time. Therefore, we can use all of this information as
fossils to study the formation and evolutionary histories of the galaxy. This
approach is called Galactic archaeology, and several on-going and future
surveys with multi-object spectrographs (MOS; e.g., APOGEE, HERMES-GALAH,
Gaia-ESO, DESI, WEAVE, 4MOST, MOONS, Subaru Prime Focus Spectrograph (PFS),
and Maunakea Spectroscopic Explorer (MSE)) are producing a vast amount of
observational data of elemental abundances. Moreover, integral field unit
(IFU) spectrographs (e.g., SAURON, SINFONI, CALIFA, SAMI, MaNGA, KMOS, MUSE,
HECTOR, and NIRSpec on JWST) allow us to measure metallicity and some
elemental abundance ratios within galaxies. It is now possible to apply the
same approach not only to our own Milky Way but also to other types of
galaxies or distant galaxies. Let us call this extra-galactic archaeology.

While the evolution of the dark matter in the standard $\Lambda$ cold dark
matter (CDM) cosmology is reasonably well understood, how galaxies form and
evolve is still much less certain because of the complexity of the baryon
physics such as star formation and feedback. The thermal energy ejected from
supernovae to the interstellar medium (ISM) suppresses star formation, while
the production of heavy elements in these supernovae enhances gas cooling.
Since these processes affect each other, galaxy formation and evolution is
complicated, and has to be solved consistently with a numerical simulation.
Feedback from central super-massive black holes (SMBHs) is also found to be
very important for explaining the observed properties of galaxies (e.g.,
Croton et al. 2006; Springel et al. 2005; Hopkins et al. 2006; Dubois et al.
2012; Taylor & Kobayashi 2014).

Since the 90’s, the development of high performance computers and
computational techniques has made it possible to simulate the formation and
evolutionary history of galaxies, not only of isolated systems (e.g., Burkert
& Hensler 1987; Katz 1992; Steinmetz & Mueller 1994; Mihos & Hernquist 1996;
Kawata & Gibson 2003; Kobayashi 2004) but also for cosmological simulations of
individual galaxies (e.g., Navarro & White 1994) or of the galaxy population
as a whole (e.g., Cen & Ostriker 1999; Kobayashi et al. 2007). Thanks to the
public release of hydrodynamical simulation codes such as Gadget (Springel et
al. 2001), RAMSES (Teyssier 2002), and AREPO (Springel 2010), it became easier
to run galaxy simulations, and even simulation data are made public (e.g.,
EAGLE, Illustris). Needless to say, the simulation results highly depend on
the input baryon physics, namely nuclear astrophysics, for predicting chemical
abundances. In this review we first discuss how the yields are constrained
with observations using a simple chemical evolution model, and then how
chemical enrichment can be calculated in more sophisticated, chemodynamical
simulations of galaxies.

##  2 Galactic chemical evolution

Galactic chemical evolution (GCE) has been calculated analytically and
numerically since 70’s (e.g., Tinsley 1980; Prantzos et al. 1993; Timmes et
al. 1995; Pagel 1997; Chiappini et al. 1997; Matteucci 2001; Kobayashi et al.
2000, hereafter K00) basically integrating the following equation:

$\frac{d(Z_{i}f_{\rm g})}{dt}=E_{\rm SW}+E_{\rm SNcc}+E_{\rm SNIa}+E_{\rm
NSM}-Z_{i}\psi+Z_{i,{\rm inflow}}R_{\rm inflow}-Z_{i}R_{\rm outflow}$ (1)

where the mass fraction of each element $i$ in gas-phase ($f_{\rm g}$ denotes
the gas fraction, or the gas mass in the system considered with a unit mass)
increases via element ejections from stellar winds ($E_{\rm SW}$), core-
collapse supernovae ($E_{\rm SNcc}$), Type Ia supernovae ($E_{\rm SNIa}$), and
neutron star mergers ($E_{\rm NSM}$). It also decreases by star formation
(with a rate $\psi$), and can be modified by inflow (with a rate $R_{\rm
inflow}$) and outflow (with a rate $R_{\rm outflow}$) of gas in/from the
system considered. It is assumed that the chemical composition of gas is
instantaneously well mixed in the system (called an one-zone model), but the
instantaneous recycling approximation is not adopted nowadays (Matteucci
2021). The model with $R_{\rm inflow}=R_{\rm outflow}=0$ is called a closed-
box model, but is not realistic in any observed galaxies. The initial
conditions are $f_{{\rm g},0}=1$ (a closed system) or $f_{{\rm g},0}=0$ (an
open system) with the chemical composition ($Z_{i,0}$) from the Big Bang
nucleosynthesis ($Z_{i,{\rm BBN}}$). External enrichment is often neglected,
assuming $Z_{i,{\rm inflow}}=Z_{i,{\rm BBN}}$. The first two terms depend only
on nucleosynthesis yields, while the third and fourth terms also depend on
modelling of the progenitor systems, which is uncertain. The last three terms
are galactic terms, and should be determined from galactic dynamics, but are
assumed with analytic formulae in GCE models. The stars formed at a given time
$t$ have the initial metallicity, which is equal to the metallicity of the ISM
from which the stars form, $Z(t)$111Later, integrated metallicity of stellar
populations in a galaxy will be described as $Z_{*}$..

###  2.1 The metal ejection terms in GCE

The ejection terms are given by integrating the contributions with various
initial mass ($m$) and metallicity ($Z$) of progenitor stars. Hence, the first
assumption is the initial mass function (IMF), $\phi$. The IMF is often
assumed to have time/metallicity-invariant mass spectrum normalized to unity
at $m_{\ell}\leq m\leq m_{u}$ as

$\phi(m)\propto m^{-x},$ (2)
$\int_{m_{\ell}}^{m_{u}}m^{-x}dm=\frac{1}{1-x}(m_{u}^{1-x}-m_{\ell}^{1-x})=1.$
(3)

Salpeter (1955) found the slope of the power-law, $x=1.35$, from observation
of nearby stars. See Tinsley (1980) for a different definition of the slope,
$\alpha=x+1$. Kroupa et al. (1993) updated this result finding three slopes at
different mass ranges, which is summarised as $x=-0.7,0.3$, and $1.3$ for
$m/M_{\odot}<\kern-11.99998pt\lower
4.73611pt\hbox{$\sim$}~{}0.08,0.08<\kern-11.99998pt\lower
4.73611pt\hbox{$\sim$}~{}m/M_{\odot}<\kern-11.99998pt\lower
4.73611pt\hbox{$\sim$}~{}0.5$, $0.05<\kern-11.99998pt\lower
4.73611pt\hbox{$\sim$}~{}m/M_{\odot}<\kern-11.99998pt\lower
4.73611pt\hbox{$\sim$}~{}150$, respectively (Kroupa 2008). However, the author
suggested that the slope of the massive end should be $x=1.7$, which has been
used in some GCE works (e.g., Romano et al. 2010). In most of results in this
paper the Kroupa IMF with $x=1.3$ is adopted for a mass range from
$m_{\ell}=0.01M_{\odot}$ to $m_{u}=120M_{\odot}$, while (Chabrier 2003)’s IMF
is often used in hydrodynamical simulations. Note that there are a few claims
for an IMF variation from observations (e.g., van Dokkum & Conroy 2010;
Gunawardhana et al. 2011; Cappellari et al. 2013).

Then the ejection terms for single stars can be calculated as follows222The
metals in Eq.4 and Eq.5 are called ‘unprocessed’ and ‘processed’ metals,
respectively. Our supernova yield table contains ‘processed’ metals only, the
AGB yield table contains net yields, so that the mass loss with star
particle’s chemical composition should also be included with Eq.4..

$E_{\rm SW}=\int_{m_{t}}^{m_{u}}\,(1-w_{m}-p_{z_{i}m,{\rm
II}})\,Z_{i}(t-\tau_{m})\,\psi(t-\tau_{m})\,\phi(m)~{}dm,$ (4) $E_{\rm
SNcc}=\int_{m_{t}}^{m_{u}}\,p_{z_{i}m,{\rm
II}}\,\psi(t-\tau_{m})\,\phi(m)~{}dm.$ (5)

The lower mass limit for integrals is the turning off mass $m_{t}$ at $t$,
which is the mass of the star with the main-sequence lifetime $\tau_{m}=t$.
$\tau_{m}$ is a lifetime of a star with $m$, and also depends on $Z$ of the
star. $w_{m}$ is the remnant mass fraction, which is the mass fraction of a WD
(for stars with initial masses of $m<\kern-11.99998pt\lower
4.73611pt\hbox{$\sim$}~{}8M_{\odot}$), a NS (for $\sim 8-20M_{\odot}$) or a
black hole (BH) (for $>\kern-11.99998pt\lower
4.73611pt\hbox{$\sim$}~{}20M_{\odot}$). $p_{z_{i}m,{\rm II}}$ denotes the
nucleosynthesis yields of core-collapse supernovae, given as a function of
$m$, $Z$, and explosion energy.

The additional ejection terms from binary systems are more complicated.
Suppose that the event rates are given, the enrichment can be calculated as

$E_{\rm Ia}=m_{\rm CO}\,p_{z_{i}m,{\rm Ia}}\,{\cal R}_{t,{\rm Ia}},$ (6)
$E_{\rm NSM}=m_{\rm NSM\,ejecta}\,p_{z_{i}m,{\rm NSM}}\,{\cal R}_{t,{\rm
NSM}}.$ (7)

${\cal R}_{t,{\rm x}}$ are the rates of SNe Ia and NSMs per unit time per unit
stellar mass formed in a population of stars with a coeval chemical
composition and age (simple stellar population, SSP), and called the delay-
time distribution (DTD). For the DTDs, simple analytic formula were also
proposed (Matteucci 2021, and references therein). However, the functions are
significantly different from what are obtained from binary population
synthesis (BPS), ignore metallicity dependences during binary evolution, and
do not includes the effects of supernova kicks for NSMs. De Donder &
Vanbeveren (2004) was the first work that combined BPS to GCE.

All of the matter in the WD ($m_{\rm CO}$) is ejected from SNe Ia (except for
a sub-class called SNe Iax, see Kobayashi et al. 2015); for Ch-mass SNe Ia,
$m_{\rm CO}=1.38M_{\odot}$. On the other hand, only a small fraction of matter
is ejected from a NSM ($m_{\rm NSM\,ejecta}\sim 0.01M_{\odot}$), depending on
the mass ratios, the equation of state of the NSs, and the spin of BHs
(Kobayashi et al. 2022). These are provided together with the nucleosynthesis
yields ($p_{z_{i}m,{\rm Ia}}$ and $p_{z_{i}m,{\rm NSM}}$).

For SNe Ia from the single-degenerate systems, Kobayashi et al. (1998)
proposed another analytic formula based on binary calculation:

${\cal R}_{t,{\rm Ia}}=b~{}\int_{\max[m_{{\rm p},\ell},\,m_{t}]}^{m_{{\rm
p},u}}\,\frac{1}{m}\,\phi(m)~{}dm~{}\int_{\max[m_{{\rm
d},\ell},\,m_{t}]}^{m_{{\rm d},u}}\,\frac{1}{m}\,\psi(t-\tau_{m})\,\phi_{\rm
d}(m)~{}dm.$ (8)

The first integral is for the primary star, and the mass range is for the
stars that can produce $\sim 1M_{\odot}$ of C+O WDs in binaries, which is set
to be $3-8M_{\odot}$ in K20. The second integral is for the secondary star,
and the mass range depends on the optical thick winds from the WD, which is
about $\sim 1M_{\odot}$ and $\sim 3M_{\odot}$ for the red-giant+WD systems and
main-sequence+WD systems, respectively, depending on $Z$.

The metallicity dependence on the SN Ia DTD is very important to reproduce the
observed [$\alpha$/Fe]–[Fe/H] relation in the solar neighbourhood (see
Kobayashi & Nomoto 2009, for the detailed parameters and the resultant DTDs).
These parameters are for Ch-mass explosions. Eq 8 can also be used for SNe Iax
and sub-Ch mass explosion triggered by slow H accretion (see Kobayashi et al.
2015, for the detailed parameters), but not for the double-degenerate systems
(which are likely to cause sub-Ch mass explosions). DTDs of these sub-classes
of SNe Ia have been predicted by various BPS models, but none of these models
can reproduce the observations because the total rate is too low and/or the
typical timescale is too short (see Kobayashi et al. 2022, for more details).

For NSMs, K20 used a metallicity-dependent DTD from a BPS (Mennekens &
Vanbeveren 2014, 2016). Kobayashi et al. (2022) used various BPS models and
also provided new analytic formulae that can reproduce the observed
[Eu/(Fe,O)]–[(Fe,O)/H] relation only with NSMs, without the r-process
associated with core-collapse supernovae. Currently, there is no BPS model
that can explain the observation only with NSMs because the rate is too low
and/or the timescale is too long (see Kobayashi et al. 2022, for more
details). Other binary systems such as novae can also contribute GCE for some
elements or isotopes (e.g., Romano et al. 2019) but the DTDs are very
uncertain (see Kemp et al. 2022, for the nova DTDs from BPS) and the yield
tables are not available.

###  2.2 Nucleosynthesis yields

Nucleosynthesis yields ($p_{z_{i}m,{\rm X}}$) are integrated over stellar
lifetimes of single stars, or delay time distributions for binaries, depending
on $Z$ (Eqs. 4-7). It is extremely important to take account of this
metallicity dependence, though it is often ignored in hydrodynamical
simulations. Because of this, it is also not possible to solve chemical
evolution as a post-process. Chemical enrichment including all elements must
be followed on-the-fly.

For core-collapse supernovae, our yields were originally calculated in K06, 3
models of which are replaced in Kobayashi et al. (2011a) (the identical table
was also used in Nomoto et al. 2013), and a new set with failed supernovae is
used in K20. This is based on the lack of observed progenitors at supernova
locations in the HST data (Smartt 2009), and on the lack of successful
explosion simulations for massive stars (Janka 2012; Burrows & Vartanyan
2021). This resulted in a 20% reduction of both the net and oxygen yields (see
Table 3 of K20). Supernova yields assumed $E_{51}\equiv E/10^{51}{\rm ergs}=1$
of explosion energy as observed for SN1987A, while hypernova yields included
the observed mass dependence of the explosion energy; $E_{51}=$ 10, 10, 20, 30
for 20, 25, 30, $40M_{\odot}$ stars. These supernova and hypernova yield
tables are provided separately, and it is recommended to assume $\epsilon_{\rm
HN}=$50% hypernova fraction for stars with $\geq 20M_{\odot}$ (Kobayashi et
al. 2006). This fraction is expected to decrease at high metallicities due to
smaller angular momentum loss, and was assumed as $\epsilon_{\rm
HN}(Z)=0.5,0.5,0.4,0.01$, and $0.01$ for $Z=0,0.001,0.004,0.02$ in Kobayashi &
Nakasato (2011).

Woosley & Weaver (1995) provided a yield table as a function of mass and
metallicity, which was tested with a GCE model by Timmes et al. (1995), and
has been widely used in chemical evolution studies. It would be useful to note
that stellar mass loss is not included in their presupernova evolution
calculations. Another problem is that their core-collapse supernova models
tend to produce more Fe than observed because of the relatively deep mass cut,
which leads to [$\alpha$/Fe]$\sim 0$ in the ejecta. Therefore, the theoretical
Fe yields were divided by a factor of 2 (Timmes et al. 1995; Romano et al.
2010), which is obviously unphysical. This artificial reduction could mimic a
shallower mass cut, but in that case the yields of other iron peak elements
formed in the same layer as Fe should also be reduced. For the first problem,
Portinari et al. (1998) obtained C+O core masses from stellar evolution models
with mass loss and adopted to the Woosley & Weaver (1995) supernova yields,
which is not physically accurate either. It is also known that Mg production
in the $40M_{\odot}$ model with $E_{51}=1$ is unreasonably small compared with
that in other models. This Mg underproduction problem also exists in Portinari
et al. (1998).

Limongi & Chieffi (2018) also provided a yield table but as a function of
stellar rotation. See Kobayashi (2022a) for the impact of stellar rotation in
GCE. Note that these yields do not reproduce the observed elemental abundances
in the solar neighbourhood, namely iron-peak elements. This is because this
yield set does not have hypernovae, and mixing during explosion is not taken
into account either.

There are also pre-supernova yields in the literature that do not include
explosive nucleosynthesis. These are not useful for GCE since during
explosions iron-peak elements are produced and $\alpha$ element yields are
also largely modified, although they can be used for studying isotopic ratios
of light elements.

The yields for AGB stars were originally calculated in Karakas (2010) and
Karakas & Lugaro (2016), but a new set with the s-process is used in K20 with
optimising the mass of the partial mixing zone. In K20, the narrow mass range
of super-AGB stars is also filled with the yields from Doherty et al. (2014);
stars at the massive end are likely to become ECSNe, as is believed to have
been the case for the progenitor of the Crab Nebula in 1054. At the low-mass
end, off-centre ignition of C flame moves inward but does not reach the
centre, which remains a hybrid C+O+Ne WD. This might become a sub-class of SNe
Ia called SN Iax. The contribution of SN Iax is negligible in the solar
neighbourhood, but can be important for dwarf spheroidal galaxies (Kobayashi
et al. 2015; Cescutti & Kobayashi 2017). There are other yields in the
literature, see Romano et al. (2010) and Romano et al. (2019) for comparison
with GCE.

For SNe Ia, the nucleosynthesis yields of near-Ch and sub-Ch mass models are
newly calculated in Kobayashi et al. (2020b), which used the same 2D code as
in Leung & Nomoto (2018) and Leung & Nomoto (2020) but with more realistic,
solar-scaled initial composition. The initial composition gives significantly
different (Ni, Mn)/Fe ratios, compared with the classical W7 model (Nomoto et
al. 1997) or more recent delayed detonation model (Seitenzahl et al. 2013).
When constraining the progenitors of SNe Ia from the observed Mn abundances in
the solar neighborhood, it is important to use the latest yields of SNe Ia.
See Blondin et al. (2022) for a complete comparison including many other
yields.

The r-process yields are taken from the following references in K20:
$8.8M_{\odot}$ model in Wanajo et al. (2013, 2D) for electron capture
supernovae, $1.2-2.4M_{\odot}$ models in Wanajo (2013) for $\nu$-driven winds,
$1.3M_{\odot}$+$1.3M_{\odot}$ model in Wanajo et al. (2014, 3D-GR) for neutron
star mergers, and $25M_{\odot}$ “b11tw1.00” model in Nishimura et al. (2015,
axisymmetric MHD) for magnetorotational supernovae. See Kobayashi et al.
(2022) for the GCE models including the NSM yields from Just et al. (2015)
with various masses of compact objects. These yields are very uncertain
depending on 1) initial conditions (e.g., $M$, $Z$, rotation, magnetic
fields), 2) the quality of the base hydrodynamical simulations (3 dimensional
general relativity or not), 3) neutrino physics, and 4) nuclear physics (e.g.,
reaction rates, fission modelling).

Figure 2:  Star formation histories (panel a), age-metallicity relations
(panel b), and metallicity distribution functions (panel c) for the solar
neighborhood (blue solid lines), halo (green short-dashed lines), halo with
stronger outflow (light-blue dotted lines), bulge (red long-dashed lines),
bulge with outflow (olive dot-short-dashed lines) and thick disk (magenta dot-
long-dashed lines). Figure is taken from Kobayashi et al. (2020a), see the
reference for the observational data sources and the model details.  Figure 3:
The same as Fig. 2 but for dwarf spheroidal galaxies, Fornax (red solid
lines), Sculptor (green long-dashed lines), Sextans (blue short-dashed lines),
and Carina (magenta dotted lines), and ultra-faint galaxies (cyan dot-dashed
lines). Figure is updated from Kobayashi et al. (2020b), see the reference for
the observational data sources and the model details.

###  2.3 The galactic terms in GCE

Star formation rate is usually modelled as $\psi\\!=\\!\frac{1}{\tau_{\rm
s}}f_{\rm g}^{\,n}$ but often $n=1$ is adopted, and thus it is proportional to
the gas fraction. ${\tau_{\rm s}}$ denotes a timescale of star formation, and
$1/{\tau_{\rm s}}$ gives the efficiency of star formation. The inflow rate is
often assumed to be exponential as, $R_{\rm
inflow}\\!\\!=\\!\frac{1}{\tau_{\rm i}}\exp\frac{-t}{\tau_{\rm i}}$, but
occasionally as $R_{\rm inflow}\\!\\!=\\!\frac{t}{\tau_{\rm
i}^{2}}\exp\frac{-t}{\tau_{\rm i}}$ in Pagel (1997) and in the solar
neighbourhood model of K20. The outflow rate is assumed to be proportional to
the star formation rate as $R_{\rm outflow}\\!\\!=\\!\frac{1}{\tau_{\rm
o}}f_{\rm g}^{\,n}\propto\psi$, which is reasonable if the outflow is driven
by supernova feedback. In addition, star formation is quenched ($\psi=0$) at a
given epoch $t_{\rm w}$, which corresponds to a galactic wind driven by the
feedback from active galactic nuclei (AGN).

The timescales are determined to match the observed metallicity distribution
function (MDF) of stars in each system modelled. Figures 2 and 3 show the
assumed star formation history and the resultant metallicity evolution and
metallicity distribution function. The parameter sets that have very similar
MDFs give almost identical tracks of elemental abundance ratios (Fig. A1 of
Kobayashi et al. 2020b). This means that, for the given MDF, elemental
abundance tracks do not depend so much on the star formation history.
Therefore GCE can be used to constrain nuclear astrophysics, if the MDF is
known, in Galactic archaeology. Without knowing the MDFs, star formation
histories are unconstrained, and elemental abundance tracks can vary depending
on the star formation histories. Provided that nuclear astrophysics is known,
elemental abundance measurements can be used to constrain the star formation
histories through GCE, which is the case for extra-galactic archaeology.

###  2.4 The origin of elements

Figure 4: The time evolution (in Gyr) of the origin of elements in the
periodic table: Big Bang nucleosynthesis (black), AGB stars (green), core-
collapse supernovae including SNe II, HNe, ECSNe, and MRSNe (blue), SNe Ia
(red), and NSMs (magenta). The amounts returned via stellar mass loss are also
included for AGB stars and core-collapse supernovae depending on the
progenitor mass. The dotted lines indicate the observed solar values. Figure
is taken from Kobayashi et al. (2020a).

Using the K20 GCE model for the solar neighbourhood, we summarize the origin
of elements in the form of a periodic table. In each box of Figure 4, the
contribution from each chemical enrichment source is plotted as a function of
time: Big Bang nucleosynthesis (black), AGB stars (green), core-collapse
supernovae including SNe II, HNe, ECSNe, and MRSNe (blue), SNe Ia (red), and
NSMs (magenta). It is important to note that the amounts returned via stellar
mass loss are also included for AGB stars and core-collapse supernovae
depending on the progenitor star mass (Eq. 4). The x-axis of each box shows
time from $t=0$ (Big Bang) to $13.8$ Gyrs, while the y-axis shows the linear
abundance relative to the Sun, $X/X_{\odot}$. The dotted lines indicate the
observed solar values, i.e., $X/X_{\odot}=1$ and $4.6$ Gyr for the age of the
Sun; the solar abundances are taken from Asplund et al. (2009), except for
$A_{\odot}({\rm O})=8.76,A_{\odot}({\rm Th})=0.22$, and $A_{\odot}({\rm
U})=-0.02$ (§2.2 of K20 for the details). The adopted star formation history
is similar to the observed cosmic star formation rate history, and thus this
figure can also be interpreted as the origin of elements in the universe,
which can be summarized as follows:

  * •

H and most of He are produced in Big Bang nucleosynthesis. As noted, the green
and blue areas also include the amounts returned to the ISM via stellar mass
loss in addition to He newly synthesized in stars. After tiny production in
Big Bang nucleosynthesis, Be and B are supposed to be produced by cosmic rays
(Prantzos et al. 1993), which are not included in the K20 model.

  * •

The Li model is very uncertain because the initial abundance and
nucleosynthesis yields are uncertain. Li is supposed to be produced also by
cosmic rays and novae, which are not included in the K20 model. The observed
Li abundances show an increasing trend from very low metallicities to the
solar metallicity, which could be explained by cosmic rays. Then the
observation shows a decreasing trend from the solar metallicities to the
super-solar metallicities, which might be caused by the reduction of the nova
rate (Grisoni et al. 2019); this is also shown in theoretical calculation with
binary population synthesis (Kemp et al. 2022), where the nova rate becomes
higher due to smaller stellar radii and higher remnant masses at low
metallicities.

  * •

49% of C, 51% of F, and 74% of N are produced by AGB stars333In extra-galactic
studies, the N production from AGB stars is referred as “secondary”, which is
confusing, and is a primary production from freshly synthesized 12C. For
massive stars, the N yield depends on the metallicity of progenitor stars for
secondary production, and can also be enhanced by stellar rotation for primary
production (Kobayashi et al. 2011a). (at $t=9.2$ Gyr). For the elements from
Ne to Ge, the newly synthesized amounts are very small for AGB stars, and the
small green areas are mostly for mass loss.

  * •

$\alpha$ elements (O, Ne, Mg, Si, S, Ar, and Ca) are mainly produced by core-
collapse supernovae, but 22% of Si, 29% of S, 34% of Ar, and 39% of Ca come
from SNe Ia. These fractions would become higher with sub-Ch-mass SNe Ia
(Kobayashi et al. 2020b) instead of 100% Ch-mass SNe Ia adopted in the K20
model.

  * •

A large fraction of Cr, Mn, Fe, and Ni are produced by SNe Ia. In classical
works, most of Fe was thought to be produced by SNe Ia, but the fraction is
only 60% in the K20 model, and the rest is mainly produced by HNe. The
inclusion of HNe is very important as it changes the cooling and star
formation histories of the universe significantly (Kobayashi et al. 2007). Co,
Cu, Zn, Ga, and Ge are largely produced by HNe. In the K20 model, 50% of stars
at $\geq 20M_{\odot}$ are assumed to explode as hypernovae, and the rest of
stars at $>30M_{\odot}$ become failed supernovae.

  * •

Among neutron-capture elements, as predicted from nucleosynthesis yields, AGB
stars are the main enrichment source for the s-process elements at the second
(Ba) and third (Pb) peaks.

  * •

32% of Sr, 22% of Y, and 44% of Zr can be produced from ECSNe, which are
included in the blue areas, even with the adopted conservative mass ranges; we
take the metallicity-dependent mass ranges from the theoretical calculation of
super-AGB stars (Doherty et al. 2015). Combined with the contributions from
AGB stars, it is possible to perfectly reproduce the observed trends, and no
extra light element primary process (LEPP) is needed (but see Prantzos et al.
2018). The inclusion of $\nu$-driven winds in GCE simulation results in a
strong overproduction of the abundances of the elements from Sr to Sn with
respect to the observations.

  * •

For the heavier neutron-capture elements, contributions from both NS-NS/NS-BH
mergers and MRSNe are necessary, and the latter is included in the blue areas.
Note again that the green areas includes the mass-loss component, i.e., not
newly produced but recycled. Note that Tc and Pm are radioactive.

Since the Sun is slightly more metal-rich than the other stars in the solar
neighborhood (see Fig. 2 of K20), the fiducial model in K20 goes through
[O/Fe]$=$[Fe/H]$=0$ slightly later compared with the Sun’s age. Thus, a
slightly faster star formation timescale ($\tau_{\rm s}=4$ Gyr instead of 4.7
Gyr) is adopted in this figure. The evolutionary tracks of [X/Fe] are almost
identical. In this model, the O and Fe abundances go though the cross of the
dotted lines, meaning [O/Fe] $=$ [Fe/H] $=0$ at 4.6 Gyr ago. This is also the
case for some important elements including N, Ca, Cr, Mn, Ni, Zn, Eu, and Th.
The remaining problems can be summarised as follows:

  * •

The contribution from rotating massive stars is not included in the K20 model,
which can probably explain in the underproduction of C and F (Kobayashi
2022a). A binary effect during stellar evolution may also increase C. These
elements could be enhanced by AGB stars as well, but the observed high F
abundance in a distant galaxy strongly supported rapid production of F from
Wolf-Rayet stars (Franco et al. 2021, at $z=4.4$ observed with ALMA).

  * •

Mg is slightly under-produced in the model, although at low metallicity the
model [Mg/Fe] is higher than observed (see Fig. 7). This may be due to a NLTE
effect (Lind et al. 2022), or due to a binary effect. Massive stars in
binaries tend to have a smaller CO core with a higher C/O ratios (Brown et al.
2001), which results in a higher Mg/O ratio. Observed Mg/O ratios suggest that
this binary effect should not be important at low metallicities
($<\kern-11.99998pt\lower 4.73611pt\hbox{$\sim$}~{}0.1Z_{\odot}$).

  * •

The underproduction of the elements around Ti is a long-standing problem since
Kobayashi et al. (2006). It was shown that these elemental abundances can be
enhanced by multi-dimensional effects (Sneden et al. 2016; see also K15 model
in K20). This is due to the lack of 3D simulations that can capture the
production of these (less abundant) elements. 2D nucleosynthesis calculation
showed the enhancement of these elements (Maeda & Nomoto 2003; Tominaga 2009).

  * •

The s-process elements are slightly overproduced even with the updated
s-process yields. Notably, Ag is over-produced by a factor of $6$, while Au is
under-produced by a factor of $5$. U is also over-produced. These problems may
require revising nuclear physics modelling (see Kobayashi 2022a, for more
discussion).

Figure 5: The [O/Fe]–[Fe/H] relations in the solar neighborhood for the models
with Ch-mass SNe Ia only (solid line), 75% Ch plus 25% sub-Ch-mass SNe Ia
(dashed line), and sub-Ch-mass SNe Ia only (dotted line). The observational
data (filled circles) are high-resolution non-local thermodynamic equilibrium
(NLTE) abundances from Zhao et al. (2016). Figure is taken from Kobayashi et
al. (2020b) with modification.

###  2.5 The [X/Fe]–[Fe/H] relations

Although the slopes in Figure 4 show the different enrichment timescales of
each element, it is easier to see the differences in the [X/Fe]–[Fe/H]
diagrams444There were attempts to plot against a different element (e.g., O)
to avoid the uncertainty of Fe yields (e.g., Cayrel et al. 2004), but it
became rather hard to understand the plots.. Figure 5 shows the
[$\alpha$/Fe]–[Fe/H] relation, which is probably the most important diagram in
GCE. O is one of the $\alpha$ elements. At the beginning of the universe, the
first stars (Population III stars) form and die, whose properties such as mass
and rotation are uncertain and have been studied using the abundance patterns
of the second generation, extremely metal-poor (EMP) stars. Secondly, core-
collapse supernovae occur, and their yields are imprinted in the abundance
patterns of Population II stars in the Galactic halo. The [$\alpha$/Fe] ratio
is high and stays roughly constant555There was a debated if the [O/Fe] ratio
increases toward the lowest metallicity or not. UV OH line showed such an
increase (e.g, Israelian et al. 1998) reaching [O/Fe] $\sim 1$ at [Fe/H]
$\sim-3$, which was later found to be due to 3D effects. with a small scatter.
This plateau value does not depend on the star formation history but does on
the IMF. Finally, SNe Ia occur, which produce more Fe than $\alpha$ elements,
and thus the [$\alpha$/Fe] ratio decreases toward higher metallicities; this
decreasing trend is seen for the Population I stars in the Galactic disk.

Figure 6: The [O/Fe]–[Fe/H] relations in the Galactic bulge (red long-dashed
line), thick disk (green short-dashed line), thin disk (blue solid line),
Fornax (magenta dot-short-dashed line), Sculptor (cyan dot-long-dashed line),
and ultra faint dwarf galaxies (olive dotted line). The star formation history
of each system is constrained from its MDF in Figs. 2 and 3. See Fig. 5 for
the observational data source for thin (navy open circles) and thick (maroon
filled circles) disk stars.

The contribution from SNe Ia depends on the progenitor binary systems, namely,
the mass of the progenitor WDs. Sub-Ch mass explosions produce less Mn and Ni,
and more Si, S, and Ar than near-Ch mass explosions (Seitenzahl et al. 2013;
Kobayashi et al. 2020b). Figure 5 shows the [O/Fe]–[Fe/H] relations with
varying the fraction of sub-Ch-mass SNe Ia. Including up to 25% sub-Ch mass
contribution to the GCE (dashed line) gives a similar relation as the K20
model (solid line), while the model with 100% sub-Ch-mass SNe Ia (dotted line)
gives too low an [O/Fe] ratio compared with the observational data. For Ch-
mass SNe Ia, the progenitor model is based on the single-degenerate scenario
with the metallicity effect due to optically thick winds from WDs (Kobayashi
et al. 1998). For sub-Ch-mass SNe Ia, the observed delay-time distribution is
used since the progenitors are the combination of mergers of two WDs in double
degenerate systems and low accretion in single degenerate systems; Kobayashi
et al. (2015)’s formula are for those in single degenerate systems only.

Because of this [$\alpha$/Fe]–[Fe/H] relation, high-$\alpha$ and low-$\alpha$
are often used as a proxy of old and young ages of stars in galaxies,
respectively. Note that, however, that this relation is not linear but is a
plateau and decreasing trend (called a ‘knee’). The location of the ‘knee’
depends on the star formation timescale, with a higher metallicity for faster
star formation history (Matteucci & Brocato 1990), e.g., [Fe/H]
$\sim-0.5,-0.8,-1,-2$ for the Galactic bulge, thick disk, thin disk, and
satellite galaxies. Figure 6 shows the GCE model results for various
environments. The bulge (with outflow), thick disk, and solar neighborhood
(thin disk) models are taken from K20 with 100% Ch-mass SNe Ia, with the star
formation histories in Fig. 2. The models for dwarf spheroidal galaxies are
taken from Kobayashi et al. (2020b) adding 100% sub-Ch mass SNe Ia on top of
the 100% of Ch mass SNe Ia, with the star formation histories in Fig. 3.
Recently, similar relations are also shown for M31 using planetary nebulae.
Since Fe is not available, Ar is used as a significant fraction of Ar is
produced by SNe Ia. The thin and thick disk dichotomy seems not exactly the
same, probably due to M31 experiencing an accretion of a relatively massive,
gas-rich satellite galaxy (Arnaboldi et al. 2022).

There are GCE models that try to explain both thin and thick disk observations
with two infalls (Chiappini et al. 1997; Spitoni et al. 2019). However, more
realistic, chemodynamical simulations show that a significant number of thick
disk stars are formed in satellite galaxies before they accrete onto the disk,
which can be better approximated with two independent GCE models in this
figure (also in Grisoni et al. 2017).

Figure 7:  Evolution of the elemental abundances [X/Fe] from C to Zn against
[Fe/H] for the models in the solar neighborhood, with only supernovae (without
AGB and super-AGB stars, blue short-dashed line), with AGB without super-AGB
stars (green dotted lines), with AGB and super-AGB stars (red solid line,
fiducial model), with ECSNe (magenta long-dashed lines), and with SNe Iax
(cyan dot-dashed lines). Figure is taken from Kobayashi et al. (2020a) and see
the reference for the observational data sources.

Figure 7 shows the evolution of elemental abundance ratios [X/Fe] against
[Fe/H] from C to Zn in the solar neighborhood. All $\alpha$ elements (O, Mg,
Si, S, and Ca) show the [$\alpha$/Fe]–[Fe/H] relations, i.e., the plateau and
decreasing trend from [Fe/H] $\sim-1$ to higher metallicities. The odd-$Z$
elements (Na, Al, and Cu) show an increasing trends toward higher
metallicities due to the metallicity dependence of the core-collapse supernova
yields. See K20 for more detailed discussion on the evolutionary trends of
each element, and Kobayashi (2022b) for the discussion with more recent
observations for Cu and Zn. In the following we focus the role of each
enrichment source.

  * •

The contribution to GCE from AGB stars (green dotted lines in Fig. 7) can be
seen mainly for C and N, and only slightly for Na, compared with the model
that includes supernovae only (blue dashed lines). Hence it seems not possible
to explain the O–Na anti-correlation observed in globular cluster stars (e.g.,
Kraft et al. 1997) with AGB stars (but see Ventura & D’Antona 2009). Although
AGB stars produce significant amounts of Mg isotopes, their inclusion does not
affect the [Mg/Fe]–[Fe/H] relation.

  * •

The contribution from super-AGB stars (red solid lines) is very small; with
super-AGB stars, C abundances slightly decrease, while N abundances slightly
increase. It would be very difficult to put a constraint on super-AGB stars
from the average evolutionary trends of elemental abundance ratios, but it
might be possible to see some signatures of super-AGB stars in the scatters of
elemental abundance ratios.

  * •

With ECSNe (magenta long-dashed lines), Ni, Cu and Zn are slightly increased.
These yields are in reasonable agreement with the high Ni/Fe ratio in the Crab
Nebula (Nomoto 1987; Wanajo et al. 2009).

  * •

No difference is seen with/without SNe Iax (cyan dot-dashed lines) in the
solar neighborhood because the progenitors are assumed to be hybrid WDs, which
has a narrow mass range in the adopted super-AGB calculations Doherty et al.
(2015, $\Delta M\sim 0.1M_{\odot}$). This mass range depends on convective
overshooting, mass-loss, and reaction rates. Even with the wider mass range in
Kobayashi et al. (2015, $\Delta M\sim 1M_{\odot}$), however, the SN Iax
contribution is negligible in the solar neighborhood, but it can be important
at lower metallicities such as in dwarf spheroidal galaxies with stochastic
chemical enrichment (Cescutti & Kobayashi 2017).

Figure 8: The [X/Fe]–[Fe/H] relations for neutron capture elements, comparing
to the models in the solar neighborhood, with s-process from AGB stars only
(blue long-dashed lines), with s-process and ECSNe (light-blue short-dashed
lines), with s-process, ECSNe, and $\nu$-driven winds (green dotted-long-
dashed lines), with s-process, ECSNe, and NS-NS mergers (olive dotted lines),
with s-process, ECSNe, and NS-NS/NS-BH mergers (orange dotted-short-dashed
lines), with s-process, ECSNe, NS-NS/NS-BH mergers, and MRSNe (red solid
lines). Figure is updated from Kobayashi et al. (2020a) with more
observational data.

Figure 8 shows the evolutions of neutron-capture elements as [X/Fe]–[Fe/H]
relations.

  * •

As known, AGB stars can produce the first (Sr, Y, Zr), second (Ba), and third
(Pb) peak s-process elements, but no heaver elements (navy long-dashed lines).
The second-peak elements are under-produced around [Fe/H] $\sim-3$, which is
eased with chemodynamical simulations, and can also be reproduced better with
rotating massive stars. However, the model with rotating massive stars results
in over-production of the first peak elements (Kobayashi 2022a).

  * •

It is surprising that ECSNe from a narrow mass range ($\Delta M\sim
0.15-0.2M_{\odot}$) can produce enough of the first peak elements; with the
combination of AGB stars, it is possible to reproduce the observational data
very well (cyan short-dashed lines). This means that no other light element
primary process (LEPP), such as rotating massive stars, is required (but see
Cescutti & Chiappini 2014). The elements from Mo to Ag seem to be
overproduced, which could be tested with the UV spectrograph proposed for the
VLT, CUBES.

  * •

Additional production from $\nu$-driven winds leads to further over-production
of these elements in the model (green dot-long-dashed lines), but this should
be studied with more self-consistent calculations of supernova explosions.

  * •

Neutron star mergers can produce lanthanides and actinides, but not enough
(olive dotted lines); the rate is too low and the timescale is too long,
according to binary population synthesis. This is not improved enough even if
NS-BH mergers are included (orange dot-short-dashed lines). It is very
unlikely this problem will be solved even if the mass-dependence of the
nucleosynthesis yields are included, unless the BH spins are unexpectedly high
(Kobayashi et al. 2022).

  * •

In the GCE model with MRSNe (magenta solid lines), it is possible to reproduce
a plateau at low metallicities for Eu, Pt, and Th, relative to Fe. However,
even with including both MRSNe and neutron star mergers, the predicted Au
abundance is more than ten times lower than observed. This underproduction is
seen not only for the solar abundance but also for low metallicity stars
although the observational data are very limited. UV spectroscopy with HST, or
NASA’s future LUVOIR, is needed for investigating this problem further.

In conclusion, an r-process associated with core-collapse supernovae, such as
MRSNe, is required. The same conclusion is obtained with other GCE models and
more sophisticated chemodynamical simulations (e.g., Haynes & Kobayashi 2019;
van de Voort et al. 2020), as well as from the observational constraints of
radioactive nuclei in the solar system (Wallner et al. 2021). See more
discussion in the next section.

It seems not to be easy to solve the missing gold problem with astrophysics;
only Au yields should be increased since Pt is already in good agreement with
the current model, and Ag is rather overproduced in the current model.
However, there are uncertainties in nuclear physics, namely in some nuclear
reaction rates and in the modelling of fission, which might be able to
increase Au yields only, without increasing Pt or Ag. It may be hard to
predict the only one stable isotope of 197Au, while Pt has several stable
isotopes. The predicted Th and U abundances are after the long-term decay, to
be compared with observations of metal-poor stars, and the current model does
not reproduce the Th/U ratio either.

###  2.6 The First Galaxies

Figure 9:  Evolution of CNO (number) abundances in GCE models for the Milky
Way (red solid and blue long-dashed lines) and for a submillimetre galaxy
(green short-dashed and magenta dotted lines), which are calculated with the
original K20 yields (red solid and green short-dashed lines) or the yields
with Wolf-Rayet stars (blue long-dashed and magenta dotted lines; see text for
the model details).

The Atacama Large Millimeter/submillimeter Array (ALMA) has opened a new
window to study elemental abundances and isotopic ratios of light elements at
very high redshifts (Zhang et al. 2018; Franco et al. 2021), which provide us
an independent constraint on stellar nucleosynthesis. In the early universe,
stellar rotation may be important at low metallicities because weaker stellar
winds result in smaller angular momentum loss.

The impact of stellar rotation on stellar evolution and nucleosynthesis has
been studied to explain Wolf-Rayet stars (e.g., Meynet & Maeder 2002; Hirschi
2007; Limongi & Chieffi 2018), and the observed N and 13C abundances at low
metallicities have been used to constrain the effect (Chiappini et al. 2006;
see also Fig. 13 of Kobayashi et al. 2011a)666However, the N/O–O/H relation
can be reproduced by inhomogeneous enrichment in chemodynamical simulations
(Fig. 26).. As discussed in the previous section, the importance of jet-like
supernova explosion for iron-peak elements (Co and Zn) and r-process elements
(via MRSNe) also indicates the importance of stellar rotation, although how to
transport the angular momentum in stellar envelopes to stellar cores is
uncertain.

Figure 9 shows the evolution of CNO abundances in GCE models for two
representative galaxies with different star formation histories: the solar
neighborhood in the Milky Way, which has continues star formation, and a
submillimetre galaxy (SMG), which underwent a star burst. The red solid lines
are the same as in Fig. 8 and Fig. 10 of K20. As discussed, $\sim 50$% C and
$\sim 75$% N comes from low- and intermediate-mass AGB stars, respectively, at
later times ([O/H] $\sim-1$ for C). The rest comes from massive stars, either
released before supernova explosions by stellar winds, or ejected at supernova
explosions. The N yield from massive stars increases while the C yield
decreases at higher metallicities. These two effects cause the increase of N/O
and C/O, and the decrease of C/N, toward high metallicities.

With rotation, both C and N yields are increased due to the interplay between
the core He-burning and the H-burning shell, triggered by the rotation-induced
instabilities (Limongi & Chieffi 2018). Since Limongi & Chieffi (2018)’s
yields do not reproduce the observations of iron-peak elements because of
their simple description of supernova explosions. Therefore, the C and N
yields from Limongi & Chieffi (2018) are combined with the K20 yields in the
models with Wolf-Rayet stars. High abundance of CNO enhances the production of
F in He convective shell, which can explain the high F abundance observed in a
submillimetre galaxy at redshift $z=4.4$ (Franco et al. 2021).

The evolution of isotopic ratios is also largely affected by stellar rotation.
Without rotation, 13C and 25,26Mg are produced from AGB stars (17O might be
overproduced in our AGB models), while other minor isotopes are more produced
from metal-rich massive stars, and thus the radios between major and minor
isotopes (e.g., 12C/13C, 16O/17,18O) generally decrease as a function of
time/metallicity (Figs. 17-19 of Kobayashi et al. (2011a) and Fig. 31 of K20;
see also Romano et al. 2019).

With rotation, 13C and 17,18O are enhanced. The evolution of isotopic ratios
and comparison to observations can be found in Fig. 7 of Kobayashi (2022a),
which gives consistent results with Romano et al. (2019)’s GCE models. 12C/13C
ratio became too high, while it is possible to explain the low 13C/18O ratio
observed in submillimetre galaxies at redshift $z\sim 2-3$ (Zhang et al. 2018)
without changing IMF.

Stellar rotation also cause the weak-s process that forms s-process elements
from the existing seeds at a much shorter timescale than AGB stars. This may
help solving the underproduction of Ba at [Fe/H] $<\kern-11.99998pt\lower
4.73611pt\hbox{$\sim$}~{}-1$ (Fig. 8), but cause an overproduction of light
s-process elements such as Sr, Y, and Zr. The impact of stellar rotation on
nucleosynthesis is not well understood yet.

###  2.7 The first chemical enrichment

Figure 10:  Elemental abundance pattern of the metal-poor C-rich DLA (filled
circles) and peculiar DLA (open circles). The red solid and green short-dashed
lines show the nucleosynthesis yields of faint core-collapse supernovae from
$25M_{\odot}$ stars with mixing fallback. The blue dotted line is for pair-
instability supernovae from $170M_{\odot}$ stars. Figure is taken from
Kobayashi et al. (2011b) with modification.

Extremely metal-poor (EMP) stars have been an extremely useful relic in the
Galactic archaeology. At the beginning of galaxy formation, stars form from a
gas cloud that was enriched only by one or very small number of supernovae
(Audouze & Silk 1995), and hence the elemental abundances of EMP stars can
offer observational evidence of particular types of supernovae in the past.
The expectation was that the first stars were so massive that they exploded as
a pair-instability supernova (PISN, Barkat et al. 1967), which causes high
[(Si,S)/O] ratios (e.g., Nomoto et al. 2013).

However, after a half century of surveys, no star has been found with an
elemental abundance pattern fitted by a PISN in the Milky Way. Instead, it is
found that quite a large fraction of massive stars become faint supernovae
(Umeda & Nomoto 2003) that give a high C/Fe ratio leaving a relatively large
black hole ($\sim 5M_{\odot}$, Nomoto et al. 2013). If there are C-enhanced
low-$\alpha$ stars, that would indicate black hole formation even from
$10-20M_{\odot}$ Pop III stars (Kobayashi et al. 2014). From faint supernovae,
the ejected iron mass is very small but the explosion energy can be high;
faint hypernovae also show Zn enhancement (Ishigaki et al. 2018).

Similar C enhancement might also be found for damped Ly-$\alpha$ system (DLA,
Kobayashi et al. 2011b). Figure 10 shows the observed abundance pattern of a
DLA at $z=2.34$ with [Fe/H] $\simeq-3$, comparing to the nucleosynthesis
yields of a faint supernova (red solid line), a faint hypernova (green dashed
line), and a PISN (blue dotted line). In recent surveys, no DLA is found with
an elemental abundance pattern consistent with a PISN either.

Figure 11: The elemental abundance of a extremely metal-poor star SMSS
J200322.54-114203.3, which has [Fe/H]$=-3.5$, comparing with mono-enrichment
from a magneto-rotational hypernovae (blue dashed line and squares) and multi-
enrichment including a neutron star merger (orange dot-dashed line and
triangles). The lower panel shows the differences, i.e., the observed values
minus model values. Figure is taken from Yong et al. (2021).

A small number of EMP stars show a relatively low $\alpha$ abundance, which
does not necessarily mean that the ISM from which the star formed was already
enriched by SNe Ia. The reasons that could cause low $\alpha$/Fe ratios were
summarised in Kobayashi et al. (2014): (1) SNe Ia, (2) less-massive core-
collapse supernovae ($<\kern-11.99998pt\lower
4.73611pt\hbox{$\sim$}~{}20M_{\odot}$), which become more important with a low
star formation rate, (3) hypernovae, although the majority of hypernovae are
expected to give normal [$\alpha$/Fe] ratios ($\sim 0.4$), and (4) PISNe,
which could be very important in the early universe. Therefore, the
[$\alpha$/Fe] ratio is not a perfect clock. It is necessary to also use other
elemental abundances (namely, Mn and neutron-capture elements) or isotopic
ratios, with higher resolution ($>40,000$) multi-object spectroscopy on 8m
class telescopes (e.g., cancelled WFMOS or planned MSE, proposed HRMOS for
VLT).

EMP stars with r-process enhancement (called rII stars), namely those with
enhancement of the heaviest elements (actinide-boost stars), can offer
confirmation of an r-process site. Recently, from the SkyMapper survey, Yong
et al. (2021) found such a star with [Fe/H] $=-3.5$, the lowest metallicty
known for rII stars. The abundance pattern is shown Figure 11, which shows a
clear detection of U and Th, the universal r-process pattern, normal $\alpha$
enhancement, but unexpectedly showing high Zn and N abundances
([Zn/Fe]$=0.72$, [N/Fe]$=1.07$). This abundance pattern cannot be explained
with a model with a neutron star merger (orange dot-dashed line), but instead
strongly supports a magnetorotational hypernova from a
$>\kern-11.99998pt\lower 4.73611pt\hbox{$\sim$}~{}25M_{\odot}$ Pop III star
(blue dashed line).

There are a few magneto-hydrodynamical simulations and post-process
nucleosynthesis that successfully showed enough neutron-rich ejecta to produce
the 3rd peak r-process elements as in the Sun (Winteler et al. 2012; Mösta et
al. 2018). The predicted iron mass is rather small (Nishimura et al. 2015;
Reichert et al. 2021), and the astronomical object may be faint. Yong et al.
(2021)’s proposal is Fe-producing, magnetorotational hypernova, which is more
luminous. The r-process elements could also be formed in the accretion disk of
a collapsar (Siegel et al. 2019), but the inter ejecta that contains Fe and Zn
must be ejected at the explosion as for a hypernova. As it is luminous it
might be possible to detect absorption features of Au and Pt in the spectra of
the transient in future.

##  3 Chemodynamical evolution of galaxies

Thanks to the development of super computers and numerical techniques, it is
now possible to simulate the formation and evolution of galaxies from
cosmological initial conditions, i.e., initial perturbation based on $\Lambda$
cold dark matter (CDM) cosmology. Star formation, gas inflow, and outflow in
Eq.(1) are not assumed but are, in principle, given by dynamics. The baryon
cycle is schematically shown in Figure 12. Due to the limited resolution, star
forming regions, supernova ejecta, and active galactic nuclei (AGN) cannot be
resolved in galaxy simulations, and thus it is necessary to model star
formation and the subsequent effects – feedback – introducing a few
parameters. Fortunately, there are many observational constraints, from which
it is usually possible to choose a best set of parameters. To ensure this, it
is also necessary to run the simulation until the present-day, $z=0$, and
reproduce a number of observed relations at various redshifts.

Figure 12:  Schematic view of chemodynamical evolution of galaxies.

Although hydrodynamics can be calculated with publicly available codes such as
Gadget, RAMSES, and AREPO, modelling of baryon physics is the key and is
different depending on the simulation teams/runs. These are various
simulations of a cosmological box with periodic boundary conditions, and
contain galaxies with a wide mass range (e.g., $10^{9-12}M_{\odot}$ in stellar
mass) at $z=0$. In order to study massive galaxies it is necessary to increase
the size of the simulation box (e.g., $>\kern-11.99998pt\lower
4.73611pt\hbox{$\sim$}~{}100$ Mpc), while in order to study internal
structures it is necessary to improve the spacial resolution (e.g.,
$<\kern-11.99998pt\lower 4.73611pt\hbox{$\sim$}~{}0.5$ kpc). The box size and
resolution are chosen depending on available computational resources. On the
other hand, zoom-in techniques allow us to increase the resolution focusing a
particular galaxy, but this also requires tuning the parameters with the same
resolution, comparing to a number of observations in the Milky Way. Needless
to say, elemental abundances are the most informative. There are a few zoom-in
simulations for Milky Way-type galaxies (e.g., Brook et al. 2012; Few et al.
2014; Grand et al. 2017; Buck et al. 2020; Font et al. 2020), but in most of
the cases the input nuclear astrophysics are too simple to make use of
elemental abundances.

Figure 13: (left) Cooling functions without photoionization. Included
processes are collisional excitation (long-dashed lines) of H0 (red) and He+
(cyan), collisional ionization (short-dashed lines) of H0 (red), He0 (green),
and He+ (cyan), standard recombination (short-dash dotted lines) of H+
(yellow), He+ (cyan), and He++ (blue), dielectric recombination (long-dashed
dotted line) of He+ (cyan), and free-free emission (dotted line). (right)
Metallicity dependence of the cooling function. The solid, long-dashed, short
dashed, dotted, short-dash dotted, and long-dash dotted lines correspond to
[Fe/H] $=0,-1,-2,-3,-4$, and $-5$, respectively. Figures are taken from
Kobayashi (2002).

###  3.1 Modelling of baryon physics

For the baryon physics, the first process to calculate is radiative cooling,
and photo-heating usually by a uniform and evolving UV background radiation
(Haardt & Madau 1996). It would be ideal to simulate radiative feedback self-
consistently but this is computationally very expensive and the available
simulations are only for high-redshifts. The cooling function depends on the
chemical composition of gas, and a metallicity-dependent cooling function
computed with the MAPPINGS III software (Sutherland & Dopita 1993; Kewley et
al. 2019) is used in our code, assuming the [$\alpha$/Fe]–[Fe/H] relation of
the solar neighborhood (Figure 13). Element-dependent cooling functions are
also available (Wiersma et al. 2009). Usually, the simulation code finds
candidate gas particles (or mesh) that can form stars in a given timestep, and
forms star particles with some certain mass.

Our star formation criteria are:

$\mbox{converging flow, }(\nabla\cdot\mbox{\boldmath$v$})_{i}<0,$ (9)
$\mbox{rapid cooling, }t_{\rm cool}<t_{\rm dyn},$ (10) $\mbox{Jeans unstable
gas, }t_{\rm dyn}<t_{\rm sound}.$ (11)

This still needs a parameter to define the star formation timescale relative
to the dynamical timescale, as

$t_{\rm sf}\equiv\frac{1}{c_{*}}t_{\rm dyn},$ (12)

where $c_{*}$ is a star formation timescale parameter. The parameter $c_{*}$
basically determines when to form stars following cosmological accretion, and
has a great impact on the size–mass relation of galaxies (Figure 14).
Alternatively, a fixed density criterion (e.g., $n_{\rm H}^{*}=0.1$ cm-3) can
be adopted (Katz 1992) although we found that it is better not to include it.

(13)

More sophisticated analytic formulae are also proposed based on smaller-scale
simulations, including the effects of turbulence and magnetic fields
(Federrath et al. 2011).

A fractional part of the mass of the gas particle turns into a new star
particle (see Kobayashi et al. 2007, for the details). Note that an individual
star particle has a typical mass of $\sim 10^{5-7}M_{\odot}$, i.e. it does not
represent a single star but an association of many stars. The masses of the
stars associated with each star particle are distributed according to an IMF.
We adopt Kroupa IMF (with $x=1.3$), which is assumed to be independent of time
and metallicity, and limit the IMF to the mass range $0.01M_{\odot}\leq m\leq
120M_{\odot}$. This assumption is probably valid for the resolution down to
$\sim 10^{4}M_{\odot}$.

Figure 14:  Dependence of the star formation timescale parameter $c_{*}$ and
the IMF slope $x$ on the galaxy scaling relations: size–mass (upper panels)
and mass–metallicity (lower panels) relations. A single slope like Salpeter
IMF assumed. The points show a set of chemodynamical simulations of early-type
galaxies that undergo various merging histories. The solid and dashed lines
indicate the observed relations. Figure is taken from Kobayashi (2005) with
modification.

Kobayashi (2004) introduced a scheme to follow the evolution of the star
particle, as a simple stellar population (SSP), which is defined as a single
generation of coeval and chemically homogeneous stars of various masses, i.e.
it consists of a number of stars with various masses but the same age and
chemical composition. The production of each element $i$ from the star
particle (with an initial mass of $m_{*}^{0}$) is calculated using a very
similar equation as Eqs.(1):

$E_{Z}(t)=m_{*}^{0}\left[E_{\rm SW}+E_{\rm SNcc}+E_{\rm SNIa}+E_{\rm
NSM}\right].$ (14)

With this SSP method, the identical equations for GCE, Eqs. 4-7 can used for
$E_{\rm SW}$, $E_{\rm SNcc}$, $E_{\rm Ia}$, and $E_{\rm NSM}$.

Similarly, the energy production from the star particle is:

$E_{e}(t)=m_{*}^{0}\left[e_{e,{\rm SW}}{\cal{R}}_{\rm SW}(t)+e_{e,{\rm
SNcc}}{\cal{R}}_{\rm SNcc}(t)+e_{e,{\rm SNIa}}{\cal{R}}_{\rm SNIa}(t)\right],$
(15)

where the event rates of SWs (${\cal{R}}_{\rm SW}$) and core-collapse
supernovae (${\cal{R}}_{\rm SNcc}$) can be calculated as

${\cal R}_{\rm
SW}=\int_{m_{t}}^{m_{u}}\,\frac{1}{m}\,\psi(t-\tau_{m})\,\phi(m)~{}dm,$ (16)
${\cal R}_{\rm SNcc}=\int_{\max[m_{{\rm SNcc},\ell},\,m_{t}]}^{\min[m_{{\rm
SNcc},u},m_{u}]}\,\frac{1}{m}\,\psi(t-\tau_{m})\,\phi(m)~{}dm.$ (17)

The lower and upper limits are the same as IMF upper limit for stellar winds
($0.01$ and $120M_{\odot}$), while it is set to be $\sim 8M_{\odot}$
(depending on $Z$) and $50M_{\odot}$ for core-collapse supernovae. The energy
production from SNe Ia is significant, and the rate is given by Eq. 8, while
the energy production from NSMs are negligible. The energy per event (in erg)
for SWs: $e_{e,{\rm SW}}=0.2\times 10^{51}({Z}/{Z_{\odot}})^{0.8}$ for
$>8M_{\odot}$ stars, Type II supernovae: $e_{e,{\rm SNII}}=1\times 10^{51}$,
hypernovae: $e_{e,{\rm HN}}=10,10,20,30\times 10^{51}$ for
$20,25,30,40M_{\odot}$ stars, respectively, and for SNe Ia: $e_{e,{\rm
SNIa}}=1.3\times 10^{51}$ erg.

Figure 15:  Dependence of the star formation timescale parameter $c_{*}$ and
the kinetic feedback fraction $f_{\rm kin}$ on the internal structure of
galaxies: surface brightness profile (left panel) and the stellar metallicity
radial gradient (right panel) for a chemodynamical simulation of an early-type
galaxy with the initial perturbation. A single slope with $x=1.1$ adopted.
Figure is taken from Kobayashi (2004) with modification.

We distribute this feedback metal mass and energy from a star particle $j$
isotropically to a fixed number, $N_{\rm FB}$, of neighbour gas particles
$\ell$, weighted by the smoothing kernel as:

$\frac{d(Z_{i,k}\,m_{{\rm
g},k})}{dt}=\sum_{j}\left[W_{kj}E_{Z,j}(t)/\sum_{\ell}^{N_{\rm FB}}W_{\ell
j}\right]$ (18)

and

$\frac{du_{k}}{dt}=(1-f_{\rm
kin})\sum_{j}\left[W_{kj}E_{e,j}(t)/\sum_{\ell}^{N_{\rm FB}}W_{\ell j}\right]$
(19)

where $u_{k}$ denotes the internal energy of a gas particle $k$, and $f_{\rm
kin}$ denotes the kinetic energy fraction (Navarro & White 1993), i.e., the
fraction of energy that is distributed as a random velocity kick to the gas
particle. Note that in the first sum the number of neighbour star particles is
not fixed to $N_{\rm FB}$. To calculate these equations, the feedback
neighbour search needs to be done twice in order to ensure proper mass and
energy conservation; first to compute the sum of weights for the
normalization, and a second time for the actual distribution.

Figure 16: Time evolution of the formation of an isolated disk galaxy in a
halo of mass of $10^{10}h^{-1}M_{\odot}$. The black points show star
particles, while the gas particles are colour-coded according to their
temperature. Each panel is $20$ kpc on a side. The upper row shows face-on
projections, the lower row gives edge-on views. Figure is taken from Kobayashi
et al. (2007).

The parameter $N_{\rm FB}$ determines the average energy distributed to gas
particles; a large value of $N_{\rm FB}$ leads to inefficient feedback as the
ejected energy radiatively cools away shortly, while a small value of $N_{\rm
FB}$ results in a only small fraction of matter affected by feedback.
Alternatively, with good resolution in zoom-in simulations, feedback neighbors
could be chosen within a fixed radius ($r_{\rm FB}$), which affects a number
of observations including the MDF (Fig. 12 of Kobayashi & Nakasato 2011). The
parameter $f_{\rm kin}$ has a great impact on surface brightness profiles and
radial metallicity gradients in galaxies, and $f_{\rm kin}=0$, i.e., purely
thermal feedback, gives the best match to the observations (Figure 15). Note
that various feedback methods are proposed such as the stochastic feedback
(Dalla Vecchia & Schaye 2008, used for EAGLE) and the mechanical feedback
(Hopkins et al. 2018a, used for FIRE). Feedback modelling can have an impact
on chemical evolution of galaxies (D. Ibrahim et al., in prep.).

Figure 16 shows a test simulation of an isolated disk galaxy with supernova
feedback. The initial condition is a rotating gas cloud with total mass of
$10^{10}h^{-1}M_{\odot}$, and non-isotropic infall of gas form a dense disc,
where stars form. As a result of the energy input, the low-density hot-gas
region expands in a bipolar flow (left panels) ejecting metals outside the
galaxy. In the disk plane, the hot gas region expands and forms a dense shell
where stars keep on forming. The energy from these stars can quickly propagate
to the surrounding low-density region. After the galactic wind forms (middle
panel), the gas density becomes so low at the centre that star formation is
terminated for a while. Because of radiative cooling, a part of the ejected,
metal-enhanced gas, however, returns, again settling in the disk where it
fuels new star formation, but this secondary star formation is not as strong
as the initial starburst. Although some small bubbles are forming in the
galaxy in this stage, not much gas and metals are ejected by them from the
disk (right panels).

Figure 17:  (Left) Spatial distribution of ISM (blue) and wind (red) particles
at $z=0$. (Right) Ejected metal fraction $f_{Z,{\rm w}}=M_{Z,{\rm
w}}/(M_{Z,{\rm g}}+M_{Z,*}+M_{Z,{\rm w}})$ against total mass. Figures are
taken from Kobayashi et al. (2007) with modification.

However, galactic winds become much weaker in more massive galaxies. Figure 17
shows the fraction of metal loss as a function of galaxy mass in a
cosmological simulation including hypernova feedback. Most of metals are
inside the galaxy at $M_{\rm tot}>\kern-11.99998pt\lower
4.73611pt\hbox{$\sim$}~{}10^{12}h^{-1}M_{\odot}$. Eq.(18) gives positive
feedback that enhances radiative cooling, while Eq.(19) gives negative
feedback that suppresses star formation. The mass ejection in Eq.(18) never
becomes zero, while the energy production in Eq.(19) becomes small after $\sim
35$ Myrs. Therefore, it is not easy to control star formation histories
depending on the mass/size of galaxies within this framework. In particular,
it is not possible to quench star formation in massive galaxies, since low-
mass stars keep returning their envelope mass into the ISM for a long
timescale, which will cool and keep forming stars.

Figure 18:  Cosmic BH accretion rate density (solid line) and BH growth rate
density (including BH formation, dashed line) with redshift. The difference of
these rates shows a BH merger rate. Figure is taken from Taylor & Kobayashi
(2014).

Therefore, in order to reproduce observed properties of massive galaxies,
additional feedback was required, and the discovery of the co-evolution of
super-massive black holes and host galaxies (Magorrian et al. 1998) provided a
solution – AGN feedback. Modelling of AGN feedback consists of (1) seed
formation, (2) growth by mergers and gas accretion, and (3) thermal and/or
kinetic feedback (see Taylor & Kobayashi 2014 for the details). In summary, we
introduced a seeding model where seed black holes originate from the formation
of the first stars, which is different from the ‘standard’ model by Springel
et al. (2005) and from most large-scale hydrodynamical simulations. In our AGN
model, seed black holes are generated if the metallicity of the gas cloud is
zero ($Z=0$) and the density is higher than a threshold density ($\rho_{\rm
g}>\rho_{\rm c}$). The growth of black holes is roughly the same as in other
cosmological simulations, and is calculated with Bondi-Hoyle accretion,
swallowing of ambient gas particles, and merging with other black holes. Since
we start from relatively small seeds, the black hole growth is driven by
mergers at $z>\kern-11.99998pt\lower 4.73611pt\hbox{$\sim$}~{}3$ (Figure 18).
On a very small scale, it is not easy to merge two black holes, and an
additional time-delay is applied in more recent simulations (P. Taylor et al.
in prep.).

Proportional to the accretion rate, thermal energy is distributed to the
surrounding gas, which is also the same as in many other simulations. In more
recent simulations, non-isotropic distribution of feedback area is used to
mimic the small-scale jet (e.g., Mukherjee et al. 2018). There are a few
parameters in our chemodynamical simulation code, but Taylor & Kobayashi
(2014) constrained the model parameters from observations, and determined the
best parameter set: $\alpha=1$, $\epsilon_{\,\rm f}=0.25$, $M_{\rm
seed}=10^{3}h^{-1}M_{\odot}$, and $\rho_{\rm c}=0.1h^{2}m_{\rm
H}\,\textrm{cm}^{-3}$. Our black holes seeds are indeed the debris of the
first stars although we explored a parameter space of
$10^{1-5}h^{-1}M_{\odot}$. This is not the only one channel for seeding, and
the direct collapse of primordial gas and/or the collapse of dense stellar
clusters ($\sim 10^{5}M_{\odot}$, e.g., Madau & Rees 2001; Woods et al. 2019)
are rarer but should be included in larger volume simulations.

Nonetheless, our model can successfully drive large-scale galactic winds from
massive galaxies (Taylor & Kobayashi 2015a) and can reproduce many
observations of galaxies with stellar masses of $\sim 10^{9-12}M_{\odot}$
(Taylor & Kobayashi 2015b, 2016, 2017). Metals are ejected from galaxies to
circum-galactic and inter-galactic medium mainly by supernova-driven winds in
low-mass galaxies (as in Fig. 17), while by AGN-driven winds by massive
galaxies (Taylor et al. 2020). The fraction of metals ejected by AGN is only a
few percent of the total produced metals, and thus it does not affect mass-
metallicity relations (Taylor & Kobayashi 2015a). The movie of our fiducial
run is available at https://www.youtube.com/watch?v=jk5bLrVI8Tw

###  3.2 Various big simulations

Figure 19: Cosmic star formation rate history of our $100h^{-1}$ Mpc
simulation, comparing to other cosmological simulations. See text for the
details. The observational data (points) are taken from Madau & Dickinson
(2014) and Driver et al. (2018). Note that other observations (e.g., Chary et
al. 2007; Rowan-Robinson et al. 2016) show higher rates at high redshifts.

There are various large volume simulations available but the input physics is
very different, which can summarised as follows. The differences in modelling
of sub-galactic physics largely affects the predicted galaxy properties, and
the parameters must be calibrated using galaxy scaling relations (e.g., Figure
14) and internal structures (e.g., Figure 15) at $z=0$, before predicting
observables at higher-redshifts using the simulation outputs. Figure 19 shows
a comparison of some of these simulation to our simulation (Kobayashi &
Taylor, in prep.). See also Wang et al. (2019) for the discussion on the black
hole seed mass.

EAGLE (Schaye et al. 2015) – Star formation is calculated with a metallicity
and temperature dependent density threshold ($n_{\rm
H}^{*}(Z)=0.1\left(\frac{Z}{0.002}\right)^{-0.64}$ cm-3 and $T<10^{4.4}$K).
Star particle mass is determined from a pressure law, obtained from the the
observed Kennicutt-Schmidt, with two free parameters ($A$ and $n$). A Chabrier
IMF is adopted. Stochastic feedback from Dalla Vecchia & Schaye (2012)
assumed. SMBHs are grown from the seeds of $10^{5}M_{\odot}$, and the AGN
feedback is thermal, independent of metallicity. The SN Ia rate decays
exponential with a timescale of $\tau=2$ Gyr, with a fixed normalization
$\nu$. The simulation box size is up to 100 Mpc3, which has $1.81\times
10^{6}M_{\odot}$ and 0.7 kpc resolution for gas.

IllustrisTNG (Pillepich et al. 2018) – It is run with a magneto-hydrodynamic
code (Pakmor et al. 2011). Star formation and pressurization of the multi-
phase are treated following Springel & Hernquist (2003), which includes a
density threshold $n_{\rm H}\simeq 0.1$ cm-3. A Chabrier IMF is adopted.
Isotropic, kinetic (wind) feedback from Springel & Hernquist (2003) is assumed
(instead of bipolar winds in Illustris). SMBHs are grown from the seeds of
$8\times 10^{5}M_{\odot}$, and kinetic wind is assumed for AGN feedback in the
case of low accretion rates as in Weinberger et al. (2017). SN Ia rate is
calculated with a “simple DTD” ($-1$ slope?) with a constant normalization
$N_{0}$. The simulation box size is up to 302.6 Mpc3, but a 110.7 Mpc3 run has
$1.4\times 10^{6}M_{\odot}$ and 0.19 kpc resolution for gas.

SIMBA (Davé et al. 2019) – This uses H2-based star formation rate from
Krumholz & Gnedin (2011), which is based on the metallicity and local column
density, assuming the efficiency of $\epsilon_{*}=0.02$ for the H2 density.
Stellar feedback is assumed to be kinetic, with decoupled wind particles,
galaxy mass dependent ($\eta(M_{*})$), and metal-loaded ($dZ=f_{\rm
SNII}\,y_{\rm SNII}/\max[\eta,1]$). SMBHs are grown from the seeds of $1\times
10^{4}h^{-1}M_{\odot}$. AGN feedback is also assumed to be kinetic for both
fast and slow accretion, and purely bipolar. In addition, on-the-fly dust
production and destruction is calculated. SNe Ia have two distinct delay-
times, 0 and 0.7 Gyr, with two constant rates. The simulation box size is up
to $100h^{-1}$ Mpc3, which has $1.82\times 10^{7}M_{\odot}$ and $0.5h^{-1}$
kpc resolution for gas.

HORIZON-AGN (Dubois et al. 2016) – This uses a density threshold of $n_{\rm
H}\simeq 0.1$ cm-3, thermal energy injection to model stellar feedback, and
two modes of AGN feedback with thermal or bipolar outflow depending on the
Eddington ratio. The simulation box size is $100h^{-1}$ Mpc3, which has
$1\times 10^{7}M_{\odot}$, and the grid is adaptively refined down to 1 proper
kpc.

Magneticum (Dolag et al. 2017) – Multiphase model for star formation from
Springel & Hernquist (2003), and isotropic winds with 350 km s-1 are adopted.
Two modes of AGN feedback with different efficiencies of energy injection
depending on the Eddington ratio are assumed. The simulation box size is
$48h^{-1}$ Mpc3, which has $3.7\times 10^{7}h^{-1}M_{\odot}$ and $0.7h^{-1}$
kpc resolution for gas.

Similarly, there are a series of cosmological zoom-in simulations of Milky
Way-type galaxies. Auriga (Grand et al. 2017) uses similar code as Illustris
but with $n=0.13$ cm-3 and $t_{\rm sf}=2.2$ Gyr, with resolution from $6\times
10^{3}$ to $4\times 10^{5}$ for gas. FIRE-2 (Hopkins et al. 2018b) is run with
the GIZMO code, which is also used for SIMBA, but with improved modelling of
small-scale physics. AGN feedback is not included. ARTEMIS (Font et al. 2020)
uses the same code as EAGLE with resolution of $2\times 10^{4}h^{-1}$ for gas.

###  3.3 Galactic archaeology

Figure 20: V-band luminosity map for the edge-on (left panels) and face-on
(right panels) views at $z=0$ in 30 kpc on a side, of our chemodynamical zoom-
in simulation of a Milky Way-type galaxy.

In the Local Group, i.e., in the Milky Way and dwarf spheroidal (dSph)
galaxies, the spatial distribution of detailed elemental abundances can be
obtained from high-resolution spectra of individual low-mass stars777Ages of
stars can also be well estimated with astroseismology recently.. This has been
done for more than half a century, and the average trends of the observations
for many elements have been used to constrain the stellar physics (Timmes et
al. 1995; Kobayashi et al. 2006; Romano et al. 2010). The recent Galactic
archaeology surveys with medium resolution MOS dramatically increased the
sample and made possible to discuss the distributions of stars along the
average trends, which requires more realistic, chemodynamical simulations of
galaxies. Meantime, observations with different lines for each element (e.g.,
Israelian et al. 1998; Sneden et al. 2016) revealed NLTE and 3D effects in the
stellar atmosphere. Detailed star-by-star analysis of wide range of high-
resolution spectra (including UV) is still required to obtain absolute values
of elemental abundances (Kobayashi 2022a, for the need of UV spectra).

Figure 20 shows the edge-on and face-on maps of our simulated Milky Way-type
galaxy, where V-band luminosities of star particles are calculated using
stellar population synthesis models from Kodama & Arimoto (1997) and all
detailed chemical enrichment is included. Although there were simulations of
isolated galaxies, Kobayashi & Nakasato (2011) was the first chemodynamical
simulations that showed the evolution elements from O to Zn in a Milky Way-
type galaxy from cosmological initial condition. Our new simulation includes
all stable elements, and a fully cosmological initial condition is applied
(Scannapieco et al. 2012). The spacial and mass resolutions are 0.5 kpc and
$3\times 10^{5}M_{\odot}$, respectively. The movie of our fiducial run is
available at https://star.herts.ac.uk/~chiaki/works/Aq-C-5-kro2.mpg

Our simulated galaxy shows very similar maps of metallicity and [$\alpha$/Fe]
ratios as observed. Figure 21 shows the observations in Gaia DR3, which
contains a half billion of stars. The data show too blue outer bulge in the
metallicity map, and in the [$\alpha$/Fe] map the patterns close to the
Ecliptic Poles are artefacts. There are both vertical and radial metallicity
gradients; metallicity is high on the plane, and becomes higher toward
galactic centre. On the other hand the [$\alpha$/Fe] ratios show a strong
vertical gradient but no radial gradient except for the galactic bulge; the
majority of the bulge stars have high [$\alpha$/Fe] ratios (Kobayashi &
Nakasato 2011), which may depend on the feedback from galactic centre. There
is no chemodynamical model that includes feedback from the central SMBH in the
Milky Way.

Figure 21: (Left panel) Metallicity map of Gaia DR3 released on 13th June
2022, in the range of [M/H]$=-0.3$ (blue) to 0.0 (red). (Right panel) The same
as the left panel but for [$\alpha$/Fe]$=-0.02$ (blue) to 0.12 (red). Figures
are taken from Gaia Collaboration et al. (2022).

The formation of the basic structure is the same as described in Kobayashi &
Nakasato (2011). The galactic bulge formed by an initial star burst, by the
assembly of gas-rich sub-galaxies beyond $z=2$. The disk formed with a longer
timescale and has grown inside-out; the disk was small in the past and becomes
larger at later times. In the $\Lambda$CDM cosmology, satellite galaxies keep
accreting, which interact and merge with the main body, but there is no major
merger after $z=2$. Otherwise, it is not possible to keep the disk structure
as in the Milky Way. These galaxy mergers also make the disk thicker.
Approximately one third of thick disk stars already formed in merging
galaxies, which are disrupted by tidal force and accreted onto the disk plane.
In the thin disk, star formation is self-regulated, and chemical enrichment
takes place following cosmological gas accretion, radial flows, and stellar
migrations; these physical processes are analysed in detail in Vincenzo &
Kobayashi (2020).

Inhomogeneous enrichment in chemodynamical simulations leads to a paradigm
shift on the chemical evolution of galaxies. As in a real galaxy, i) the star
formation history is not a simple function of radius, ii) the ISM is not
homogeneous at any time, and iii) stars migrate, which are fundamentally
different from one-zone or multi-zones GCE models. As a consequence, (1) there
is no tight age–metallicity relation, namely for stars formed in merging
galaxies. It is possible to form extremely metal-poor stars at a later time,
from accretion of nearly primordial gas, or in isolated chemically-primitive
regions. (2) Enrichment sources with long time-delays such as AGB stars, SNe
Ia, and NSMs can appear at low metallicities. This effect can naturally
explain the observed N/O–O/H relation (Vincenzo & Kobayashi 2018a), but is not
sufficient to explain observed [Eu/(O,Fe)] ratios only with NSMs (Haynes &
Kobayashi 2019). (3) There is a significant scatter in elemental abundance
ratios at a given time/metallicity, as shown in Figs. 16-18 of Kobayashi &
Nakasato (2011).

Figure 22: The [O/Fe]–[Fe/H] relation in the solar neighborhood of our
simulated Milky Way-type galaxy at $z=0$. Figure is the same as Fig. 10 of
Kobayashi & Nakasato (2011) but with a newer simulation. The colour indicates
the number of stars in logarithmic scale. The observational data sources are:
Bensby et al. (open and filled triangles, 2004) for thick and thin disk stars,
respectively, Spite et al. (open circles, 2005), Zhao et al. (filled circles,
2016), Amarsi et al. (crosses, 2019).

Figure 22 shows the frequency distribution of [O/Fe] ratios in the solar
neighborhood of our simulated galaxy. The same [O/Fe]–[Fe/H] relation as in
Fig. 5 is seen: the plateau at [O/Fe] $\sim 0.5$ caused by core-collapse
supernovae, and the decrease of [O/Fe] from [Fe/H] $\sim-1$ to $\sim 0$ due to
the delayed enrichment from SNe Ia. The difference is that this chemodynamical
model predicts not a line but a distribution (contours). A similar figure was
also shown in Kobayashi & Nakasato (2011), which predicted a bimodal
distribution of [O/Fe] ratio at a given [Fe/H], concluding that “this may be
because the mixing of heavy elements among gas particles is not included in
our model”. However, with the APOGEE survey (Hayden et al. 2015) clearly
showed the bimodality of [$\alpha$/Fe], which was in fact already seen in
earlier works with a much smaller sample but careful analysis (Fuhrmann 1998;
Bensby et al. 2003). The advantage of the APOGEE survey was its large dynamic
range and the change of the bimodality depending on the location within the
Galaxy is also clearly shown. Kobayashi (2016) showed a similar change in the
simulated Galaxy, and Vincenzo & Kobayashi (2020) showed a comparison to
APOGEE DR16. Our simulated galaxy predicts this bimodality should exist not
only for [$\alpha$/Fe] ratios but also for most of elements from He to U; some
elements even show trimodality (Kobayashi 2022a; Kobayashi in prep.).

At [Fe/H] $<\kern-11.99998pt\lower 4.73611pt\hbox{$\sim$}~{}-1$ Figure 22 also
shows a significant number of stars with low [$\alpha$/Fe] ratios, which is
caused by local enrichment from SNe Ia. Note that the SN Ia rate becomes
almost zero from the stars below [Fe/H] $=-1.1$ in the adopted SN Ia model.
Nonetheless, stars formed in less-dense regions can have low [$\alpha$/Fe]
caused by SNe Ia, or by low-mass Type II supernovae ($13-15M_{\odot}$); this
effect is important also for dSph galaxies. The number of such low-$\alpha$
stars increases with sub-Ch mass SNe Ia (Kobayashi, in prep.).

The origin of the scatter/bimodality is discussed in Kobayashi (2014). In
chemodynamical simulations, it is possible to trace back the formation place
of star particles. As a result, stars formed in merging galaxies found to have
old age, low metallicity, high [$\alpha$/Fe] ratios, and relatively low
[(Na,Al,Cu)/Fe] ratios, while stars formed in-situ show a tighter age-
metallicity relation and low [$\alpha$/Fe] ratios. There is no age-metallicity
relation for the stars that migrated, as expected. In cosmological zoom-in
simulations, stars tend to migrate outwards by a few kpc, which flatten the
metallicity gradient by $\sim 0.05$ dex/kpc. Radial flows steepen the gradient
as chemically enriched gas moves inwards, but the average velocity is found to
be only $\sim 0.7$ km/s (Vincenzo & Kobayashi 2020).

Figure 23: Distribution of [Eu/Fe] against [Fe/H] for the star particles in
the solar neighbourhood of our simulated Milky Way-type galaxy at $z=0$. The
panels in order show: control, ECSNe+$\nu$-driven winds, NSMs, and MRSNe. The
observational data sources are: Hansen et al. (2016, red squares), Roederer et
al. (2014, orange circles), Zhao et al. (2016, green triangles), and Buder et
al. (2018, GALAH DR2, cyan contours). The contours show 10, 50, and 100 stars
per bin. The colour bar shows the linear number per bin of simulation data.
Figure is taken from Haynes & Kobayashi (2019).

Figure 23 shows the [Eu/Fe] ratios in the solar neighborhood of our simulated
galaxy but with switching the r-process sites. With only AGB stars (panel a),
Eu is not sufficiently produced. Different from one-zone GCE models, AGB
contribution can be seen with a large scatter at low metallicity. With ECSNe
or $\nu$-driven winds (panel b), Eu production is not increased enough. With
NSMs (panel c), it is possible to reproduce the solar Eu/Fe ratios at the
solar metallicity, but the scatter is too large at [Fe/H]
$<\kern-11.99998pt\lower 4.73611pt\hbox{$\sim$}~{}0$. Unlike one-zone GCE
models, there are a small number of stars that have high [Eu/Fe] ratios at low
metallicities due to the inhomogeneous enrichment. The scatter can be much
reduced with MRSNe (panel d) as they occur at a very short timescale. Since
only a small fraction of core-collapse supernovae produce both Eu and Fe, the
small scatter still remains, consistent with observations.

Very high [Eu/Fe] ratios were reported for ultrafaint dSphs (Ji et al. 2016),
which might be caused by local enrichment with a NSM with no Fe production
(and no Zn enhancement unlike Yong et al. (2021)’s star). The effect of
supernova kick in binary neutron star systems may help (van de Voort et al.
2022). Although supernova kick is included in the delay-time distributions
from BPSs, it is difficult to include it in the sub-galactic scale. In dSph
galaxies, star formation takes place slowly, chemical enrichment proceeds
inefficiently, and thus the inhomogeneous enrichment effect becomes even more
important. Higher resolution of simulations with a new formalism/code will be
required.

###  3.4 Extra-galactic archaeology

Figure 24: (Left panel) The first NIRCam/JWST image of a cluster of galaxy
SMAC0723 at $z=0.39$ but showing much higher redshift galaxies too, released
on 11th July 2022. Figure is taken from https://webbtelescope.org. (Right
panel) Emission lines that can be obtained with NIRSpec/JWST as a function of
wavelength and redshift.

Elemental abundances and isotopic ratios provide additional constrains on the
timescales of formation and evolution of galaxies. In galaxies beyond the
Local Group, it is not usually possible to resolve stars or HII regions, and
metallicites of stellar populations888Multiple lines should be used to break
the age-metallicity degeneracy. or ISM are estimated from absorption or
emission lines in integrated spectra. Thanks to MOS surveys since SDSS the
sample is greatly increased, and thanks to IFU surveys since SAURON the
spacial (projected) distributions of metallicities are also obtained. Note
that the methods for obtaining the absolute values of physical quantities are
still debated (Worthey et al. 1992; Conroy 2013; Maiolino & Mannucci 2019;
Kewley et al. 2019). For stellar populations, $\alpha$/Fe ratios have been
used to estimate the formation timescale of early-type galaxies (Thomas et al.
2005; Kriek et al. 2016), while Fe is not accessible in star-forming galaxies,
and instead CNO abundances can be used (Vincenzo & Kobayashi 2018b). The James
Webb Space Telescope (JWST) will push these observations toward higher
redshifts, and the wavelength coverage of which will allow us to measure CNO
abundance simultaneously (Figure 24). More elements are available for X-ray
hot gas or quasar absorption line systems, although neutron capture elements
are out of reach. Isotopic ratios of light elements and some light element
abundances (e.g., F) are also estimated with Atacama Large Millimeter Array
(ALMA).

Figure 25: The time evolution of our cosmological simulation in a periodic box
$50h^{-1}$ Mpc on a side. We show the projected stellar $V$-luminosity (upper
panels) and gas metallicity $\log Z_{\rm g}/Z_{\odot}$ (lower panels).

Figure 25 shows the evolution of stellar luminosities and metallicities in a
cosmological simulation. More massive galaxies tend to have higher
metallicities, which is called mass–metallicity relations (MZR). As already
discussed, the origin of these relations are mass-dependent galactic winds by
supernova feedback (Kobayashi et al. 2007; Taylor et al. 2020). More intense
star formation also happens in the simulated massive galaxies. An IMF
variation is also possible (Kobayashi 2010). Metal-loss by AGN feedback is not
significant (Taylor & Kobayashi 2015a).

Various ways are possible to calculate the abundance ratio of the stellar
population (or ISM) of a galaxy:

$\displaystyle\log
X/Y=\log\left(\sum_{N}M_{X}/\sum_{N}M_{Y}\right),\log\left(\sum_{N}M\frac{M_{X}}{M_{Y}}/\sum_{N}M\right),\left(\sum_{N}M\,\log\frac{M_{X}}{M_{Y}}\right)/\sum_{N}M,$
$\displaystyle\log\left(\sum_{N}\frac{LM_{X}}{M}/\sum_{N}\frac{LM_{Y}}{M}\right),\log\left(\sum_{N}L\frac{M_{X}}{M_{Y}}/\sum_{N}L\right),\mathrm{or}\left(\sum_{N}L\,\log\frac{M_{X}}{M_{Y}}\right)/\sum_{N}L,$
(20)

where $N$ is the number of stars/star particles (or HII regions/gas
particles). The first three are mass-weighted, but in order to compare with
observational data, it should be weighted by luminosity $L$ (or star formation
rate). For stellar populations, we use V-band luminosity where Mgb line
(5143-5206 Å) is located, and these definitions give slightly different
results. For the ISM, the star formation rate of each gas particle can be used
for weighting since it correlates with H$\alpha$ emissions. This weighting
gives very different results as there are many non-star forming gas particles
in galaxies (see Fig. 19 of Kobayashi et al. 2007). It is not obvious which is
the best definition.

Within galaxies, central parts are more metal-rich than outskirts, which is
called metallicity radial gradients. The origin of these gradients are the
inside-out formation, where star formation starts from the centre with a
higher efficiency and more chemical enrichment cycles. The star formation is
more intense at the centre in the simulated galaxies. The star formation
duration is not necessarily longer in the centre if quenching happens also
inside-out by AGN feedback. Radial flows steepen the gradients, while galaxy
mergers flatten the gradients (Kobayashi 2004). Re-distribution of metals by
AGN feedback may be important for gas-phase metallicity gradients (Taylor &
Kobayashi 2017). Simulations predict a significant scatter of gradients at a
given mass (Kobayashi 2004; Taylor & Kobayashi 2017), while the average of the
gradients become steeper at higher redshifts (Pilkington et al. 2012;
Kobayashi & Taylor, in prep.).

The metallicity gradient is defined in a log-scale for stellar populations but
in a linear scale for gas-phase, and their example values are:

$\Delta\log Z_{*}/\Delta\log r\sim 0.3$ (21)

for nearby early-type galaxies and

$\Delta Z_{\rm g}/\Delta r\sim 0.05\mathrm{~{}[dex/kpc]}$ (22)

for nearby disks, where $r$ is the projected, galactocentric radius. These
definitions give a good fit to slit observations of nearby galaxies (e.g.,
Zaritsky et al. 1994; Kobayashi & Arimoto 1999) that cover a large radius
($r<\kern-11.99998pt\lower 4.73611pt\hbox{$\sim$}~{}2r_{\rm e}$), where
$r_{\rm e}$ is the effective radius at which a half of stellar luminosity is
enclosed. Often, the central part ($r<\kern-11.99998pt\lower
4.73611pt\hbox{$\sim$}~{}1$ kpc) is excluded from the fitting because it shows
a flattening due to a limited seeing in observations, or due to a limited
spatial resolution in simulations. The outer part also show a sudden decrease
in metallicity, or an increase due to a satellite galaxies, which is also
excluded from the fitting. However, recent analysis of IFU surveys stellar
metallicities are also measured against linear $r$, and it is important to
clarify the definition.

Figure 26: (Left) Evolution of gas-phase CNO abundance ratios in our simulated
disk galaxies with different star formation timescales, taken from a
cosmological simulation. Figure is taken from Vincenzo & Kobayashi (2018b).
(Right) Evolution of the gas-phase N/O–O/H relation of our simulated galaxies
in a cosmological simulation. The gray bar indicates the compilation of
observational data. Figure is taken from Vincenzo et al. (2019).

Our cosmological simulations also predict elemental abundances and isotopic
ratios. The left panel of Figure 26 shows theoretical predictions of CNO
abundance ratios for simulated disk galaxies that have different star
formation timescales. These galaxies are chosen from a cosmological
simulation, which is run by a Gadget-3 based code that includes detailed
chemical enrichment (Kobayashi et al. 2007). In the simulation, C is mainly
produced from low-mass stars ($<\kern-11.99998pt\lower
4.73611pt\hbox{$\sim$}~{}4M_{\odot}$), N is from intermediate-mass stars
($>\kern-11.99998pt\lower 4.73611pt\hbox{$\sim$}~{}4M_{\odot}$) as a primary
process (Kobayashi et al. 2011a), and O is from massive stars
($>\kern-11.99998pt\lower 4.73611pt\hbox{$\sim$}~{}13M_{\odot}$). C and N are
also produced by massive stars; the N yield depends on the metallicity as a
secondary process, and can be greatly enhanced by stellar rotation (as for F).

In the nearby universe, the N/O–O/H relation is known for stellar and ISM
abundances, which shows a plateau (N/O $\sim-1.6$) at low metallicities and a
rapid increase toward higher metallicities. This was interpreted as the
necessity of rotating massive stars by Chiappini et al. (2006). However, this
should be studied with hydrodynamical simulations including detailed chemical
enrichment, and Kobayashi (2014) first showed the N/O–O/H relation in a
chemodynamical simulation. Vincenzo & Kobayashi (2018a) showed that both the
global relation, which is obtained for average abundances of the entire
galaxies, and the local relation, which is obtained for spatially resolved
abundances from IFU data, can be reproduced by the inhomogeneous enrichment
from AGB stars. Since N yield increases at higher metallicities, the global
relation originates from the mass–metallicity relation of galaxies, while the
local relation is caused by radial metallicity gradients within galaxies.
Moreover, the right panel of Figure 26 shows a theoretical prediction on the
time evolution of the N/O–O/H relation, where galaxies evolve along the
relation. Recent observation with KMOS on VLT confirmed a near redshift-
invariant N/O-O/H relation (Hayden-Pawson et al. 2021, KLEVER survey).

Figure 27: Stellar [O/Fe] ratios as a function of central velocity dispersion,
which is a proxy of galaxy mass, for our simulated galaxies with (black
asterisks) and without (red diamonds) AGN feedback at $z=0$. The arrows
indicate the differences in two cosmological simulations from the same initial
conditions. The solid lines indicate the observed relations. Figure is taken
from Taylor & Kobayashi (2015b).

Figure 27 shows [$\alpha$/Fe] ratios of stellar populations as a function of
galaxy mass for our simulated galaxies with and without AGN feedback. Without
AGN feedback (red diamonds), since star formation lasts longer in massive
galaxies with a deeper potential well, [$\alpha$/Fe] ratios become lower in
massive galaxies, which is the opposite compared with the observations (solid
lines). With AGN feedback (black asterisks) star formation can be suppressed
before the SN Ia enrichment becomes dominant in massive galaxies, so that
[$\alpha$/Fe] ratios can stay high. However, the scatter of [$\alpha$/Fe]
ratios at a given galaxy mass is still larger than observed.

[$\alpha$/Fe] ratios become higher at higher redshifts in the simulations, but
not as much as observed (e.g., Kriek et al. 2016). Some cosmological
simulations claimed that they could reproduce the observed [$\alpha$/Fe]–mass
relation at $z=0$, but with an ad-hoc SN Ia model. Some even introduced a
variation in the IMF or in the binary fraction. However, the [$\alpha$/Fe]
problem should be discussed with a chemodynamical model that is based on
nuclear astrophysics and that can reproduce the observed elemental abundances
both in the Milky Way and dSph galaxies. We do not have such a model yet.

##  4 Conclusions and Future Prospects

Thanks to the long-term collaborations between nuclear and astrophysics, we
have good understanding on the origin of elements (except for the elements
around Ti and a few neutron-capture elements such as Au). Inhomogeneous
enrichment is extremely important for interpreting the elemental abundance
trends. It can reproduce the observed N/O–O/H relation only with AGB stars and
supernovae (Fig. 26), but not the observed r-process abundances only with NSMs
(Fig. 23); an r-process associated with core-collapse supernovae such as
magneto-rotational hypernovae is required (Fig. 11), although the explosion
mechanism is unknown. It is necessary to run chemodynamical simulations from
cosmological initial conditions, including detailed chemical enrichment based
on nuclear astrophysics.

The impact of stellar rotation, binaries, and magnetic fields during stellar
evolution, and the multi-dimensional effects of supernova explosions in
nucleosynthesis should be investigated further. If Wolf-Rayet stars are
producing heavy elements on a very short timescale, it might be hard to find
very metal-poor or Population III (and dust-free) galaxies at very high
redshifts even with JWST. In that case the elemental abundances can be quite
different in very high redshift galaxies (e.g., high (C,N)/O; Fig. 9).
Although there is no observational evidence, pair-instability supernovae could
also cause very different abundance pattern (e.g., high (Si,S)/O; Fig. 10).
Finally, some metals are locked in a solid state - it is also important to
calculate element-by-element dust formation, growth, and destruction, as well
as the detailed chemical enrichment.

Galactic archaeology is a powerful approach for reconstructing the formation
history of the Milky Way and its satellite galaxies. APOGEE and HERMES-GALAH
surveys have provided homogeneous datasets of many elemental abundances that
can be statistically compared with chemodynamical simulations. Future surveys
with WEAVE and 4MOST will provide more. Having said that, the number of EMP
stars will not be increased so much in these surveys, and a target survey such
as the SkyMapper EMP survey is also needed, in particular for constraining the
early chemical enrichment from the first stars. It is important to increase
spectral resolution and wavelength coverage (including UV), to obtain more
accurate abundances of more elements (namely neutron-capture elements).

Reflecting the difference in the formation timescale, elemental abundances
depend on the location within galaxies. Although this dependence has been
explored toward the Galactic bulge by APOGEE, the dependence at the outer disk
is still unknown, which requires 8m-class telescopes such as the PFS on Subaru
telescope. Despite the limited spectral resolution of the PFS, $\alpha$/Fe and
a small number of elements will be available. The PFS will also be able to
explore the $\alpha$/Fe bimodality in M31; it is not yet known if M31 has a
similar $\alpha$/Fe dichotomy or not.

The next step will be to apply the Galactic archaeology approach to external
or distant galaxies. Although it became possible to map metallicity, some
elemental abundances, and kinematics within galaxies with IFU, the sample
and/or spacial resolution are still limited even with JWST. Integrated
physical quantities over galaxies, or stacked quantities at a given mass bin,
will also be useful, which can be done with the same MOS developed for
Galactic archaeology (although optimal spectral resolutions and wavelength
coverages are different). Spectroscopic surveys with 8m-class telescopes will
be useful, and it is a matter of urgency to establish the analysis methods to
obtain absolute values of metallicities and elemental abundances from the
observational data. In addition, ALMA has opened a new window for elemental
abundances and isotopic ratios in high-redshift galaxies. Extra-galactic
archaeology will become popular in coming years.

###### Acknowledgements.

We thank D. Yong, F. Vincenzo, A. Karakas, M. Lugaro, N. Tominaga, S.-C.
Leung, M. Ishigaki, K. Nomoto, L. Kewley, R. Maiolino, A. Bunker for fruitful
discussion, and V. Springel for providing Gadget-3. CK acknowledge funding
from the UK Science and Technology Facility Council through grant
ST/M000958/1, ST/R000905/1, ST/V000632/1. The work was also funded by a
Leverhulme Trust Research Project Grant on “Birth of Elements”.

## References

  * Amarsi et al. (2019) Amarsi, A. M., Nissen, P. E., & Skúladóttir, Á. 2019, A&A, 630, A104, doi: 10.1051/0004-6361/201936265
  * Arcones et al. (2007) Arcones, A., Janka, H. T., & Scheck, L. 2007, A&A, 467, 1227, doi: 10.1051/0004-6361:20066983
  * Arnaboldi et al. (2022) Arnaboldi, M., Bhattacharya, S., Gerhard, O., et al. 2022, arXiv e-prints, arXiv:2208.02328.  https://arxiv.org/abs/2208.02328
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARAA, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222
  * Audouze & Silk (1995) Audouze, J., & Silk, J. 1995, ApJ, Letter, 451, L49, doi: 10.1086/309687
  * Barkat et al. (1967) Barkat, Z., Rakavy, G., & Sack, N. 1967, Physical Review Letters, 18, 379, doi: 10.1103/PhysRevLett.18.379
  * Bensby et al. (2003) Bensby, T., Feltzing, S., & Lundström, I. 2003, A&A, 410, 527, doi: 10.1051/0004-6361:20031213
  * Bensby et al. (2004) —. 2004, A&A, 415, 155, doi: 10.1051/0004-6361:20031655
  * Blondin et al. (2022) Blondin, S., Bravo, E., Timmes, F. X., Dessart, L., & Hillier, D. J. 2022, A&A, 660, A96, doi: 10.1051/0004-6361/202142323
  * Brook et al. (2012) Brook, C. B., Stinson, G. S., Gibson, B. K., et al. 2012, MNRAS, 426, 690, doi: 10.1111/j.1365-2966.2012.21738.x
  * Brown et al. (2001) Brown, G. E., Heger, A., Langer, N., et al. 2001, New Astronomy, 6, 457, doi: 10.1016/S1384-1076(01)00077-X
  * Buck et al. (2020) Buck, T., Obreja, A., Macciò, A. V., et al. 2020, MNRAS, 491, 3461, doi: 10.1093/mnras/stz3241
  * Buder et al. (2018) Buder, S., Asplund, M., Duong, L., et al. 2018, MNRAS, 478, 4513, doi: 10.1093/mnras/sty1281
  * Burkert & Hensler (1987) Burkert, A., & Hensler, G. 1987, in Nuclear Astrophysics, ed. W. Hillebrandt, R. Kuhfuss, E. Mueller, & J. W. Truran, Vol. 287, 159, doi: 10.1007/BFb0016576
  * Burrows & Vartanyan (2021) Burrows, A., & Vartanyan, D. 2021, Nature, 589, 29, doi: 10.1038/s41586-020-03059-w
  * Busso et al. (1999) Busso, M., Gallino, R., & Wasserburg, G. J. 1999, ARAA, 37, 239, doi: 10.1146/annurev.astro.37.1.239
  * Cappellari et al. (2013) Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2013, MNRAS, 432, 1862, doi: 10.1093/mnras/stt644
  * Cayrel et al. (2004) Cayrel, R., Depagne, E., Spite, M., et al. 2004, A&A, 416, 1117, doi: 10.1051/0004-6361:20034074
  * Cen & Ostriker (1999) Cen, R., & Ostriker, J. P. 1999, ApJ, Letter, 519, L109, doi: 10.1086/312123
  * Cescutti & Chiappini (2014) Cescutti, G., & Chiappini, C. 2014, A&A, 565, A51, doi: 10.1051/0004-6361/201423432
  * Cescutti & Kobayashi (2017) Cescutti, G., & Kobayashi, C. 2017, A&A, 607, A23, doi: 10.1051/0004-6361/201731398
  * Chabrier (2003) Chabrier, G. 2003, The Publications of the Astronomical Society of the Pacific, 115, 763, doi: 10.1086/376392
  * Chary et al. (2007) Chary, R., Berger, E., & Cowie, L. 2007, ApJ, 671, 272, doi: 10.1086/522692
  * Chiappini et al. (2006) Chiappini, C., Hirschi, R., Meynet, G., et al. 2006, A&A, 449, L27, doi: 10.1051/0004-6361:20064866
  * Chiappini et al. (1997) Chiappini, C., Matteucci, F., & Gratton, R. 1997, ApJ, 477, 765, doi: 10.1086/303726
  * Conroy (2013) Conroy, C. 2013, ARAA, 51, 393, doi: 10.1146/annurev-astro-082812-141017
  * Cristallo et al. (2011) Cristallo, S., Piersanti, L., Straniero, O., et al. 2011, ApJS, 197, 17, doi: 10.1088/0067-0049/197/2/17
  * Croton et al. (2006) Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11, doi: 10.1111/j.1365-2966.2005.09675.x
  * Dalla Vecchia & Schaye (2008) Dalla Vecchia, C., & Schaye, J. 2008, MNRAS, 387, 1431, doi: 10.1111/j.1365-2966.2008.13322.x
  * Dalla Vecchia & Schaye (2012) —. 2012, MNRAS, 426, 140, doi: 10.1111/j.1365-2966.2012.21704.x
  * Davé et al. (2019) Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827, doi: 10.1093/mnras/stz937
  * De Donder & Vanbeveren (2004) De Donder, E., & Vanbeveren, D. 2004, New Astronomy Reviews, 48, 861, doi: 10.1016/j.newar.2004.07.001
  * Doherty et al. (2014) Doherty, C. L., Gil-Pons, P., Lau, H. H. B., Lattanzio, J. C., & Siess, L. 2014, MNRAS, 437, 195, doi: 10.1093/mnras/stt1877
  * Doherty et al. (2015) Doherty, C. L., Gil-Pons, P., Siess, L., Lattanzio, J. C., & Lau, H. H. B. 2015, MNRAS, 446, 2599, doi: 10.1093/mnras/stu2180
  * Dolag et al. (2017) Dolag, K., Mevius, E., & Remus, R.-S. 2017, Galaxies, 5, 35, doi: 10.3390/galaxies5030035
  * Driver et al. (2018) Driver, S. P., Andrews, S. K., da Cunha, E., et al. 2018, MNRAS, 475, 2891, doi: 10.1093/mnras/stx2728
  * Dubois et al. (2012) Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, MNRAS, 420, 2662, doi: 10.1111/j.1365-2966.2011.20236.x
  * Dubois et al. (2016) Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948, doi: 10.1093/mnras/stw2265
  * Federrath et al. (2011) Federrath, C., Sur, S., Schleicher, D. R. G., Banerjee, R., & Klessen, R. S. 2011, ApJ, 731, 62, doi: 10.1088/0004-637X/731/1/62
  * Few et al. (2014) Few, C. G., Courty, S., Gibson, B. K., Michel-Dansac, L., & Calura, F. 2014, MNRAS, 444, 3845, doi: 10.1093/mnras/stu1709
  * Font et al. (2020) Font, A. S., McCarthy, I. G., Poole-Mckenzie, R., et al. 2020, MNRAS, 498, 1765, doi: 10.1093/mnras/staa2463
  * Franco et al. (2021) Franco, M., Coppin, K. E. K., Geach, J. E., et al. 2021, Nature Astronomy, 5, 1240, doi: 10.1038/s41550-021-01515-9
  * Frischknecht et al. (2016) Frischknecht, U., Hirschi, R., Pignatari, M., et al. 2016, MNRAS, 456, 1803, doi: 10.1093/mnras/stv2723
  * Fuhrmann (1998) Fuhrmann, K. 1998, A&A, 338, 161
  * Gaia Collaboration et al. (2022) Gaia Collaboration, Recio-Blanco, A., Kordopatis, G., et al. 2022, arXiv e-prints, arXiv:2206.05534.  https://arxiv.org/abs/2206.05534
  * Grand et al. (2017) Grand, R. J. J., Gómez, F. A., Marinacci, F., et al. 2017, MNRAS, 467, 179, doi: 10.1093/mnras/stx071
  * Grichener et al. (2021) Grichener, A., Kobayashi, C., & Soker, N. 2021, arXiv e-prints, arXiv:2112.08301.  https://arxiv.org/abs/2112.08301
  * Grisoni et al. (2019) Grisoni, V., Matteucci, F., Romano, D., & Fu, X. 2019, MNRAS, 489, 3539, doi: 10.1093/mnras/stz2428
  * Grisoni et al. (2017) Grisoni, V., Spitoni, E., Matteucci, F., et al. 2017, MNRAS, 472, 3637, doi: 10.1093/mnras/stx2201
  * Gunawardhana et al. (2011) Gunawardhana, M. L. P., Hopkins, A. M., Sharp, R. G., et al. 2011, MNRAS, 415, 1647, doi: 10.1111/j.1365-2966.2011.18800.x
  * Haardt & Madau (1996) Haardt, F., & Madau, P. 1996, ApJ, 461, 20, doi: 10.1086/177035
  * Hansen et al. (2016) Hansen, T. T., Andersen, J., Nordström, B., et al. 2016, A&A, 586, A160, doi: 10.1051/0004-6361/201527235
  * Hayden et al. (2015) Hayden, M. R., Bovy, J., Holtzman, J. A., et al. 2015, ApJ, 808, 132, doi: 10.1088/0004-637X/808/2/132
  * Hayden-Pawson et al. (2021) Hayden-Pawson, C., Curti, M., Maiolino, R., et al. 2021, arXiv e-prints, arXiv:2110.00033.  https://arxiv.org/abs/2110.00033
  * Haynes & Kobayashi (2019) Haynes, C. J., & Kobayashi, C. 2019, MNRAS, 483, 5123, doi: 10.1093/mnras/sty3389
  * Herwig (2005) Herwig, F. 2005, ARAA, 43, 435, doi: 10.1146/annurev.astro.43.072103.150600
  * Hirschi (2007) Hirschi, R. 2007, A&A, 461, 571, doi: 10.1051/0004-6361:20065356
  * Hopkins et al. (2006) Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2006, ApJS, 163, 1, doi: 10.1086/499298
  * Hopkins et al. (2018a) Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018a, MNRAS, 477, 1578, doi: 10.1093/mnras/sty674
  * Hopkins et al. (2018b) —. 2018b, MNRAS, 480, 800, doi: 10.1093/mnras/sty1690
  * Ishigaki et al. (2018) Ishigaki, M. N., Tominaga, N., Kobayashi, C., & Nomoto, K. 2018, ApJ, 857, 46, doi: 10.3847/1538-4357/aab3de
  * Israelian et al. (1998) Israelian, G., García López, R. J., & Rebolo, R. 1998, ApJ, 507, 805, doi: 10.1086/306351
  * Janka (2012) Janka, H.-T. 2012, Annual Review of Nuclear and Particle Science, 62, 407, doi: 10.1146/annurev-nucl-102711-094901
  * Ji et al. (2016) Ji, A. P., Frebel, A., Chiti, A., & Simon, J. D. 2016, Nature, 531, 610, doi: 10.1038/nature17425
  * Just et al. (2015) Just, O., Bauswein, A., Ardevol Pulpillo, R., Goriely, S., & Janka, H. T. 2015, MNRAS, 448, 541, doi: 10.1093/mnras/stv009
  * Karakas (2010) Karakas, A. I. 2010, MNRAS, 403, 1413, doi: 10.1111/j.1365-2966.2009.16198.x
  * Karakas & Lattanzio (2014) Karakas, A. I., & Lattanzio, J. C. 2014, PASA, 31, e030, doi: 10.1017/pasa.2014.21
  * Karakas & Lugaro (2016) Karakas, A. I., & Lugaro, M. 2016, ApJ, 825, 26, doi: 10.3847/0004-637X/825/1/26
  * Katz (1992) Katz, N. 1992, ApJ, 391, 502, doi: 10.1086/171366
  * Kawata & Gibson (2003) Kawata, D., & Gibson, B. K. 2003, MNRAS, 340, 908, doi: 10.1046/j.1365-8711.2003.06356.x
  * Kemp et al. (2022) Kemp, A. J., Karakas, A. I., Casey, A. R., Kobayashi, C., & Izzard, R. G. 2022, MNRAS, 509, 1175, doi: 10.1093/mnras/stab3103
  * Kewley et al. (2019) Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARAA, 57, 511, doi: 10.1146/annurev-astro-081817-051832
  * Kobayashi (2002) Kobayashi, C. 2002, PhD Thesis, University of Tokyo
  * Kobayashi (2004) —. 2004, MNRAS, 347, 740, doi: 10.1111/j.1365-2966.2004.07258.x
  * Kobayashi (2005) —. 2005, MNRAS, 361, 1216, doi: 10.1111/j.1365-2966.2005.09248.x
  * Kobayashi (2010) Kobayashi, C. 2010, in American Institute of Physics Conference Series, Vol. 1240, Hunting for the Dark: the Hidden Side of Galaxy Formation, ed. V. P. Debattista & C. C. Popescu, 123–126, doi: 10.1063/1.3458465
  * Kobayashi (2014) Kobayashi, C. 2014, in Setting the scene for Gaia and LAMOST, ed. S. Feltzing, G. Zhao, N. A. Walton, & P. Whitelock, Vol. 298, 167–178, doi: 10.1017/S1743921313006339
  * Kobayashi (2016) Kobayashi, C. 2016, in The General Assembly of Galaxy Halos: Structure, Origin and Evolution, ed. A. Bragaglia, M. Arnaboldi, M. Rejkuba, & D. Romano, Vol. 317, 57–63, doi: 10.1017/S1743921315009783
  * Kobayashi (2022a) —. 2022a, arXiv e-prints, arXiv:2203.01980.  https://arxiv.org/abs/2203.01980
  * Kobayashi (2022b) —. 2022b, Experimental Astronomy, doi: 10.1007/s10686-022-09862-9
  * Kobayashi & Arimoto (1999) Kobayashi, C., & Arimoto, N. 1999, ApJ, 527, 573, doi: 10.1086/308092
  * Kobayashi et al. (2014) Kobayashi, C., Ishigaki, M. N., Tominaga, N., & Nomoto, K. 2014, ApJ, Letter, 785, L5, doi: 10.1088/2041-8205/785/1/L5
  * Kobayashi et al. (2020a) Kobayashi, C., Karakas, A. I., & Lugaro, M. 2020a, ApJ, 900, 179, doi: 10.3847/1538-4357/abae65
  * Kobayashi et al. (2011a) Kobayashi, C., Karakas, A. I., & Umeda, H. 2011a, MNRAS, 414, 3231, doi: 10.1111/j.1365-2966.2011.18621.x
  * Kobayashi et al. (2020b) Kobayashi, C., Leung, S.-C., & Nomoto, K. 2020b, ApJ, 895, 138, doi: 10.3847/1538-4357/ab8e44
  * Kobayashi & Nakasato (2011) Kobayashi, C., & Nakasato, N. 2011, ApJ, 729, 16, doi: 10.1088/0004-637X/729/1/16
  * Kobayashi & Nomoto (2009) Kobayashi, C., & Nomoto, K. 2009, ApJ, 707, 1466, doi: 10.1088/0004-637X/707/2/1466
  * Kobayashi et al. (2015) Kobayashi, C., Nomoto, K., & Hachisu, I. 2015, ApJ, Letter, 804, L24, doi: 10.1088/2041-8205/804/1/L24
  * Kobayashi et al. (2007) Kobayashi, C., Springel, V., & White, S. D. M. 2007, MNRAS, 376, 1465, doi: 10.1111/j.1365-2966.2007.11555.x
  * Kobayashi et al. (2011b) Kobayashi, C., Tominaga, N., & Nomoto, K. 2011b, ApJ, Letter, 730, L14, doi: 10.1088/2041-8205/730/2/L14
  * Kobayashi et al. (2000) Kobayashi, C., Tsujimoto, T., & Nomoto, K. 2000, ApJ, 539, 26, doi: 10.1086/309195
  * Kobayashi et al. (1998) Kobayashi, C., Tsujimoto, T., Nomoto, K., Hachisu, I., & Kato, M. 1998, ApJ, Letter, 503, L155, doi: 10.1086/311556
  * Kobayashi et al. (2006) Kobayashi, C., Umeda, H., Nomoto, K., Tominaga, N., & Ohkubo, T. 2006, ApJ, 653, 1145, doi: 10.1086/508914
  * Kobayashi et al. (2022) Kobayashi, C., Mandel, I., Belczynski, K., et al. 2022, arXiv e-prints, arXiv:2211.04964.  https://arxiv.org/abs/2211.04964
  * Kodama & Arimoto (1997) Kodama, T., & Arimoto, N. 1997, A&A, 320, 41.  https://arxiv.org/abs/astro-ph/9609160
  * Kraft et al. (1997) Kraft, R. P., Sneden, C., Smith, G. H., et al. 1997, AJ, 113, 279, doi: 10.1086/118251
  * Kriek et al. (2016) Kriek, M., Conroy, C., van Dokkum, P. G., et al. 2016, Nature, 540, 248, doi: 10.1038/nature20570
  * 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.  https://arxiv.org/abs/0708.1164
  * Kroupa et al. (1993) Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545, doi: 10.1093/mnras/262.3.545
  * Krumholz & Gnedin (2011) Krumholz, M. R., & Gnedin, N. Y. 2011, ApJ, 729, 36, doi: 10.1088/0004-637X/729/1/36
  * Leung & Nomoto (2018) Leung, S.-C., & Nomoto, K. 2018, ApJ, 861, 143, doi: 10.3847/1538-4357/aac2df
  * Leung & Nomoto (2020) —. 2020, ApJ, 888, 80, doi: 10.3847/1538-4357/ab5c1f
  * Limongi & Chieffi (2018) Limongi, M., & Chieffi, A. 2018, ApJS, 237, 13, doi: 10.3847/1538-4365/aacb24
  * Lind et al. (2022) Lind, K., Nordlander, T., Wehrhahn, A., et al. 2022, A&A, 665, A33, doi: 10.1051/0004-6361/202142195
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARAA, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Madau & Rees (2001) Madau, P., & Rees, M. J. 2001, ApJ, Letter, 551, L27, doi: 10.1086/319848
  * Maeda & Nomoto (2003) Maeda, K., & Nomoto, K. 2003, ApJ, 598, 1163, doi: 10.1086/378948
  * Magorrian et al. (1998) Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285, doi: 10.1086/300353
  * Maiolino & Mannucci (2019) Maiolino, R., & Mannucci, F. 2019, The Astronomy and Astrophysics Review, 27, 3, doi: 10.1007/s00159-018-0112-2
  * Matteucci (2001) Matteucci, F. 2001, The chemical evolution of the Galaxy, Vol. 253, doi: 10.1007/978-94-010-0967-6
  * Matteucci (2021) —. 2021, The Astronomy and Astrophysics Review, 29, 5, doi: 10.1007/s00159-021-00133-8
  * Matteucci & Brocato (1990) Matteucci, F., & Brocato, E. 1990, ApJ, 365, 539, doi: 10.1086/169508
  * Mennekens & Vanbeveren (2014) Mennekens, N., & Vanbeveren, D. 2014, A&A, 564, A134, doi: 10.1051/0004-6361/201322198
  * Mennekens & Vanbeveren (2016) —. 2016, A&A, 589, A64, doi: 10.1051/0004-6361/201628193
  * Meynet & Maeder (2002) Meynet, G., & Maeder, A. 2002, A&A, 390, 561, doi: 10.1051/0004-6361:20020755
  * Mihos & Hernquist (1996) Mihos, J. C., & Hernquist, L. 1996, ApJ, 464, 641, doi: 10.1086/177353
  * Mösta et al. (2018) Mösta, P., Roberts, L. F., Halevi, G., et al. 2018, ApJ, 864, 171, doi: 10.3847/1538-4357/aad6ec
  * Mukherjee et al. (2018) Mukherjee, D., Bicknell, G. V., Wagner, A. Y., Sutherland, R. S., & Silk, J. 2018, MNRAS, 479, 5544, doi: 10.1093/mnras/sty1776
  * Navarro & White (1993) Navarro, J. F., & White, S. D. M. 1993, MNRAS, 265, 271, doi: 10.1093/mnras/265.2.271
  * Navarro & White (1994) Navarro, J. F., & White, S. D. M. 1994, MNRAS, 267, 401, doi: 10.1093/mnras/267.2.401
  * Nishimura et al. (2015) Nishimura, N., Takiwaki, T., & Thielemann, F.-K. 2015, ApJ, 810, 109, doi: 10.1088/0004-637X/810/2/109
  * Nomoto (1987) Nomoto, K. 1987, ApJ, 322, 206, doi: 10.1086/165716
  * Nomoto et al. (1997) Nomoto, K., Iwamoto, K., Nakasato, N., et al. 1997, Nuclear Physics A, 621, 467, doi: 10.1016/S0375-9474(97)00291-1
  * Nomoto et al. (2013) Nomoto, K., Kobayashi, C., & Tominaga, N. 2013, ARAA, 51, 457, doi: 10.1146/annurev-astro-082812-140956
  * Pagel (1997) Pagel, B. E. J. 1997, Nucleosynthesis and Chemical Evolution of Galaxies
  * Pakmor et al. (2011) Pakmor, R., Bauer, A., & Springel, V. 2011, MNRAS, 418, 1392, doi: 10.1111/j.1365-2966.2011.19591.x
  * Pilkington et al. (2012) Pilkington, K., Few, C. G., Gibson, B. K., et al. 2012, A&A, 540, A56, doi: 10.1051/0004-6361/201117466
  * Pillepich et al. (2018) Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077, doi: 10.1093/mnras/stx2656
  * Portinari et al. (1998) Portinari, L., Chiosi, C., & Bressan, A. 1998, A&A, 334, 505.  https://arxiv.org/abs/astro-ph/9711337
  * Prantzos et al. (2018) Prantzos, N., Abia, C., Limongi, M., Chieffi, A., & Cristallo, S. 2018, MNRAS, 476, 3432, doi: 10.1093/mnras/sty316
  * Prantzos et al. (1993) Prantzos, N., Casse, M., & Vangioni-Flam, E. 1993, ApJ, 403, 630, doi: 10.1086/172233
  * Reichert et al. (2021) Reichert, M., Obergaulinger, M., Eichler, M., Aloy, M. Á., & Arcones, A. 2021, MNRAS, 501, 5733, doi: 10.1093/mnras/stab029
  * Roederer et al. (2014) Roederer, I. U., Preston, G. W., Thompson, I. B., et al. 2014, AJ, 147, 136, doi: 10.1088/0004-6256/147/6/136
  * Romano et al. (2010) Romano, D., Karakas, A. I., Tosi, M., & Matteucci, F. 2010, A&A, 522, A32, doi: 10.1051/0004-6361/201014483
  * Romano et al. (2019) Romano, D., Matteucci, F., Zhang, Z.-Y., Ivison, R. J., & Ventura, P. 2019, MNRAS, 490, 2838, doi: 10.1093/mnras/stz2741
  * Rowan-Robinson et al. (2016) Rowan-Robinson, M., Oliver, S., Wang, L., et al. 2016, MNRAS, 461, 1100, doi: 10.1093/mnras/stw1169
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161, doi: 10.1086/145971
  * Scannapieco et al. (2012) Scannapieco, C., Wadepuhl, M., Parry, O. H., et al. 2012, MNRAS, 423, 1726, doi: 10.1111/j.1365-2966.2012.20993.x
  * Schaye et al. (2015) Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521, doi: 10.1093/mnras/stu2058
  * Seitenzahl et al. (2013) Seitenzahl, I. R., Cescutti, G., Röpke, F. K., Ruiter, A. J., & Pakmor, R. 2013, A&A, 559, L5, doi: 10.1051/0004-6361/201322599
  * Siegel et al. (2019) Siegel, D. M., Barnes, J., & Metzger, B. D. 2019, Nature, 569, 241, doi: 10.1038/s41586-019-1136-0
  * Smartt (2009) Smartt, S. J. 2009, ARAA, 47, 63, doi: 10.1146/annurev-astro-082708-101737
  * Sneden et al. (2016) Sneden, C., Cowan, J. J., Kobayashi, C., et al. 2016, ApJ, 817, 53, doi: 10.3847/0004-637X/817/1/53
  * Spite et al. (2005) Spite, M., Cayrel, R., Plez, B., et al. 2005, A&A, 430, 655, doi: 10.1051/0004-6361:20041274
  * Spitoni et al. (2019) Spitoni, E., Silva Aguirre, V., Matteucci, F., Calura, F., & Grisoni, V. 2019, A&A, 623, A60, doi: 10.1051/0004-6361/201834188
  * Springel (2010) Springel, V. 2010, MNRAS, 401, 791, doi: 10.1111/j.1365-2966.2009.15715.x
  * Springel et al. (2005) Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361, 776, doi: 10.1111/j.1365-2966.2005.09238.x
  * Springel & Hernquist (2003) Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289, doi: 10.1046/j.1365-8711.2003.06206.x
  * Springel et al. (2001) Springel, V., Yoshida, N., & White, S. D. M. 2001, New Astronomy, 6, 79, doi: 10.1016/S1384-1076(01)00042-2
  * Steinmetz & Mueller (1994) Steinmetz, M., & Mueller, E. 1994, A&A, 281, L97.  https://arxiv.org/abs/astro-ph/9312010
  * Sutherland & Dopita (1993) Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253, doi: 10.1086/191823
  * Taylor & Kobayashi (2014) Taylor, P., & Kobayashi, C. 2014, MNRAS, 442, 2751, doi: 10.1093/mnras/stu983
  * Taylor & Kobayashi (2015a) —. 2015a, MNRAS, 452, L59, doi: 10.1093/mnrasl/slv087
  * Taylor & Kobayashi (2015b) —. 2015b, MNRAS, 448, 1835, doi: 10.1093/mnras/stv139
  * Taylor & Kobayashi (2016) —. 2016, MNRAS, 463, 2465, doi: 10.1093/mnras/stw2157
  * Taylor & Kobayashi (2017) —. 2017, MNRAS, 471, 3856, doi: 10.1093/mnras/stx1860
  * Taylor et al. (2020) Taylor, P., Kobayashi, C., & Kewley, L. J. 2020, MNRAS, 496, 4433, doi: 10.1093/mnras/staa1904
  * Teyssier (2002) Teyssier, R. 2002, A&A, 385, 337, doi: 10.1051/0004-6361:20011817
  * Thomas et al. (2005) Thomas, D., Maraston, C., Bender, R., & Mendes de Oliveira, C. 2005, ApJ, 621, 673, doi: 10.1086/426932
  * Timmes et al. (1995) Timmes, F. X., Woosley, S. E., & Weaver, T. A. 1995, ApJS, 98, 617, doi: 10.1086/192172
  * Tinsley (1980) Tinsley, B. M. 1980, Fundamentals of Cosmic Physics, 5, 287
  * Tominaga (2009) Tominaga, N. 2009, ApJ, 690, 526, doi: 10.1088/0004-637X/690/1/526
  * Umeda & Nomoto (2003) Umeda, H., & Nomoto, K. 2003, Nature, 422, 871, doi: 10.1038/nature01571
  * van de Voort et al. (2022) van de Voort, F., Pakmor, R., Bieri, R., & Grand, R. J. J. 2022, MNRAS, 512, 5258, doi: 10.1093/mnras/stac710
  * van de Voort et al. (2020) van de Voort, F., Pakmor, R., Grand, R. J. J., et al. 2020, MNRAS, 494, 4867, doi: 10.1093/mnras/staa754
  * van Dokkum & Conroy (2010) van Dokkum, P. G., & Conroy, C. 2010, Nature, 468, 940, doi: 10.1038/nature09578
  * Ventura & D’Antona (2009) Ventura, P., & D’Antona, F. 2009, A&A, 499, 835, doi: 10.1051/0004-6361/200811139
  * Vincenzo & Kobayashi (2018a) Vincenzo, F., & Kobayashi, C. 2018a, MNRAS, 478, 155, doi: 10.1093/mnras/sty1047
  * Vincenzo & Kobayashi (2018b) —. 2018b, A&A, 610, L16, doi: 10.1051/0004-6361/201732395
  * Vincenzo & Kobayashi (2020) —. 2020, MNRAS, 496, 80, doi: 10.1093/mnras/staa1451
  * Vincenzo et al. (2019) Vincenzo, F., Kobayashi, C., & Yuan, T. 2019, MNRAS, 488, 4674, doi: 10.1093/mnras/stz2065
  * Wallner et al. (2021) Wallner, A., Froehlich, M. B., Hotchkis, M. A. C., et al. 2021, Science, 372, 742, doi: 10.1126/science.aax3972
  * Wanajo (2013) Wanajo, S. 2013, ApJ, Letter, 770, L22, doi: 10.1088/2041-8205/770/2/L22
  * Wanajo et al. (2013) Wanajo, S., Janka, H.-T., & Müller, B. 2013, ApJ, Letter, 767, L26, doi: 10.1088/2041-8205/767/2/L26
  * Wanajo et al. (2009) Wanajo, S., Nomoto, K., Janka, H. T., Kitaura, F. S., & Müller, B. 2009, ApJ, 695, 208, doi: 10.1088/0004-637X/695/1/208
  * Wanajo et al. (2014) Wanajo, S., Sekiguchi, Y., Nishimura, N., et al. 2014, ApJ, Letter, 789, L39, doi: 10.1088/2041-8205/789/2/L39
  * Wang et al. (2019) Wang, E. X., Taylor, P., Federrath, C., & Kobayashi, C. 2019, MNRAS, 483, 4640, doi: 10.1093/mnras/sty3491
  * Weinberger et al. (2017) Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291, doi: 10.1093/mnras/stw2944
  * Wiersma et al. (2009) Wiersma, R. P. C., Schaye, J., & Smith, B. D. 2009, MNRAS, 393, 99, doi: 10.1111/j.1365-2966.2008.14191.x
  * Winteler et al. (2012) Winteler, C., Käppeli, R., Perego, A., et al. 2012, ApJ, Letter, 750, L22, doi: 10.1088/2041-8205/750/1/L22
  * Woods et al. (2019) Woods, T. E., Agarwal, B., Bromm, V., et al. 2019, PASA, 36, e027, doi: 10.1017/pasa.2019.14
  * Woosley & Weaver (1995) Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181, doi: 10.1086/192237
  * Worthey et al. (1992) Worthey, G., Faber, S. M., & Gonzalez, J. J. 1992, ApJ, 398, 69, doi: 10.1086/171836
  * Yong et al. (2021) Yong, D., Kobayashi, C., Da Costa, G. S., et al. 2021, Nature, 595, 223, doi: 10.1038/s41586-021-03611-2
  * Zaritsky et al. (1994) Zaritsky, D., Kennicutt, Robert C., J., & Huchra, J. P. 1994, ApJ, 420, 87, doi: 10.1086/173544
  * Zhang et al. (2018) Zhang, Z.-Y., Romano, D., Ivison, R. J., Papadopoulos, P. P., & Matteucci, F. 2018, Nature, 558, 260, doi: 10.1038/s41586-018-0196-x
  * Zhao et al. (2016) Zhao, G., Mashonkina, L., Yan, H. L., et al. 2016, ApJ, 833, 225, doi: 10.3847/1538-4357/833/2/225