ar
X
iv
:1
01
2.
51
44
v2
[
as
tro
-p
h.G
A]
1
1 F
eb
20
11
To appear in the Astrophysical Journal
Chemodynamical Simulations of the Milky Way Galaxy
Chiaki KOBAYASHI1 and Naohito NAKASATO2,3,4
1 Research School of Astronomy & Astrophysics, The Australian National University, Cotter Rd.,
Weston ACT 2611, Australia; chiaki@mso.anu.edu.au
2 School of Computer Science and Engineering, The University of Aizu, Aizu-Wakamatsu City,
Fukushima, 965-8580, Japan
3 Center for Computational Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8577, Japan
4 Computational Astrophysics Laboratory, RIKEN, Wako, Saitama 351-0198, Japan
ABSTRACT
We present chemodynamical simulations of a Milky Way-type galaxy using a self-
consistent hydrodynamical code that includes supernova feedback and chemical enrich-
ment, and predict the spatial distribution of elements from Oxygen to Zinc. In the
simulated galaxy, the kinematical and chemical properties of the bulge, disk, and halo
are consistent with the observations. The bulge formed from the assembly of subgalax-
ies at z >∼ 3, and has higher [α/Fe] ratios because of the small contribution from Type Ia
Supernovae. The disk formed with a constant star formation over 13 Gyr, and shows a
decreasing trend of [α/Fe] and increasing trends of [(Na,Al,Cu,Mn)/Fe] against [Fe/H].
However, the thick disk stars tend to have higher [α/Fe] and lower [Mn/Fe] than thin
disk stars. We also predict the frequency distribution of elemental abundance ratios as
functions of time and location, which can be directly compared with galactic archeology
projects such as HERMES.
Subject headings: Galaxy: abundances — Galaxy: evolution — Galaxy: formation —
methods: numerical — stars: supernovae
1. Introduction
While the evolution of the dark matter is reasonably well understood, the evolution of the bary-
onic component is much less certain because of the complexity of the relevant physical processes,
such as star formation and feedback. With the commonly employed, schematic star formation
criteria alone, the predicted star formation rates (SFRs) are higher than what is compatible with
the observed luminosity density. Thus feedback mechanisms are in general invoked to reheat gas

– 2 –
and suppress star formation. We include the feedback from stellar winds, core-collapse supernovae
(normal Type II Supernovae (SNe II) and hypernovae (HNe)), and Type Ia Supernovae (SNe Ia)
in our hydrodynamical simulations. Supernovae inject not only thermal energy but also heavy ele-
ments into the interstellar medium (ISM), which can enhance star formation. Chemical enrichment
must be solved as well as energy feedback. Supernova feedback is also important for solving the
angular momentum problem (Steinmetz & Navarro 1999; Sommer-Larsen, Go¨lz, & Portinari 2003)
and the missing satellite problem (Moore et al. 1999), and for explaining the existence of heavy
elements in intracluster medium (Renzini et al. 1993) and intergalactic medium, and possibly the
mass-metallicity relation of galaxies (Kobayashi et al. 2007, hereafter K07).
Since different heavy elements are produced from different supernovae with different timescales,
elemental abundance ratios can provide independent information on “age” (e.g., Kobayashi & Nomoto
2009, hereafter KN09). Therefore, stars in a galaxy can be used as fossils where the chemical en-
richment history of the galaxy is imprinted, and this approach is called as galactic archeology. In
the Milky Way Galaxy, the detailed structures of kinematics and metallicities have been obtained
with SDSS (the Sloan Digital Sky Survey) and RAVE (the Radial Velocity Experiment), and will
be extended with Pan-STARRS and SkyMapper. In the next decade, high-resolution multi-object
spectroscopy (e.g., APOGEE with SDSS and HERMES on the AAT) and space astrometry mission
(e.g., GAIA) will provide 6D-kinematics and elemental abundances for a million stars in the Local
Group. In order to untangle the formation and evolution history of the galaxy from observational
data, a “realistic” model that includes star formation and chemical enrichment is required.
Although there are several simulations of the Milky Way-type galaxy with very high resolutions
(e.g., Diemand, Kuhlen, & Madau 2007), there are few studies with hydrodynamics and chemical
enrichment (e.g., Scannapieco et al. 2009). In most of the previous work, however, the adopted
assumptions are not appropriate for comparison with observations; the effect of mass-loss from
stars is ignored, the timescales of SNe Ia is given by only one parameter, the nucleosynthesis yields
of supernovae are outdated, and so on. Recently, several new findings have been obtained from the
collaboration between observations, stellar astrophysics, and galactic chemical evolution modeling.
From the observed light curves and spectra of nearby core-collapse supernovae, it has been shown
that there are energetic supernovae, HNe, ejecting more than ten times larger explosion energy
(E51 >∼ 10) and a certain amount of iron (e.g., Nomoto et al. 2006) as well as α elements (O, Mg,
Si, S, Ca, and Ti). Such HNe play a major role in reproducing the observed [Zn/Fe] trend in the
solar neighborhood, where [Zn/Fe] is ∼ 0 over −2 <∼ [Fe/H] <∼ 0 and possibly increases toward
lower [Fe/H] (Kobayashi et al. 2006, hereafter K06). From the observations of supernova rates
in various types of galaxies (Sullivan et al. 2006; Mannucci et al. 2006), it has been shown that
there is a young population of SNe Ia in addition to the old population that is usually found in
present-day elliptical galaxies. The metallicity effect on the occurrence of SNe Ia, which has been
proposed by Kobayashi et al. (1998, hereafter K98), plays a major role in reproducing the observed
[(α,Mn)/Fe]-[Fe/H] relations (KN09).
Including the up-to-date knowledge of chemical evolution of galaxies, we provide chemody-

– 3 –
namical simulations of a Milky Way-type galaxy from CDM initial conditions. In §2, we summarize
our simulation code and chemical enrichment sources. In §3, we show the differences between the
chemical properties of the disk and bulge, namely the metallicity distribution function, the age-
metallicity relation, and the [α/Fe]-[Fe/H] relation, and predict the [X/Fe]-[Fe/H] relation. We also
propose a new diagram of [α/Fe]-[Mn/Fe] as a sequence of the SN Ia contribution. §4 gives our
conclusions.
2. Model
2.1. Hydrodynamics
The details of our GRAPE-SPH code are described in Kobayashi (2004) and can be summarized
as follows.
i) The gravity is computed with the tree method using the special purpose computer GRAPE
(GRAvity PipE) system at the National Astronomical Observatory of Japan. For hydro-
dynamics, the Smoothed Particle Hydrodynamical (SPH) method is adopted, and the SPH
formulation is the almost same as Navarro & White (1993). The GRAPE-SPH code is highly
adaptive in space and time by means of individual smoothing lengths and individual timesteps,
and it has very high performance (a week for one simulation with N ∼ 50, 000).
ii) Radiative cooling is computed with the metal-dependent cooling functions, which are gen-
erated with the MAPPINGS III software (Sutherland & Dopita 1993) as a function of [Fe/H].
[O/Fe] is fixed with the observed [O/Fe]-[Fe/H] relation in the solar neighborhood.
iii) Our star formation criteria are the same as Katz (1992); 1) convergent, 2) cooling, and 3)
Jeans unstable. The SFR is determined from the Schmidt law; the star formation timescale is
proportional to the dynamical timescale (tsf ≡ 1c tdyn), where c = 0.1 is chosen from the size-
luminosity relation of elliptical galaxies (Kobayashi 2005). In addition, we find that c = 0.1
gives a better fit to the elemental abundance ratios in the Milky Way Galaxy than c = 0.02
(see Fig. 12). We also adopt the probability criterion (Katz 1992). In our simulations, the
SFR depends both on the local density and metallicity, which is different from other simplified
models such as one-zone models.
iv) If a gas particle satisfies the star formation criteria, a part of the mass of the gas particle
turns into a star particle. Since one star particle has the mass of 105−6M⊙, the star particle
is not a single star but an association of many stars. The masses of individual stars span
according to the initial mass function (IMF). We adopt the Salpeter IMF that is invariant
to time and metallicity with a slope x = 1.35 for 0.07 − 120M⊙.
v) For the feedback of energy and heavy elements, we do not adopt the instantaneous recy-
cling approximation. Via stellar winds, SNe II/HNe, and SNe Ia, thermal energy and heavy

