CHEMODYNAMICAL SIMULATIONS OF THE MILKY WAY GALAXY

The details of our GRAPE-SPH code are described in Kobayashi (2004) and can be summarized as follows.

• 1.  
  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 hydrodynamics, the Smoothed Particle Hydrodynamical (SPH) method is adopted, and the SPH formulation is the almost the same as Navarro & White's (1993). The GRAPE-SPH code is highly adaptive in space and time by means of individual smoothing lengths and individual time steps, and it has very high performance (a week for one simulation with N ∼ 50, 000).
• 2.  
  Radiative cooling is computed with the metal-dependent cooling functions, which are generated 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.
• 3.  
  Our star formation criteria are the same as Katz's (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 ($t_{\rm sf} \equiv \frac{1}{c}t_{\rm dyn}$), 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 Figure 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.
• 4.  
  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 10⁵–10⁶ M_☉, 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–120 M_☉.
• 5.  
  For the feedback of energy and heavy elements, we do not adopt the instantaneous recycling approximation. Via stellar winds, SNe II/HNe, and SNe Ia, thermal energy and heavy 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 r_FB. 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 (N_FB) in galaxy-scale simulations; the metallicity of gas particles does not depend very much on the local density if N_FB is fixed. We find that r_FB = 1 kpc gives a better fit to the elemental abundance ratios than r_FB = 0.5 and 3 kpc, or N_FB = 50 and 100 (see Figure 12). Note that the diffusion between gas particles is supposed to enhance the mixing of heavy elements but is not included in our model.
• 6.  
  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 function 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.

The scheme to include the detailed chemical enrichment into SPH simulations has been proposed by Kobayashi (2004) and updated as follows.

Hypernovae (HNe). The explosion mechanism of core-collapse supernovae is still uncertain, although a few groups have presented feasible calculations in exploding core-collapse supernovae (Kitaura et al. 2006; Marek & Janka 2009; Bruenn et al. 2009). However, the ejected explosion energy and ⁵⁶Ni mass (which decays to ⁵⁶Fe) 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 ⩾ 20 M_☉) have more than 10 times larger explosion energy (E₅₁ ≳ 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 progenitor's mass range is 8–50 M_☉ but the metal ejection is only from 10–50 M_☉. Assuming that a large fraction of supernovae with M ⩾ 20 M_☉ 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 two-dimensional calculation of nucleosynthesis (Maeda & Nomoto 2003). The observed increasing trend toward the lower metallicity may suggest 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 ⩾ 20 M_☉), _HN = 0.5, 0.5, 0.4, 0.01, and 0.01 for Z = 0.0.001, 0.004, 0.02, and 0.05,⁵ 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 et al. 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 ∼1 M_☉ and near main-sequence (MS) companions with initial masses of ∼3 M_☉ (Hachisu et al. 1999).

In our simulations, based on the SD scenario, the lifetime distribution function of SNe Ia is calculated with Equation (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 MDF 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 metallicity 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 WD (Iwamoto et al. 1999).

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–8 M_☉ (depending on metallicity) produce light elements such as C and N, while the contribution of heavier elements are negligible in the galactic chemical evolution (Kobayashi 2010; Kobayashi et al. 2010a). Rotating massive stars (≳40 M_☉) could also produce C and N (e.g., Meynet & Maeder 2002) but are not included.

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 H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, Ω_Λ = 0.7, Ω_b = 0.04, n = 1, and σ₈ = 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 ∼10¹² M_☉ 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⁻¹M^−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 does not seem 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 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 until z = 1.28 including gas dynamics and basic chemical enrichment (Figure 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 Figure 1) for detailed simulations.

Zoom In Zoom Out Reset image size

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 Figures 2–11 and 15, and with the standard initial condition in Figures 12, 13, 14, and 16–20.

Zoom In Zoom Out Reset image size

Figure 2 shows the time evolution of the projected density of dark matter, gas, V-band luminosity, 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. According 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) × 10⁹ M_☉ 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 ongoing 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 and the alignment 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.

Zoom In Zoom Out Reset image size

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 et al. 2002). Note that it is not feasible to resolve the central black hole, nuclear bulge, and bar in our simulations. The total masses inside 1, 10, and 100 kpc are 9.4 × 10⁹ M_☉, 2.2 × 10¹¹ M_☉, and 1.1 × 10¹² M_☉, 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 kpc ≲ r ≲ 20–60 kpc (Merrifield 1992; Honma et al. 2007; Levine et al. 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.5kpc ⩽ r ⩽ 8.5 kpc and the height of |z| ⩽ 0.5 kpc for the solar neighborhood, r ⩽ 1 kpc for the bulge, and 5kpc ⩽ 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_Θ, V_Z). The number of particles significantly drops around the velocity dispersion $\sigma \equiv \sqrt{\vphantom{A^A}\smash{\hbox{${V_\Phi ^2+V_\Theta ^2+V_Z^2}$}}} \sim 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; Ló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.

Zoom In Zoom Out Reset image size

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 et al. 2006).