– 4 –
elements are ejected from an evolved star particle as a function of time, and distributed to
the surrounding gas particles within a sphere of feedback radius rFB. The mass transfer is
weighted by the SPH kernel from the gas particle to the dying star particle. This feedback
scheme is better than that implementing a fixed number of feedback neighbor particles (NFB)
in galaxy-scale simulations; the metallicity of gas particles does not depend very much on the
local density if NFB is fixed. We find that rFB = 1 kpc gives a better fit to the elemental
abundance ratios than rFB = 0.5 and 3 kpc, or NFB = 50 and 100 (see Fig. 12). Note that
the diffusion between gas particles is supposed to enhance the mixing of heavy elements but
is not included in our model.
vi) The photometric evolution of one star particle is identical to the evolution of a simple
stellar population (SSP). SSP spectra are taken from Kodama & Arimoto (1997) as a func-
tion of age and metallicity. The photometric evolution of the galaxy is calculated from the
summation of the SSP spectra of star particles with various ages and metallicities.
2.2. Chemical enrichment
The scheme to include the detailed chemical enrichment into SPH simulations has been pro-
posed by Kobayashi (2004) and updated as follows.
Hypernovae (HNe) — The explosion mechanism of core-collapse supernovae is still uncer-
tain, although a few groups have presented feasible calculations in exploding core-collapse super-
novae (Kitaura, Janka & Hillebrandt 2006; Marek & Janka 2009; Bruenn et al. 2009). However,
the ejected explosion energy and 56Ni mass (which decays to 56Fe) can be directly estimated from
the observations, i.e., the light curve and spectra fitting of individual supernova. As a result, it is
found that many core-collapse supernovae (M ≥ 20M⊙) have more than ten times larger explosion
energy (E51 >∼ 10) and produce a significant amount of iron (Nomoto et al. 2006).
We adopt our nucleosynthesis yields of core-collapse supernovae from K06, depending on the
progenitor mass, metallicity (Z = 0−0.05), and the explosion energy (normal SNe II and HNe). The
progenitors mass range is 8− 50M⊙ but the metal ejection is only from 10− 50M⊙. Assuming that
a large fraction of supernovae with M ≥ 20M⊙ are HNe, the evolution of the elemental abundance
ratios from oxygen to zinc is in excellent agreement with observations in the solar neighborhood,
bulge, halo, and thick disk. In particular, to solve the [Zn/Fe] problem, the HN fraction ǫHN
has to be as large as ∼ 50% (K06). Note that Ti is underabundant at −3 <∼ [Fe/H]<∼ 0, which
will be solved with the 2D calculation of nucleosynthesis (Maeda & Nomoto 2003). The observed
increasing trend toward the lower metallicity may suggest that higher energy under the assumption
of inhomogeneous enrichment.
At high metallicity, since neutron-rich isotopes 66−70Zn are produced, ǫHN can be as small
as 1%. In this work, we assume the metal-dependent efficiency of HNe (M ≥ 20M⊙), ǫHN =

– 5 –
0.5, 0.5, 0.4, 0.01, and 0.01 for Z = 0.0.001, 0.004, 0.02, and 0.051, which gives better agreement
with the observed present HN rate (Podsiadlowski et al. 2004). Pair-instability supernovae, which
produce much more Fe, more [S/Fe], and less [Zn/Fe], should not contribute in the galactic chemical
evolution (Kobayashi, Tominaga, & Nomoto 2010b).
Type Ia Supernovae (SNe Ia) — The progenitors of the majority of SNe Ia are most likely
the Chandrasekhar (Ch) mass white dwarfs (WDs). For the evolution of accreting C+O WDs
toward the Ch mass, two scenarios have been proposed; one is the double-degenerate scenario, i.e.,
the merging of double C+O WDs with a combined mass surpassing the Ch mass limit. However,
it has been theoretically suggested that it leads to accretion-induced collapse rather than SNe Ia
(Nomoto & Kondo 1991), and the lifetimes are too short to reproduce the chemical evolution in
the solar neighborhood (KN09). The other is the single-degenerate (SD) scenario, i.e., the WD
mass grows by accretion of hydrogen-rich matter via mass transfer from a binary companion. The
mass accretion rate is very limited to trigger carbon deflagration (Nomoto 1982), but the allowed
parameter space of binary systems can be significantly increased by the WD wind effect if the
metallicity is higher than [Fe/H] ∼ −1 (K98). Two progenitor systems are found for red-giant
(RG) companions with initial masses of ∼ 1M⊙ and near main-sequence (MS) companions with
initial masses of ∼ 3M⊙ (Hachisu et al. 1999).
In our simulations, based on the SD scenario, the lifetime distribution function of SNe Ia is
calculated with Eq.[2] in KN09, taking into account the metallicity dependence of the WD winds
(K98) and the mass-stripping effect on the binary companion stars (KN09). The lifetimes of SNe
Ia are determined from the lifetimes of companion stars, and thus the lifetime distribution function
shows double peaks; one is for the MS+WD systems with timescales of ∼ 0.1 − 1 Gyr, which
are dominant in star-forming galaxies. The other is for the RG+WD systems with ∼ 1 − 20 Gyr
timescales, which are dominant in early-type galaxies. Although the metallicity effect has not been
confirmed from the supernova survey yet, it is more strongly required in the presence of the young
population of SNe Ia, to be consistent with the chemical evolution of the Milky Way Galaxy. The
absolute rates of each component depend on two binary parameters, which are determined by the
requirement to meet the observed metallicity distribution function and [α/Fe]-[Fe/H] relations in
the solar neighborhood, and the observed supernova rates in various types of galaxies. This SN
Ia model gives much better agreement with the observed elemental abundance ratios (from O to
Zn) than models with the double-degenerate scenario or metal-independent exponential lifetime
functions (KN09).
For SNe Ia, we take the nucleosynthesis yields from Nomoto et al. (1997b), where the metal-
licity dependence is not included. Ni is overproduced at [Fe/H] >∼ − 1, which will be solved
by tuning the propagation speed of the burning front and the central density of the white dwarf
(Iwamoto et al. 1999).
1ǫHN = 0.5 is adopted independent of metallicity in K06, K07, and KN09.

– 6 –
Stellar winds — Needless to say, the envelope mass and pre-existing heavy elements are
returned by stellar winds from all dying stars. In addition, Asymptotic Giant Branch (AGB) stars
with initial masses between about 0.8−8M⊙ (depending on metallicity) produce light elements such
as C and N, while the contribution of heavier elements are negligible in the galactic chemical evo-
lution (Kobayashi 2010; Kobayashi, Karakas, & Umeda 2010a). Rotating massive stars (>∼ 40M⊙)
could also produce C and N (e.g., Meynet & Maeder 2002), but are not included.
2.3. Initial Condition
The initial condition is generated in the same way as in Nakasato & Nomoto (2003) and
Kobayashi (2004, 2005). The CDM initial fluctuation is generated by the COSMICS package
(Bertschinger 1995). The cosmological parameters are set to be H0 = 70 km s−1 Mpc−1, Ωm = 0.3,
ΩΛ = 0.7, Ωb = 0.04, n = 1, and σ8 = 0.9. We generate a periodic boundary condition with the
lattice size 8 Mpc having a top-hat perturbation of amplitude 3σ in radius 1.4 Mpc at the starting
redshift of z = 24. Then, particles in a spherical region with a comoving radius that contains
∼ 1012M⊙ are picked up. Finally, because the simulated field is not large enough to generate tidal
torque, the initial angular momentum is given to the system in rigid rotation with a constant spin
parameter as large as λ ≡ J |E|1/2G−1M−5/2 = 0.1 (Nakasato & Nomoto 2003; Kobayashi et al.
2003). Such a large value is required to reproduce disk galaxies, although rare in cosmological
simulations (e.g., Warren et al. 1992).
It seems not easy to form the Galactic disk in the CDM Universe as suggested by previous works
(e.g., Toth & Ostriker 1992; Kauffmann & White 1993, but see Koda et al. 2009; Hopkins et al.
2009). To form a disk galaxy, it is necessary that the system does not undergo the major mergers
with mass ratios larger than 0.2 and at higher redshifts than z = 2, which destroy the disk structure
and the radial metallicity gradient (Kobayashi 2004). Even in some cases without major mergers,
disk structures cannot form (e.g., Scannapieco et al. 2009), although this argument depends on the
modeling of feedback and numerical resolution. We have run simulations with 150 initial conditions
Table 1. Initial and Final Quantities. (1) comoving radius, (2) initial number of all particles, (3)
total mass, (4) mass of an SPH particle, (5) final number of all particles, (6) final baryon fraction
within 100 kpc, and (7) final stellar fraction within 100 kpc.
(1) (2) (3) (4) (5) (6) (7)
Model rinit Ntot,init Mtot Msph Ntot,final b f∗
Standard 1.5 Mpc 55,534 0.8 × 1012M⊙ 3.8× 106M⊙ 180,293 0.14 0.23
Wider region 2 Mpc 117,525 1.7 × 1012M⊙ 3.8× 106M⊙ 316,991 0.15 0.25
Higher resolution 1.5 Mpc 209,133 0.8 × 1012M⊙ 1.0× 106M⊙ 445,259 0.12 0.29

– 7 –
until z = 1.28 including gas dynamics and basic chemical enrichment (Fig. 1). The disk structure
with a scale length of 3 − 4 kpc is seen in 48 runs, and Milky Way type galaxies form only in 5
runs, from which we select one realization (143rd in Fig. 1) for detailed simulations.
In order to check the resolution and boundary effects, we consider three cases; the total number
of particles, the total mass (the half for gas and the rest for dark matter), and the mass resolution
for each case are summarized in Table 1. The initial conditions with the wider region and higher
resolution are generated with the same initial perturbation as in the initial condition of the standard
model. Because of star formation, the particle number increases by a factor of 2− 3 at z = 0. The
final baryon and gas fractions of the galaxy are also summarized in Table 1 and do not vary very
much among these cases. For the higher resolution run, the calculation is stopped at z = 0.12. The
gravitational softening length is set to 0.5 and 0.32 kpc for the low and high resolutions, respectively.
We obtain qualitatively similar results for all three cases, and consequently, we mainly show the
results with the wider initial condition in Figs. 2-11 and 15, and with the standard initial condition
in Figs. 12, 13, 14, and 16-20.
3. Results
3.1. Star Formation History
Figure 2 shows the time evolution of the projected density of dark matter, gas, V-band lumi-
nosity, and luminosity-weighted stellar metallicity. After the start of the simulation, the system
expands according to the Hubble flow. The CDM initial fluctuations grow into the structures of
nodes and filaments, and small collapsed halos are realized both in dark matter and gas. In the
halos, the gas is allowed to cool radiatively, and star formation takes place since z ∼ 15. Accord-
ing to the hierarchical clustering of dark halos, subgalaxies merge to form large galaxies, which
induces the initial starburst. Under the CDM picture, any galaxy forms through the successive
merging of subgalaxies with various masses. In this simulated galaxy, the bulge is formed by the
initial starburst that is induced by the assembly of gas-rich sub-galaxies with stellar masses of
∼ 5 − 10 × 109M⊙ and gas fractions of 0.2 − 0.4 at z >∼ 3. Because of the angular momentum,
the gas accretes onto the plane forming a rotationally supported disk that grows from inside out.
In the disk, star formation takes place with a longer timescale, which is maintained not by the
slow gas accretion, but by the self-regulation due to supernova feedback. Many satellite galaxies
successively come in and disrupt, but there is no major merger event after z ∼ 2, which is necessary
to retain the disk structure. Metallicity gradients, increasing toward higher density regions, are
generated both in the gas phase and stars from z ∼ 5 onward.
Figure 3 shows the B-band (upper panels) and K-band (lower panels) luminosity maps at the
present-day for the face-on (left panels) and edge-on views (right panels). In the B-band, the disk
component is visible because of the on-going star formation, while the bulge component is dominant
in the K-band. The galaxy center and the disk plane are determined from the center of gravity

– 8 –
and the aliment of the angular momentum vector, respectively. In the following, r and z denote
the radius on the disk plane and the height from the disk plane, respectively.
The bulge has a de Vaucouleurs surface brightness profile with an effective radius of ∼ 1.5 kpc,
and the disk has an exponential profile with a scale length of ∼ 5 kpc. In the surface brightness
profile, an excess is seen at r ∼ 12 kpc, which corresponds to the remnant of satellite galaxies. The
bulge in the simulated galaxy seems to be larger than the Galactic bulge, and these properties may
be more consistent with M31 (e.g., Klypin, Zhao, & Somerville 2002). Note that it is not feasible
to resolve the central blackhole, nuclear bulge, and bar in our simulations. The total masses inside
1, 10, and 100 kpc are 9.4 × 109M⊙, 2.2 × 1011M⊙, and 1.1 × 1012M⊙, respectively. The baryon
fraction is ∼ 0.65 at r < 2 kpc, and then decreases to ∼ 0.2 at r > 8 kpc. The rotation velocity of
the disk rapidly increases to reach ∼ 230 km/s at r ∼ 3 kpc, stays constant until r <∼ 20 kpc, and
then gradually decreases until r ∼ 40 kpc. This is roughly consistent with the observed rotation
curve of the Milky Way Galaxy, which shows a plateau at the local standard of rest (180 − 250
km/s) for 3 − 6 <∼ r <∼ 20 − 60 kpc (Merrifield 1992; Honma et al. 2007; Levine, Heiles, & Blitz
2008).
In the following, we define the three major components simply from the location of the stars
at the present-day: the radius of 7.5 ≤ r ≤ 8.5 kpc and the height of |z| ≤ 0.5 kpc for the solar
neighborhood, r ≤ 1 kpc for the bulge, and 5 ≤ r ≤ 10 kpc for the halo. Thick disk stars are defined
from the kinematics using the Toomre diagram in the cylindrical space motion (VΦ, VΘ, VZ). The
number of particles significantly drops around the velocity dispersion σ ≡
√
V 2Φ + V 2Θ + V 2Z ∼ 150
km/s and σ ∼ 250 km/s. Thus, we define thick disk particles in the range of σ = 160 − 260
km/s. This range is different from the observations (Bensby et al. 2003; Ruchti et al. 2010), but
gives a consistent mass fraction of the thick disk in the simulated galaxy. With this kinematical
definition, the thick disk is smaller and thicker than the thin disk. The surface density profiles
(top panel) and density profiles (bottom panel) at r = 7− 9 kpc for the thin disk (solid lines) and
thick disk (dot-dashed lines) are shown in Figure 4. The scale radius is 4.9 and 4.1 kpc, the scale
height is 1.9 and 3.2 for the thin and thick disks, respectively. These are larger than the 2MASS
(Two Micron All Sky Survey) observations (Lo´pez-Corredoira et al. 2002), which may be due to
the lack of resolution. The mass, mass-to-light ratio, age, metallicity, and abundance ratio of each
component are summarized in Table 2.
Figure 5 shows the radial and vertical metallicity gradients of stars at present (top panel)
and at t = 8 Gyr (bottom panel). The radial gradient of stars becomes flatter at higher distances
from the plane, which seems to be consistent with the SDSS observations (Morrison 2010). At
|z| = 0− 0.5 kpc (solid line), the stellar radial gradient is ∆[M/H]∗/∆r = −0.025 dex/kpc. This is
flatter than the radial gradient of the ISM (dotted line), which is ∆[M/H]g/∆r = −0.049 dex/kpc
at |z| = 0 − 0.5 kpc. The metallicity radial gradients are steeper at higher redshifts. At t = 8
Gyr, ∆[M/H]∗/∆r = −0.062 dex/kpc and ∆[M/H]g/∆r = −0.094 dex/kpc for stars and the ISM,
respectively. The time evolution of the ISM gradient is consistent with the observations with
planetary nebulae (Macial, Lago, & Costa 2006).

– 9 –
The resultant star formation histories (Fig. 6) are different for different components. The
initial starburst is induced by the assembly, and most of stars in the bulge (dashed line) and halo
(dotted line) are formed at this stage. In the disk (solid line), the SFR is regulated mainly by
supernova feedback. Figure 7 shows the cumulative functions of ages of stars in each component
in the present-day galaxy. Most of stars in the present bulge (dashed line) have formed in the first
2 Gyr, 80% of bulge stars are older than 10 Gyr, and 90% of bulge stars are older than 8 Gyr. In
the disk (solid line), star formation takes place continuously, and 50% of solar-neighborhood stars
are younger than 8 Gyr. In the halo (dotted line), most of stars are as old as bulge, although the
numerical resolution of this simulation is not enough to discuss the halo in detail. However, there
does not seem to be very much difference in the SFR between the inner and outer halos. These
star formation histories can be seen for other initial conditions that form disk galaxies at present
in our simulations, and also for different initial conditions e.g., those in Scannapieco et al. (2009)
and Sa´nchez-Bla´zquez et al. (2009).
3.2. Inhomogeneous Chemical Enrichment
In chemodynamical simulations, the ISM is not homogeneous at any time. First of all, there
is a local variation in star formation and metal flow by the inflow and outflow of the ISM. Second,
the mixing of heavy elements is caused by supernova feedback. With our feedback scheme (§2.1),
star particles obtain heavy elements from the gas particles from which the stars form. Gas particles
obtain heavy elements only when they pass through within a feedback radius of dying star particles.
The amount of heavy elements that the gas particles obtain is physically uncertain, and depend
on the modeling of feedback. Third, the elemental abundance ratios of the metals produced by
supernovae depend on the metallicity and mass of progenitors (§2.2), i.e., the metallicity and age
of star particles.
In our simulations, the metallicity of the first enriched stars reaches [Fe/H] ∼ −3. At later
times, the star forming region becomes denser, and both metal richer and poorer stars than [Fe/H]
∼ −3 appear. Different from one-zone chemical evolution models, the following phenomena occur
in the case of inhomogeneous enrichment: i) The age-metallicity relation is weak. ii) SNe Ia can
affect the elemental abundance ratios at [Fe/H] <∼ − 1 even with the metallicity inhibition of SNe
Ia. iii) The scatter of elemental abundance ratios becomes large if the supernova yield depends on
progenitor metallicity. The details are described in the following sections. As discussed in §2.1, we
run simulations with several feedback parameters and schemes, and show the best results in this
paper.

–
10
–
Table 2. Mean Values in the Solar Neighborhood, Bulge, Halo, and Thick Disk in the Simulated Milky Way-type Galaxy: (1)
stellar mass, (2-4) B, V, and K-band mass-to-light ratios, (5-6) mass and luminosity weighted age of stars, (7-8) mass and
luminosity weighted metallicity of stars, (9) metallicity of gas, (10-11) mass and luminosity weighted [O/Fe], (12) [O/Fe] of
gas, (13) gas fraction, and (14) baryon fraction.
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14)
Component M∗ [M⊙] M/LB M/LV M/LK age [Gyr] ageL [Gyr] [M/H] [M/H]L [M/H]g [O/Fe] [O/Fe]L [O/Fe]g fg b
solar neighborhood 8.72× 108 1.90 2.22 0.95 7.07 3.27 -0.19 -0.20 -0.23 -0.07 -0.16 -0.21 0.42 0.63
bulge 4.53× 109 6.92 6.15 1.51 11.17 9.03 -0.01 -0.03 0.14 0.01 -0.10 -0.32 0.22 0.62
halo 1.26× 109 7.92 6.28 1.69 12.02 12.07 -0.27 -0.38 -1.63 0.24 0.25 -0.85 0.02 0.06
thick disk 1.78× 108 5.57 4.53 1.27 8.61 7.49 -0.23 -0.30 - -0.04 -0.08 - - 0.70