Zoom In Zoom Out Reset image size

The resultant star formation histories (Figure 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 the 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 the stars are as old as the 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 Sánchez-Blázquez et al. (2009).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 (Section 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 (Section 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: (1) the age–metallicity relation is weak. (2) SNe Ia can affect the elemental abundance ratios at [Fe/H] ≲−1 even with the metallicity inhibition of SNe Ia. (3) 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 Section 2.1, we run simulations with several feedback parameters and schemes, and show the best results in this paper.

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 a given time is caused by the inhomogeneity of chemical enrichment in our chemodynamical model (Section 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).

Zoom In Zoom Out Reset image size

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.

The difference in the chemical enrichment timescales results in a difference in the elemental 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 Section 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.

Zoom In Zoom Out Reset image size

Figure 10 shows the [O/Fe]–[Fe/H] relations, and the other α-elements show the same trends. 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 (Section 2.2).

Zoom In Zoom Out Reset image size

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 (Figure 10). This may be because the mixing of heavy elements among gas particles is not included in our model (Section 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 Section 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 Section 3.6 for more discussion). Second, there is an intrinsic variation in [α/Fe] of supernova yields depending on the progenitor mass and metallicity, the effect of which is included in our chemodynamical model. Low-mass SNe II (10–13 M_☉) 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 black hole) for SNe II and HNe.

Zoom In Zoom Out Reset image size

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 (r_FB = 0.5 kpc; dotted line), and a fixed feedback number (N_FB = 50; dot-dashed line). The smaller initial condition (Table 1) is used here. The standard model gives a much better agreement with the peak value of the observations (arrow) than the other models.

Zoom In Zoom Out Reset image size

Using the 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 process 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.15 M_☉ for 13–15 M_☉). Now the iron mass is estimated from the light curve and spectral fitting, and is 0.07 M_☉ for SNe II with <20 M_☉. The high [α/Fe] values are reduced by including HNe. Figure 13 shows the distribution functions at −4⩽ [Fe/H] ⩽−1.9 for the models with (solid line) and without (dashed line) HN feedback using the smaller initial condition. With HNe, the scatter around [O/Fe] ∼0.5 is as small as the observations (Cayrel et al. 2004).

Zoom In Zoom Out Reset image size

In the bulge (lower panel of Figure 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 et al. 2007; Lecureur et al. 2007), except for the O observation in Fulbright et al. (2007; see Section 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 (1) local enrichment by the stellar winds of dying stars or SNe II, (2) intrinsic variation of supernovae depending on energy and mixing, or (3) the selective mass-loss of the gas enriched by SNe Ia. The bulge wind might be driven by SNe Ia, and most of the 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 (1) and (2) may be preferable.

The 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 HN (solid line) gives much better agreement with the observations in the solar neighborhood. 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 HN 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. 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.

Zoom In Zoom Out Reset image size

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 Nordströ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. The G-dwarf problem also seems to exist 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.

Zoom In Zoom Out Reset image size

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–270 M_☉ stars. Stars with ≲100 M_☉ and ≳300 M_☉ are usually supposed to collapse to black holes, 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.

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 16 is for the solar neighborhood. The dots are for the observations of individual stars, 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 et al. (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 three-dimensional corrections of −0.24 for O in C04, and the non-local thermodynamic equilibrium (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; Andrievsky et al. 2010, hereafter A10) with the NLTE corrections. For O, S, Cu, and Zn, we employ other data sources as noted below.

• 1.  
  α 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 three-dimensional 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.'s (2004) data with the forbidden line [O i] 6300 Å since an empirical relation using [O i] and O i triplet at 7775 Å 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.
• 2.  
  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 Section 3.2, in our chemodynamical simulations, there is only a weak 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.
• 3.  
  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 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 a higher energy and/or a higher fraction of HNe where more Co and Zn are synthesized.
• 4.  
  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 a large contribution of HNe. In detail, there is a small oscillating trend; [Zn/Fe] is 0 at [Fe/H] ∼0, increases to 0.2 at [Fe/H] ∼−0.5, decreases to 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 found by Saito et al. (2009). Theoretically, Zn production depends on many parameters; ⁶⁴Zn 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.
• 5.  
  Manganese—Mn is a characteristic element and is produced more 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 for both SNe II and SNe Ia. Feltzing et al. (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 metallicities, the [Mn/Fe] scatter, if it is real, should stem from the SN II yields depending on the remnant mass and the electron fraction Y_e 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 (black holes; Kobayashi et al. 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. (2009). This means the normalization should be changed both for observations and simulations. When the observational data are analyzed with three-dimensional 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 et al. (2007) and Zoccali et al. (2008). As shown in Figure 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] ≲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 et al. (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 similar regions 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 Figure 16, the metallicity dependence of the 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 et al. (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 neighborhood 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 disk 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 thin-disk like nor bulge like (see Table 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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: (1) 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]. (2) The thick disk stars (small open circles) populate the following region, [α/Fe] ∼0.2–0.4 and [Mn/Fe] ∼−0.4 to −0.2. (3) 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.