– 11 –
3.3. Age-Metallicity Relations
The chemical enrichment timescale is also different for the different components. Figure 8
shows the age-metallicity relations in the present-day galaxy. In the solar neighborhood (upper
panel), the average metallicity reaches [Fe/H] ∼ 0 at t ∼ 2 Gyr and does not show strong evolution
for t >∼ 2 Gyr. The scatter in metallicity at given time is caused by the inhomogeneity of chemical
enrichment in our chemodynamical model (§3.2). As a result, both the average and scatter are
in good agreement with the observations (dots) in spite of the uncertainties in the observational
estimates of the ages. At t <∼ 5 Gyr, a week trend is seen in our simulation, which seems to be
consistent with the RAVE observations (Freeman et al. 2010).
In the bulge (lower panel), star formation takes place more quickly, and thus the chemical
enrichment timescale is much shorter than in the disk. The age-metallicity relation shows a more
rapid increase than in the disk. The maximum metallicity reaches super solar ([Fe/H] ∼ 1) at t ∼ 2
Gyr. Although the SFR becomes small after ∼ 5 Gyr, a few stars form at >∼ 5 Gyr. These have
super-solar metallicity in general and the average metallicity does not show time evolution.
3.4. [α/Fe]-[Fe/H] Relations
The difference in the chemical enrichment timescales results in a difference in the elemen-
tal abundance ratios, since different elements are produced by different supernovae with different
timescales. The best known clock is the α-elements to iron ratio [α/Fe] since SNe Ia produce more
iron than α elements with longer timescales than SNe II. Nevertheless it should be noted that
low-mass SNe II also provide relatively low [α/Fe] because of their smaller envelope mass compared
to more massive SNe II. This mass dependence is important in a system with a low SFR such as
dwarf spheroidal galaxies.
Figure 9 shows the cumulative functions of [O/Fe] of stars in each component in the present-
day galaxy. In the solar neighborhood (solid line), 80% of stars have [O/Fe] < 0.3, which indicates a
significant contribution from SNe Ia. In the bulge (dashed line) and halo (dotted line), the chemical
enrichment timescale is so short that 60% and 80% of stars have [O/Fe] > 0.3, respectively, although
there is a significant number of low [α/Fe] stars. The majority of low [α/Fe] stars also have high
[Mn/Fe] because of the SN Ia contribution. They are formed at t >∼ 1 Gyr and have non-rotating
(vΘ/σ <∼ 0.3) kinematics. A small fraction of low [α/Fe] stars, however, have low [Mn/Fe] and are
formed at t <∼ 1 Gyr. They represent the effect of low-mass SNe II in the inhomogeneous enrichment
(the third reason in §3.2). Such low [α/Fe] stars seem to be observed in SEGUE (Beers 2010). The
fraction of these stars is important for discussions of the formation history of the Galactic halo. In
order to distinguish these two enrichment sources, elemental abundance ratios of iron-peak elements
such as Mn are critically important.
Figure 10 shows the [O/Fe]-[Fe/H] relations, and the other α-elements show the same trends.

– 12 –
In the solar neighborhood (upper panel), at the beginning, only SNe II contribute, and [α/Fe] shows
a plateau ([α/Fe]∼ 0.5). Around [Fe/H] ∼ −1, SNe Ia start to occur, which produce more iron
than α elements. This delayed enrichment of SNe Ia causes the decrease in [α/Fe] with increasing
[Fe/H]. This trend is in great agreement with the observations (dots). Theoretically, the decreasing
point of [α/Fe] has been used to put a constraint on the lifetimes of SNe Ia. Kobayashi et al. (1998)
and Kobayashi & Nomoto (2009) showed that [α/Fe] decreases too early in the double-degenerate
scenario and in the single degenerate without metallicity effect since the average lifetime is too
short. The slope of ∆[α/Fe]/∆[Fe/H] depends on the binary parameters in our SN Ia model, but
we should note that it is consistent with the observed supernova rates in various types of galaxies
(§2.2).
A significant scatter is seen in Figure 10, and more clearly, the distribution functions of [O/Fe]
at given [Fe/H] in the solar neighborhood are shown in Figure 11. At [Fe/H] <∼ −1, [O/Fe] is almost
constant, and slightly larger for higher metallicity (the peak [O/Fe] = 0.6, 0.5, and 0.4 at [Fe/H]
= −3,−2, and −1, respectively) because of the mass dependence of the yields of SNe II. The scatter
of [O/Fe] is caused by the inhomogeneity of chemical enrichment in our chemodynamical model.
At [Fe/H] >∼ − 1, the majority of [O/Fe] is lower than 0.2 and the peak [O/Fe] is −0.15. Around
[Fe/H] ∼ −1, the scatter looks a bit larger than observed (Fig. 10). This may be because the
mixing of heavy elements among gas particles is not included in our model (§2.1). At [Fe/H] <∼ −1,
the scatter is caused from the following two reasons. First, although SNe Ia cannot occur at [Fe/H]
≤ −1.1 in our SN Ia progenitor scenario, SNe Ia can contribute in the abundances of stars at [Fe/H]
<∼ − 1. As discussed in §3.2, in the inhomogeneous enrichment, it is natural that the enrichment
sources have [Fe/H] > −1.1 but the stars forming from the ejecta have [Fe/H] <∼ − 1 because of the
large amount of hydrogen in protostellar clouds. In this case, [Mn/Fe] is as high as ∼ 0 (see §3.6 for
more discussion). Secondly, there is an intrinsic variation in [α/Fe] of supernova yields depending
on the progenitor mass and metallicity, of which effect is included in our chemodynamical model.
Low-mass SNe II (10 − 13M⊙) provide smaller [α/Fe] than massive SNe II. In this case, [Mn/Fe]
is lower than 0. The [α/Fe] scatter in the simulation could be larger, since there is also a variation
depending on the explosion energy and the remnant mass (neutron star and blackhole) for SNe II
and HNe.
From the comparison of the [α/Fe] scatter between simulations and observations, it is possible
to put constraints on the unsolved physics of supernova explosion and mixing of the ISM. Figure
12 shows the distribution functions of [O/Fe] for stars in the solar neighborhood at −0.5 ≤ [Fe/H]
≤ −0.3 for the standard model (solid line) and the models with a longer star formation timescale
(c = 0.02, dashed line), a smaller feedback radius (rFB = 0.5 kpc, dotted line), and a fixed feedback
number (NFB = 50, dot-dashed line). The smaller initial condition (Tab. 1) is used here. The
standard model gives much better agreement with the peak value of the observations (arrow) than
the other models.
Using stochastic chemical evolution model, Argast et al. (2002) claimed that the scatter of
[α/Fe] caused by supernova yields is too large compared to observations and an efficient mixing pro-

– 13 –
cess is required. This is because the nucleosynthesis yields adopted in that paper (Tsujimoto et al.
1995; Thielemann et al. 1996; Nomoto et al. 1997a) were different from currently accepted values.
The low [α/Fe] values resulted from the larger amount of iron assumed for low-mass SNe II (0.15M⊙
for 13 − 15M⊙). Now the iron mass is estimated from the light curve and spectral fitting, and is
0.07M⊙ for SNe II with < 20M⊙. The high [α/Fe] values are reduced by including hypernovae.
Figure 13 shows the distribution functions at −4 ≤ [Fe/H] ≤ −1.9 for the models with (solid line)
and without (dashed line) hypernova feedback using the smaller initial condition. With hypernovae,
the scatter around [O/Fe] ∼ 0.5 is as small as the observations (Cayrel et al. 2004).
In the bulge (lower panel of Fig. 10), the [α/Fe]-[Fe/H] relation is very much different. The
chemical enrichment timescale is so short that the metallicity reaches super solar before SNe Ia
contribute. Thus, the [α/Fe] plateau continues to [Fe/H] ∼ +0.3. This is roughly consistent
with the observations (Fulbright, McWilliam, & Rich 2007; Lecureur et al. 2007), except for the O
observation in Fulbright, McWilliam, & Rich (2007, see §3.6 for more discussion). In this simulated
galaxy, some new stars are still forming in the bulge, which in general have super solar metallicity
and low [α/Fe] because of the large contribution from SNe Ia. Although a small fraction of stars
with the age of ∼ 3 Gyr (t ∼ 10.3 Gyr) do have high [α/Fe] ([O/Fe] = 0.3), which is caused by
the inhomogeneous enrichment, stars younger than 1 Gyr have −0.5 ≤ [O/Fe] ≤ 0 and 0 ≤ [Fe/H]
≤ 0.8.
Cunha et al. (2007), however, showed several young metal-rich stars with high [α/Fe] in the
Galactic center, which suggests additional physical processes such as (i) local enrichment by the
stellar winds of dying stars or SNe II, (ii) intrinsic variation of supernovae depending on energy
and mixing, or (iii) the selective mass-loss of the gas enriched by SNe Ia. The bulge wind might be
driven by SNe Ia, and most of Fe enriched gas might be selectively blown away before these young
stars are formed. If the observed scatter in [O/H], [O/Fe], and [Ca/Fe] are real, then the scenarios
(i) and (ii) may be preferable.
3.5. Metallicity Distribution Functions
The metallicity distribution function (MDF) of G-dwarf stars is a stringent constraint on
chemodynamical models. Figure 14 shows the present MDF for three different models with the
smaller initial condition. The model with the feedback from the UV background radiation (UVB)
and hypernova (solid line) gives much better agreement with the observations in the solar neigh-
borhood. However, the number of metal-poor stars is still larger than observed, in other words,
there is the G-dwarf problem. In one-zone chemical evolution models, this problem is easily solved
by assuming a long timescale of gas infall (e.g., Tinsley 1980), but we cannot control the infall
by hand in chemodynamical simulations. The mass accretion timescale is mainly determined from
the hierarchical clustering of dark halos. Although hypernova feedback could help in generating
outflow, the mass accretion history is not affected very much. In our models, 77% of dark matter
(45% of total mass) is accreted inside 10 kpc in the first 1 Gyr, then the rest is gradually accreted.

– 14 –
There is a room to change the star formation and chemical enrichment timescales. With the UVB,
the initial starburst is suppressed by a factor of 1.5, which results in the longer duration of star
formation. With HNe, the SFR is a little suppressed at t >∼ 3 Gyr, and the metal ejection is
significantly changed.
Figure 15 shows the present MDF in the different components of the simulated galaxy for the
best model with the wider initial condition. In the solar neighborhood (solid line), the MDF is
roughly consistent with the previous observations (Edvardsson et al. 1993; Wyse & Gilmore 1995).
However, the Geneva-Copenhagen survey showed a narrower MDF (Holmberg et al. 2007, updated
from Nordstro¨m et al. 2004). If this is valid, the G-dwarf problem becomes more serious. In the
bulge (dashed line), McWilliam & Rich (1994) and Sadler et al. (1996) showed a broad MDF, from
which Nakasato & Nomoto (2003) discussed the metal-rich population in the bulge. Also in our
simulations, the metal-rich peak at [Fe/H] ∼ 0.4 is caused by the delayed enrichment of SNe Ia.
However, Zoccali et al. (2008) showed a narrower MDF with a large spectroscopic sample. There
seems to be the G-dwarf problem also for the Galactic bulge. The lack of metal-rich tail in the
MDF may suggest that the star formation is truncated in the bulge, possibly by an SNe Ia driven
bulge wind.
The other solutions of the G-dwarf problem are, as summarized in Tinsley (1980), pre-enrichment
and/or variable IMF. If the gas is enriched by the first stars before it accretes onto the Galaxy, the
number of metal-poor stars can be small. A bulge wind driven by SNe Ia could also work as pre-
enrichment in the disk. If the IMF at low-metallicity is top-heavy as suggested by the simulations of
primordial star formation (Bromm & Larson 2004), larger amounts of metals and heating sources
are provided. If no low-mass stars can be formed at low-metallicity, the number of metal-poor
G-dwarf stars obviously becomes zero. We should note, however, that the observed abundance
patterns of the extremely metal-poor (EMP) stars rule out large contributions from pair instability
supernovae that result from ∼ 140− 270M⊙ stars. Stars with <∼ 100M⊙ and >∼ 300M⊙ are usually
supposed to collapse to blackholes, but may be able to explode as core-collapse supernovae, which
do not conflict with the observations (Ohkubo et al. 2006). Or, simply, simulations with higher
resolution may solve the G-dwarf problem.
3.6. [X/Fe]-[Fe/H] Diagrams
Finally, using our chemodynamical simulation, we predict the frequency distributions of the
elements from O to Zn as a function of time and location. Figures 16-18 show the mass density of
stars in the [X/Fe]-[Fe/H] diagrams at present. These results can be statistically compared with
the next generation of observations from high resolution and multi-object spectrographs such as
HERMES. These instruments enable us not only to study the formation and evolution history of
the Galaxy but also to update stellar evolution and supernova physics.
Figure 16 is for the solar neighborhood. The dots are for the observations of individual stars,

– 15 –
and they almost perfectly overlap with the contours of our simulation. In order to minimize the
systematic error among various methods of observational data analysis, we consider only data from
Reddy et al. (2003, hereafter R03), Reddy, Lambert, & Prieto (2006), Reddy & Lambert (2008) at
[Fe/H] <∼ − 1, Fulbright (2000) at −2 <∼ [Fe/H] <∼ − 1, and Cayrel et al. (2004, hereafter C04)
and Honda et al. (2004, hereafter H04) at [Fe/H] <∼ − 2. As in K06, the 3D corrections of −0.24
for O in C04, and the NLTE corrections of +0.1 and +0.5, respectively, for Mg and Al in H04 are
applied. For Mg, Na, Al, and K, C04’s data are updated by Andrievsky et al. (2007, 2008, 2010,
hereafter A10) with the NLTE corrections. For O, S, Cu, and Zn, we employ other data sources as
noted below.
• α elements — O, Mg, Si, S, and Ca show the same relation; a plateau at [Fe/H] <∼ −1, and then
a decrease due to the delayed enrichment of SNe Ia (see KN09 for the SN Ia progenitor models).
For O and S, an increasing trend toward lower [Fe/H] reported by Israelian et al. (1998) and
Israelian & Reboro (2001), respectively, is affected by the NLTE and 3D effects, and is not
seen in our simulation. For O, there is an offset between the lines used in the abundance
analysis. To compare with the observations at low metallicity, we plot Bensby et al. (2004)’s
data with the forbidden line [OI] 6300A˚ since an empirical relation using [OI] and OI triplet at
7775A˚ is used in R03. For S, we plot data from Takada-Hidai et al. (2005) and Nissen et al.
(2007) at low metallicity, and Chen et al. (2002) at high metallicity, where more lines are
used than in R03. In R03, the scatter of [S/Fe] is large and the decreasing trend toward low
metallicity is not clearly seen. At [Fe/H] >∼ − 1, [Si/Fe] may be slightly larger, and [Ca/Fe]
may be slightly smaller than the observations, which are due to the SN Ia yields.
• Odd-Z elements — Na, Al, and Cu show a decreasing trend toward lower metallicity since the
nucleosynthesis yields of these elements strongly depend on the metallicity of progenitor stars
(see K06 for the details). Na and Al also show a decreasing trend toward higher metallicity
due to SNe Ia, which is shallower than that of α elements. [Na/Fe] at [Fe/H] >∼ − 1 may
be slightly larger than observed. The large scatter is caused by the metallicity dependence
of yields. As discussed in §3.2, in our chemodynamical simulations, there is only a week
correlation between the metallicities of supernova progenitors and those of the surrounding
ISM. If the ejected metals are diluted, the ISM metallicity and the metallicity of stars that
form from the ISM become lower than the progenitor metallicity. The scatter starts increasing
at [Fe/H] ∼ −2 in this simulation, which looks consistent with the observations. Observations
of an unbiased sample with the NLTE analysis will put a constraint on the mixing process in
the chemodynamical simulations.
• Iron-peak elements — The scatter of Cr and Co is very small since there is not very much
variety in the SN II yields and [(Cr, Co)/Fe] is almost ∼ 0 for SNe Ia. For Cr, since a
dependence on the temperature of the observed stars has been reported, we plot the Cr II
abundance from H04 for low metallicity, which is consistent with our simulation. For Co, an
increasing trend toward lower metallicity is seen in the observations, which is not realized

– 16 –
in the simulation. This is due to the input of supernova yields, and cannot be solved with
different star formation histories or initial conditions. This could be solved if we adopt higher
energy and/or higher fraction of HNe where more Co and Zn are synthesized.
• Zinc — Zn is one of the most important elements for supernova physics. [Zn/Fe] is about ∼ 0
for a wide range of metallicities, which can only be generated by such large contribution of
HNe. In detail, there is a small oscillating trend; [Zn/Fe] is 0 at [Fe/H] ∼ 0, increases to be 0.2
at [Fe/H] ∼ −0.5, decreases to be 0 at [Fe/H] ∼ −2, then increases toward lower metallicity.
This is characteristic in our SN Ia model (KN09), and is consistent with the observations
as found by Saito et al. (2009). Theoretically, Zn production depends on many parameters;
64Zn is synthesized in the deepest region of HNe, while neutron-rich isotopes of zinc 66−70Zn
are produced by neutron-capture processes, which are larger for higher metallicity massive
SNe II. The scatter at −1.5 <∼ [Fe/H] <∼ − 0.5 is larger than observed, which is caused by the
dependence of the Zn yield on the metallicity and mass. The increasing trend toward lower
metallicity in the observations (Primas et al. 2000; Nissen et al. 2007) could be generated
with the same effects as for Co.
• Manganese — Mn is a characteristic element and is more produced by SNe Ia relative to
iron. From [Fe/H] ∼ −1, [Mn/Fe] shows an increasing trend toward higher metallicity, which
is caused by the delayed enrichment of SNe Ia. Mn is also an odd-Z element, and the Mn
yield depends on the metallicity both for SNe II and SNe Ia. Feltzing, Fohlman, & Bensby
(2007) showed a bit steeper slope at [Fe/H] > 0, which would be generated by the metallicity
dependence of SN Ia yields. At [Fe/H] <∼ − 2, the scatter in the simulation is a bit smaller
than observed. Since SNe Ia do not contribute at such low metallicity, the [Mn/Fe] scatter, if
it is real, should stem from the SN II yields depending on the remnant mass and the electron
fraction Ye that is related to the mechanism of supernova explosion. In the supernova ejecta,
Mn and Cr are synthesized in the relatively outer regions compared to Cr, Zn and the majority
of Fe. Therefore, relatively high [(Mn, Cr)/Fe] and low [(Co, Zn)/Fe] are obtained with
faint SNe with relatively large remnant masses (blackholes; Kobayashi, Tominaga, & Nomoto
2010b).
In our figures, we take the solar abundance from Anders & Grevesse (1989), which is different
for some elements from the solar abundance updated by Asplund et al. (2010). This means the
normalization should be changed both for observations and simulations. When the observational
data are analyzed with 3D and NLTE effects, the agreement in Figures 16-18 may be or may not
be modified. For example, one of the NLTE models shows [Mn/Fe] ∼ 0, [Cr/Fe] ∼ 0, and steeper
increase in [Co/Fe] toward lower metallicity (Bergemann & Gehren 2008, 2010, Bergemann et al.
2010).
Figure 17 is for the bulge, where the observational data are taken from Fulbright, McWilliam, & Rich
(2007) and Zoccali et al. (2008). As shown in Fig. 10, the α-element plateau continues to [Fe/H]
∼ 0.3, then the second component with a mild decreasing trend are generated at −0.5 <∼ [α/Fe]

– 17 –
<∼ 0 from the SN Ia contribution. Very similar trends are seen for O, Mg, Si, S, and Ca. This
is roughly consistent with the observations, although the double peak distributions are not seen
in the observations. For O, Fulbright, McWilliam, & Rich (2007) showed a steeper decrease than
Mg at [Fe/H] >∼ 0. Theoretically, it is not easy to produce such a different trend of O from Mg
since these two elements are synthesized in the similar region and ejected to the ISM in a similar
way. The contribution from AGB stars is negligible. Such a decrease in O/Mg might require some
additional effects that are not included in our stellar evolution models such as strong stellar winds
or a process that causes the change in the C/O ratio. For other α elements, explosion energy
and remnant mass may cause a variation in (O, Mg)/(Si, S) and (O, Mg)/Ca, respectively. These
should be discussed with future observations with a larger sample and with nucleosynthesis yields
that cover the parameter spaces. The separation between the two populations also depends on the
IMF, and the Salpeter IMF is adopted in this simulation. The IMF will be constrained from the
elemental abundance ratios.
The odd-Z elements and also Zn show the increasing component from [Fe/H] ∼ −3 to ∼ 0,
and the second flat component without any trend at [(Na, Al, Cu, Zn)/Fe] ∼ 0 at 0 <∼ [Fe/H]
<∼ 1. Different from Fig. 16, the metallicity dependence of Na, Al, and Cu yields are not seen very
much in the bulge since the chemical enrichment timescale is short. [Mn/Fe] shows the increasing
component from [Fe/H] ∼ −2 to ∼ 0, which continues to the second flat component at [Mn/Fe]
∼ 0 at 0 <∼ [Fe/H] <∼ 1. This is also due to the smaller contribution from SNe Ia in the bulge. The
observations also show the high [Na/Fe] and [Al/Fe] at [Fe/H] >∼ − 1, and the increasing trend of
[Mn/Fe], although the double peak distributions are not visible.
Figure 18 is for the thick disk, where the observational data are taken from Prochaska et al.
(2000), Bensby et al. (2004), Reddy, Lambert, & Prieto (2006), and Reddy & Lambert (2008).
When we select the thick disk stars from kinematics in the solar neighborhood, the [α/Fe] plateau
is dominant, which is similar to the bulge. This is also consistent with the RAVE observations
(Ruchti et al. 2010). We should note, however, that the split trends observed in the solar neigh-
borhood for the [α/Fe] trends associated with the thick and thin disks are not clearly seen in the
predicted distributions. The formation timescale of the thick disk is 3 − 4 Gyr in our simulation,
which is shorter than the thin disk but longer than the bulge. In Figures 7, 9, and 15, the functions
of the thick disk are also shown (dot-dashed lines). For thick disks stars, the age is older and [O/Fe]
is higher than thin disk stars, and the peak metallicity is 0.2 dex lower than the thin disk.
Figure 19 shows the distribution functions of [X/Fe] in the solar neighborhood (solid line),
bulge (dashed line), and thick disk (dotted line). α elements show the double peak distributions,
and the [α/Fe]-rich peak is at the same [α/Fe] for the solar neighborhood, bulge, and thick disk. On
the other hand, Na, Al, Cu, and Zn show the single peak distributions in the solar neighborhood,
but the double peak distributions in the bulge. On the average, Na, Al, Cu, and Zn are much more
enhanced in the bulge, and (Na, Al, Cu, Zn)/α is larger than in the solar neighborhood. In the thick
disk, the [(Na, Al, Cu, Zn)/Fe] distributions are rather similar to the solar neighborhood because
the metallicity is not as high as in the bulge. The stellar population of the thick disk is neither

– 18 –
thin-disk like nor bulge like (see Tab. 3). For the thick disk stars, the star formation timescale is
as short as in the bulge, but the chemical enrichment efficiency is not as high as in the bulge. This
is because half of the thick disk stars have already formed in satellite galaxies before they accrete
onto the disk, and the metals have been ejected from the satellite galaxies by the galactic winds.
3.7. [α/Fe]-[Mn/Fe] Diagram
In Figure 20, for the solar neighborhood of the present-day galaxy, [Mn/Fe] is plotted against
[α/Fe]=([O/Fe]+[Mg/Fe])/2, which clearly shows the sequence of the SN Ia contribution. With
SNe II and HNe only, [α/Fe] is as high as ∼ 0.5, and [Mn/Fe] is as low as ∼ −0.5. With more
SNe Ia, [α/Fe] decreases, while [Mn/Fe] increases. The three populations of the observed stars are
along this trend; i) the EMP stars (O in C04, large open circles; Mg in H04, filled pentagons) are
found in the left-bottom region with high [α/Fe] and low [Mn/Fe]. ii) The thick disk stars (small
open circles) populate the following region, [α/Fe] ∼ 0.2 − 0.4 and [Mn/Fe] ∼ −0.4 to −0.2. iii)
The thin disks stars (small closed circles) occur at [α/Fe] ∼ 0.1 and [Mn/Fe] ∼ −0.1, which are
formed from the ISM largely enriched by SNe Ia. In other words, it is possible to select thick disk
stars only from the elemental abundance ratios.
The scatter should come from the variety in the SN II yields. There is a larger scatter in
H04 than C04, which may be partially caused by the observational error; the right-bottom star
CS22952-015 is (-0.1, -0.34) in H06, but (0.36, -0.33) in C04 and A10. On the other side, the
left-top star BS16085-050 is at (0.7, -0.01), which has to be checked. If it is real, the enrichment
source may be faint SNe.
4. Conclusions
We present the chemodynamical simulations of a Milky Way-type galaxy using a self-consistent
hydrodynamical code with supernova feedback and chemical enrichment. In our nucleosynthesis
yields of core-collapse supernovae, the light curve and spectra fitting of individual supernova are
used to estimate the mass of the progenitor, explosion energy, and ejected iron mass. A large con-
Table 3. Stellar Populations in Each Component.
Component age [(O, Mg, Si, S, Ca)/Fe] [(Na, Al, Cu)/Fe] [Zn/Fe] [Mn/Fe]
solar neighborhood young low ∼ 0 ∼ 0 ∼ 0
bulge old high high ∼ 0.2 low
thick disk old high ∼ 0 ∼ 0 low

– 19 –
tribution from hypernovae is required from the observed abundance of Zn ([Zn/Fe] ∼ 0) especially
at [Fe/H] <∼ − 1. In our progenitor model of SNe Ia, based on the single degenerate scenario, the
SN Ia lifetime distribution spans a range of 0.1 − 20 Gyr with the double peaks at ∼ 0.1 and 1
Gyr. Because of the metallicity effect of white dwarf winds, the SN Ia rate is very small at [Fe/H]
<∼ − 1, which plays an important role in chemical evolution of galaxies.
In the simulated galaxy, the kinematical and chemical properties of the bulge, disk, and halo
are consistent with the observations. The bulge formed from the assembly of subgalaxies at z >∼ 3;
80% of bulge stars are older than ∼ 10 Gyr, and 60% have [O/Fe] > 0.3. The disk formed with
constant star formation over 13 Gyr; 50% of solar-neighborhood stars are younger than ∼ 8 Gyr,
80% have [O/Fe] < 0.3. When we define the thick disk from kinematics, the thick disk stars tend
to be older and have higher [α/Fe] than the thin disk stars. The formation timescale of the thick
disk is 3− 4 Gyr.
Because the star formation history is different for different components, the age-metallicity
relation and the metallicity distribution function are also different. The age-metallicity relation
shows a more rapid increase in the bulge than in the disk. In both cases, the average metallicity
does not show strong evolution at t >∼ 2 Gyr, as in the observations. The scatter is originated from
the inhomogeneity of chemical enrichment in our chemodynamical model. The observed metallicity
distribution function is better reproduced with UV background radiation and hypernovae, but is
still problematic. The bulge wind induced by SNe Ia seems to be a good solution to reduce the
numbers of metal-rich stars in the bulge and of metal-poor stars in the disk.
The difference in the chemical enrichment timescales results in the difference in the elemen-
tal abundance ratios, since different elements are produced by different supernovae with different
timescales. We also predict the frequency distribution of elemental abundance ratios as functions
of time and location, which will be statistically compared with a large homogeneous sample from
galactic archeology surveys such as HERMES, when they become available.
• Because of the delayed enrichment of SNe Ia, α elements (O, Mg, Si, S, and Ca) show a
plateau at [Fe/H] ∼ −1, and then the decreasing trend against [Fe/H], where [Mn/Fe] also
shows the increasing trend. Odd-Z elements (Na, Al, and Cu) show the increasing trend at
[Fe/H] <∼ − 1 because of the metallicity dependence of nucleosynthesis yields. These are in
excellent agreement with the available observations.
• In the bulge, the star formation timescale is so short that the [α/Fe] plateau continues to
[Fe/H] ∼ +0.3. Because of the smaller contribution from SNe Ia, the majority of stars
shows high [α/Fe] and low [Mn/Fe]. [(Na, Al, Cu, Zn)/Fe] are also high because of the high
metallicity in the bulge.
• The stellar population of the thick disk is neither disk-like nor bulge-like as summarized in
Table 3. For thick disk stars, [α/Fe] is higher, and [Mn/Fe] is lower than thin disk stars
because of the short formation timescale. However, [(Na, Al, Cu, Zn)/Fe] are lower than

– 20 –
bulge stars because of the lower chemical enrichment efficiency. This is because half of the
thick disk stars have already formed in satellite galaxies before they accrete onto the disk,
and the metals have been ejected from the satellite galaxies by the galactic winds.
We thank the National Astronomical Observatory of Japan for the GRAPE system, where
most of the simulations in this paper are performed. We also thank K. Nomoto, M. Mori, K. C.
Freeman, G. Da Costa, and D. Yong for fruitful discussions.
REFERENCES
Anders, E., & Grevesse, N. 1989, Geochimica et Cosmochimica Acta, 53, 197
Andrievsky, S. M., Spite, M., Korotin, S. A., Spite, F., Bonifacio, P., Cayrel, R., Franc¸ois, P., &
Hill, V. 2010, A&A, 509, 88
Andrievsky, S. M., Spite, M., Korotin, S. A., Spite, F., Bonifacio, P., Cayrel, R., Hill, V., &
Franc¸ois, P. 2007, A&A, 464, 1081
Andrievsky, S. M., Spite, M., Korotin, S. A., Spite, F., Bonifacio, P., Cayrel, R., Hill, V., &
Franc¸ois, P. 2008, A&A, 481, 481
Argast, D., Samland, M., Thielemann, F.-K., & Gerhard, O. E., 2002, A&A, 388, 842
Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
Beers, T. 2010, in IAU Symposium 265 Chemical Abundances in the Universe: Connecting First
Stars to Planets, eds. K. Cunha et al., p.453
Bensby, T., Feltzing, S., & Lundstro¨m, I. 2003, A&A, 410, 527
Bensby, T., Feltzing, S., & Lundstro¨m, I. 2004, A&A, 415, 155
Bergemann, M., & Gehren, T. 2008, A&A, 492, 823
Bergemann, M., & Gehren, T. 2010, in IAU Symposium 265 Chemical Abundances in the Universe:
Connecting First Stars to Planets, eds. K. Cunha et al., p.348
Bergemann, M. Pickering, J. C., & Gehren, T. 2010, MNRAS, 401, 1334
Bertschinger, E. 1995, http://arcturus.mit.edu/cosmics/
Bromm, V., & Larson, R. B. 2004, ARA&A, 42, 79
Bruenn, S. W., Mezzacappa, A., Hix, W. R., Blondin, J. M., Marronetti, P., Messer, O. E. B.,
Dirk, C. J., & Yoshida, S. 2009, AIP Conference Proceedings, Vol. 1111, p. 593
Cayrel, R., et al. 2004, A&A, 416, 1117
Chen, Y. Q., Nissen, P. E., Zhao, G., & Asplund, M. 2002, A&A, 390, 225
Cunha, K., Sellgren, K., Smith, V. V.. Ramirez, S. V.. Blum, R. D., & Terndrup, D. M., 2007,
ApJ, 669, 1011
Diemand, J., Kuhlen, M., & Madau, P. 2007, ApJ, 657, 262

– 21 –
Edvardsson, B., Andersen, J., Gustafsson, B., Lambert, D. L., Nissen, P. E., & Tomkin, J. 1993,
A&A, 275, 101
Feltzing, S., Fohlman, M., & Bensby, T. 2007, A&A, 467, 665
Freeman, K. C., et al. 2010, in American Astronomical Society Meeting, Vol. 42, p.478
Fulbright, J. P. 2000, AJ, 120, 1841
Fulbright, J. P., McWilliam, A., & Rich, R. M. 2007, ApJ, 661, 1152
Gratton, R. G., et al. 2003, A&A, 404, 187
Hachisu, I., Kato, M., Nomoto, K., & Umeda, H. 1999, ApJ, 519, 314
Holmberg, J., Nordstro¨m, B., & Andersen, J. 2007, A&A, 475, 519
Honda, S., et al. 2004, ApJ, 607, 474
Honma, M., et al. 2007, PASJ, 59, 889
Hopkins, P. F., et al. 2009, MNRAS, 397, 802
Israelian, G., Garc´ıa-Lo´pez, R. J., & Rebolo, R. 1998, ApJ, 207, 805
Israelian, G., & Reboro, R. 2001, ApJ, 557, L43
Iwamoto, K., Brachwitz, F., Nomoto, K., Kishimoto, N., Umeda, H., Hix, W. R., & Thielemann,
F.-K. 1999, ApJS, 125, 439
Katz, N. 1992, ApJ, 391, 502
Kauffmann, G., & White, S. D. M. 1993, MNRAS, 261, 921
Kitaura, F. S., Janka, H.-Th., & Hillebrandt, W. 2006, A&A, 450, 345
Klypin, A., Zhao, H., & Somerville, R. S. 2002, ApJ, 573, 597
Kobayashi, C., 2004, MNRAS, 347, 740
Kobayashi, C., 2005, MNRAS, 361, 1216
Kobayashi, C. 2010, in IAU Symposium 265 Chemical Abundances in the Universe: Connecting
First Stars to Planets, eds. K. Cunha et al., p.336
Kobayashi, C., Nakasato, N., & Nomoto, K. 2003, in IAU Symposium 208 Astrophysical Super-
computing using Particle Simulations, 2001, eds. J. Makino & P. Hut (Dordrecht: Kluwer),
p.419
Kobayashi, C., & Nomoto, K. 2009, ApJ, 707, 1466 (KN09)
Kobayashi, C., Karakas, I. A., & Umeda, H. 2010a, MNRAS, submitted
Kobayashi, C., Tominaga, N., & Nomoto, K. 2010b, ApJ, submitted, arXiv:1101.1227
Kobayashi, C., Springel, V., & White, S. D. M. 2007, MNRAS, 376, 1465 (K07)
Kobayashi, C., Tsujimoto, T., Nomoto, K., Hachisu, I, & Kato, M. 1998, ApJ, 503, L155
Kobayashi, C., Umeda, H., Nomoto, K., Tominaga, N., & Ohkubo, T. 2006, ApJ, 653, 1145 (K06)
Koda, J., Milosavljevic´, M., & Shapiro, P. R. 2009, ApJ, 696, 254
Kodama, T., & Arimoto, N. 1997, A&A, 320, 41
Lecureur, A., et al. 2007, A&A, 465, 799

– 22 –
Levine, E. S., Heiles, C., & Blitz, L. 2008, ApJ, 679, 1288
Lo´pez-Corredoira, M., Cabrera-Lavers, A., Garzo´n, F., & Hammersley, P. L. 2002, A&A, 394, 883
Macial, W. J., Lago, L. G., & Costa, R. D. D. 2006, A&A, 453, 587
Maeda, K., & Nomoto, K. 2003, ApJ, 598, 1163
Mannucci, F., Della Valle, M., & Panagia, N. 2006, MNRAS, 370, 773
Marek, A., & Janka, H.-Th. 2009, ApJ, 694, 664
McWilliam, A., & Rich, R. M. 1994, ApJS, 91, 749
Merrifield, M. R. 1992, AJ, 103, 1552
Meynet, G., & Maeder, A. 2002, A&A, 390, 561
Moore, B., Ghigna, S., Governato, F., Lake, G., Quinn, T., Stadel, J., & Tozzi, P. 1999, ApJ, 524,
19
Morrison, H. 2010, talk at the conference on The Chemical Enrichment of the Milky Way Galaxy,
at Ringberg Castle, Germany
Nakasato, N., & Nomoto, K. 2003, ApJ, 588, 842
Navarro, J. F., & White, S. D. M. 1993, MNRAS, 265, 271
Nissen, P. E., Akerman, C., Asplund, M., Fabbian, D., Kerber, F., Ka¨ufl, H. U., & Pettini, M.
2007, A&A, 469, 319
Nomoto, K. 1982, ApJ, 253, 798
Nomoto, K., & Kondo, Y. 1991, ApJ, 367, L19
Nomoto, K., Tominaga, N., Umeda, H., Kobayashi, C, & Maeda, K. 2006, Nuclear Physics A, 777,
424
Nomoto, K., et al. 1997a, Nuclear Physics, A616, 79c
Nomoto, K., et al. 1997b, Nuclear Physics, A621, 467c
Nordstro¨m, B., et al. 2004, A&A, 418, 989
Ohkubo, T., Umeda, H., Maeda, K., Nomoto, K., Suzuki, T., Tsuruta, S., & Rees, M. J. 2006,
ApJ, 645, 1352
Podsiadlowski, P., Mazzali, P. A., Nomoto, K., Lazzati, D., & Cappellaro, E. 2004, ApJ, 607, L17
Primas, F., Reimers, D., Wisotzki, L., Reetz, J., Gehren, T., & Beers, T. C. 2000, in The First
Stars, ed. A. Weiss, T. Abel, &, V. Hill (Berlin:Springer), 51
Prochaska, J., Naumov, S. O., Carney, B, W,, McWilliam, A., & Wolfe, A. M. 2000, AJ, 120, 2513
Reddy, B. E., & Lambert, D. L. 2008, MNRAS, 391, 95
Reddy, B. E., Lambert, D. L., & Prieto, C. A. 2006, MNRAS, 367, 1329
Reddy, B. E., Tomkin, J., Lambert, D. L., & Prieto, C. A. 2003, MNRAS, 340, 304
Renzini, A., Ciotti, L., D’Ercole, A., & Pellegrini, S. 1993, ApJ, 419, 52
Ruchti, G. R., et al. 2010, ApJ, 721, L92
Sadler, E, M., Rich, R. M., & Terndrup, D. M. 1996, AJ, 112, 171

– 23 –
Saito, Y., Takada-Hidai, M., Honda, S., & Takeda, Y. 2009, PASJ, 61, 549
Sa´nchez-Bla´zquez, P., Courty, S., Gibson, B. K., & Brook, C. B. 2009, MNRAS, 398, 591
Scannapieco, C., White, S. D. M., Springel, V., & Tissera, P. B. 2009, MNRAS, 396, 696
Sommer-Larsen, J., Go¨lz, M., & Portinari, L. 2003, ApJ, 596, 47
Steinmetz, M., & Navarro, J. F. 1999, ApJ, 513, 555
Sullivan, M., et al. 2006, ApJ, 648, 868
Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 235
Takada-Hidai, M., et al. 2005, PASJ, 57, 347
Thielemann, F.-K., Nomoto, K., & Hashimoto, M. 1996, ApJ, 460, 408
Tinsley, B. M. 1980, Fundamentals of Cosmic Phisics Vol.5, p.287
Toth, G., & Ostriker, J. P. 1992, ApJ, 389, 5
Tsujimoto, T., Nomoto, K., Yoshii, Y., Hashimoto, M., Yanagida, S., & Thielemann, F.-K. 1995,
MNRAS, 277, 945
Warren, M. S., Quinn, P. J., Salmon, J. K., & Zurek, W. H. 1992, ApJ, 399, 405
Wyse, R. F. G., & Gilmore, G. 1995, AJ, 110, 2771
Zoccali, M., et al. 2008, A&A, 486, 177
This preprint was prepared with the AAS LATEX macros v5.2.

– 24 –
Fig. 1.— Projected stellar density maps for the edge-on views at z = 1.28 for 150 galaxies from
various initial conditions.

– 25 –
Fig. 2.— Time evolution of the projected density of dark matter (first row), gas (second row),
V-band luminosity (third row), and the luminosity-weighted stellar metallicity (forth row) in the
range of [M/H]= −1 to 0.1 at z = 8, 5, 3, 2.3, 1.2, and 0. The leftmost panels are 200 kpc and the
other panels are 20 kpc on a side. The rightmost panels are for the face-on views at z = 0.

– 26 –
Fig. 3.— B-band (upper panels) and K-band (lower panels) luminosities for the edge-on (left
panels) and face-on (right panels) views at z = 0 in 20 kpc on a side.

– 27 –
Fig. 4.— Surface density profile (top panel) and density profile (bottom panel) at r = 7 − 9 kpc
for total (dotted lines), thick disk (dot-dashed lines), and thin disk (solid lines) at present.

– 28 –
Fig. 5.— Radial metallicity gradients of stars at the height |z| = 0−0.5 (solid line), 0.5−1 (dashed
line), 1− 1.5 (dot-dashed line), and 5− 10 (triple-dot-dashed line) kpc at present (top panel) and
at 8 Gyr (bottom panel). The dotted lines show the radial gradients of gas at |z| = 0 − 0.5 kpc.
The metallicities are mass-weighted.

– 29 –
Fig. 6.— Star formation histories in the solar neighborhood (solid line), bulge (dashed line), and
inner (dotted line) and outer (dot-dashed line) halo at z = 0.

– 30 –
Fig. 7.— Cumulative functions of the ages of stars in the solar neighborhood (solid line), bulge
(dashed line), halo (dotted line), and thick disk (dot-dashed line) at z = 0.

– 31 –
Fig. 8.— Age-metallicity relations in the solar neighborhood (upper panel) and bulge (lower panel)
at z = 0. The contours show the frequency distribution of stars in the simulated galaxy, and
red is for the highest frequency. The white dots show the observations in the solar neighborhood
(Holmberg et al. 2007).

– 32 –
Fig. 9.— Cumulative functions of [O/Fe] for stars in the solar neighborhood (solid line), bulge
(dashed line), halo (dotted line), and thick disk (dot-dashed line) at z = 0.

– 33 –
Fig. 10.— [O/Fe]-[Fe/H] relations in the solar neighborhood (upper panel) and bulge (lower panel)
at z = 0. The contours show the frequency distribution of stars in the simulated galaxies, and red is
for the highest frequency. The white dots show the observations of stars: for the solar neighborhood,
Edvardsson et al. (1993), small open circles; thin and thick disk stars in Bensby et al. (2004), filled
and open triangles, respectively; dissipative component in Gratton et al. (2003), filled squares; and
Cayrel et al. (2004), large open circles. For the bulge, Fulbright, McWilliam, & Rich (2007), filled
stars; and Lecureur et al. (2007), open circles.

– 34 –
Fig. 11.— Distribution functions of [O/Fe] for stars in the solar neighborhood (z = 0) at [Fe/H]
= −3.5 to −2.5 (solid line), −2.5 to −1.5 (dashed line), −1.5 to −0.5 (dotted line), −0.5 to 0.5
(dot-dashed line), and 0.5 to 1.5 (triple-dot-dashed line).

– 35 –
Fig. 12.— Distribution functions of [O/Fe] for stars in the solar neighborhood (z = 0) at −0.5 ≤
[Fe/H] ≤ −0.3 for the standard model (solid line) and the models with the star formation timescale
c = 0.02 (dashed line), the smaller feedback radius rFB = 0.5 kpc (dotted line), and the fixed
feedback number NFB = 50 (dot-dashed line). The arrow shows the peak value of observations.

– 36 –
Fig. 13.— Distribution functions of [O/Fe] for stars in the solar neighborhood (z = 0) at ≤ −4.0
[Fe/H] ≤ −1.9 for the models with (solid line) and without (dashed line) hypernova feedback. The
arrow shows the peak value of observations.

– 37 –
Fig. 14.— Metallicity distribution functions of stars in the solar neighborhood at z = 0 for the
models without UV background radiation and hypernovae (dotted line), with UV background
radiation but without hypernovae (dashed line), and with UV background radiation and hypernovae
(solid line). The dots with errorbars are for the observations taken from Edvardsson et al. (1993,
open circles) and Holmberg et al. (2007, filled circles).

– 38 –
Fig. 15.— Metallicity distribution functions of stars in the solar neighborhood (solid line), bulge
(dashed line), and thick disk (dot-dashed line) at z = 0. The observational data sources are:
Edvardsson et al. (1993), open circles; Holmberg et al. (2007), filled circles; McWilliam & Rich
(1994), open squares; and Zoccali et al. (2008), filled squares.

– 39 –
Fig. 16.— [X/Fe]-[Fe/H] relations in the solar neighborhood at z = 0. The contours show the
frequency distribution of stars in the simulated galaxies, where red is for the highest frequency. The
observational data (white dots) are taken from Cayrel et al. (2004), large open circles; Honda et al.
(2004), filled pentagons; Fulbright (2000), crosses; Reddy et al. (2003, 2006) and Reddy & Lambert
(2008) for thin (small filled circles) and thick (small open circles) disk stars. For O, Bensby et al.
(2004) for thin (filled triangles) and thick (open triangles) disk stars. For Si, S, and Zn, Chen et al.
(2002), filled squares; Takada-Hidai et al. (2005), filled pentagon; Nissen et al. (2007), open squares.
For Cu and Zn, Primas et al. (2000), asterisks.

– 40 –
Fig. 17.— Same as Fig. 16 but for the bulge. The observational data (white dots) are taken from
Fulbright, McWilliam, & Rich (2007), filled stars; and Lecureur et al. (2007), open circles.

– 41 –
Fig. 18.— Same as Fig. 16 but for the thick disk. The observational data (white dots) are
taken from Prochaska et al. (2000), six-pointed stars; Bensby et al. (2004), open triangles; and
Reddy, Lambert, & Prieto (2006) and Reddy & Lambert (2008), open circles.

– 42 –
Fig. 19.— Distribution functions of the elemental abundance ratios relative to iron [X/Fe] in the
solar neighborhood (solid line), bulge (dashed line), and thick disk (dot-dashed line) at z = 0.

– 43 –
Fig. 20.— [α/Fe]-[Mn/Fe] relation in the solar neighborhood at z = 0. The contours show the
frequency distribution of stars in the simulated galaxies, where red is for the highest frequency. The
observational data (white dots) are taken from Cayrel et al. (2004), large open circles; Honda et al.
(2004), filled pentagons; Reddy et al. (2003, 2006) and Reddy & Lambert (2008) for thin (small
filled circles) and thick (small open circles) disk stars.