{'fraction_non_alphanumeric': 0.037242617425075106, 'fraction_numerical': 0.019766248992452554, 'mean_word_length': 4.422030194676202, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
GCD+: A New Chemodynamical Approach to Modeling SNe and Chemical Enrichment

We have developed a new galactic chemo-dynamical evolution code, called GCD+, for studies of galaxy formation and evolution. This code is based on our original three-dimensional tree N-body/smoothed particle hydrodynamics code which includes self-gravity, hydrodynamics, radiative cooling, star formation, supernova feedback, and metal enrichment. GCD+ includes a new Type II (SNe II) and Ia (SNe Ia) supernovae model taking into account the lifetime of progenitor stars, and chemical enrichment from intermediate mass stars. We apply GCD+ to simulations of elliptical galaxy formation, and examine the colour-magnitude relation (CMR), the Kormendy relation, and the [Mg/Fe]--magnitude relation of simulation end-products. GCD+ is a useful and unique tool which enables us to compare simulation results with the observational data directly and quantitatively. Our simulation confirm the results of Kawata (2001) who uses a simpler chemo-dynamical evolution code. We newly find that radiative cooling becomes more efficient and thus the gas infall rate increases, with decreasing mass of galaxies, which contributes to the slope of the CMR. In addition, the sophisticated treatments of both SNe II and SNe Ia in GCD+ show that feedback from SNe Ia plays a crucial role in the evolution of elliptical galaxies. We conclude that the feedback effect of SNe Ia should not be ignored in studying the evolution of elliptical galaxies.


INTRODUCTION
The abundances of elements heavier than helium ("metals") in interstellar gas and in stars vary systematically from place to place within galaxies. Since heavy elements are believed to have been synthesized in stars, the metal abundances in present-day galaxies should offer a record through which we may trace the history of star formation within galaxies. Since the pioneering work of Tinsley (1972), studies of chemical evolution have succeeded in connecting observed metal abundances to the history of star formation within galaxies (e.g. Arimoto & Yoshii 1987;Matteucci & Tornambè 1987;Timmes, Woosley, & Weaver 1995;Gibson 1997;Chappini, Matteucci, & Gratton 1997;Tantalo et al. 1998b). Moreover, recently developed semi-analytic models (White & Rees 1987;White & Frenk 1991) including chemical evolution make it possible to relate the observed chemical properties of galaxies to cosmological theory (e.g. Kauffmann & ⋆ E-mail: dkawata,bgibson@astro.swin.edu.au Charlot 1998; Cole et al. 2000;Okamoto & Nagashima 2001;Somerville, Primack, & Faber 2001;Beasley et al. 2002). These techniques require low computational costs and allow large coverage of parameter space. Although the pure chemical evolution models and the semi-analytic models seem to succeed in explaining the observational properties of galaxies at various redshifts, they inevitably assume a phenomenological model involving a number of parameters to describe the processes of gas accretion rates, radiative cooling, star formation, and supernova feedback within a galactic halo. In addition, structures of the end-products can not be discussed, because there is no information about the dynamical history.
On the other hand, recent advances in computer technology and numerical methods have made it possible to calculate the dynamical and chemical evolution of galaxies selfconsistently (e.g. Carraro, Lia, & Chiosi 1998;Mosconi et al. 2001;Chiosi & Carraro 2002). Three dimensional chemodynamical evolution codes have been applied to the following studies over the past five years: Steinmetz & Müller (1995) studied the chemo-dynamical evolution of disk galaxies and succeeded in distinguishing the chemical properties between halo, bulge, and disk stars (see also Bekki & Chiba 2000. Raiteri, Villata, & Navarro (1996) and Berczik (1999) took into account metal enrichment by Type Ia supernovae (SNe Ia) as well as Type II supernovae (SNe II), and reproduced the correlation between [O/Fe] and [Fe/H] for stars in the solar neighbourhood. Such self-consistent calculation of chemical evolution allows one to analyse the photometric properties of simulation end-products in combination with stellar population synthesis, and to obtain the luminosity of the end-products without the assumption of an arbitrary mass to light ratio. Taking advantage of such technique, Steinmetz & Navarro (1999), Navarro & Steinmetz (2000), and Koda, Sofue, & Wada (2000) discussed the zero-point of the Tully-Fisher relation of disk galaxies. Bekki & Shioya (1998, 2000 studied the dynamical and photometric properties of elliptical galaxies formed by a merger of two disk galaxies. Mori et al. (1997) showed that the galactic wind caused by supernovae (SNe) feedback leads to a positive colour gradient observed in some dwarf elliptical galaxies. Using cosmological simulation, Gnedin (1998), Cen & Ostriker (1999), Tissera et al. (2001), and Theuns et al. (2002) studied metal enrichment of the intergalactic medium (see also Aguirre et al. 2001, as a work coupling semi-analytic model to hydro simulations), and Nagamine et al. (2001) examined the luminosity function and colour distribution function of galaxies at different redshifts. Since radiative cooling depends on the metallicity, dynamical evolution is quite sensitive to the chemical evolution (Käellander & Hultman 1998;Kay et al. 2000). Thus, we must calculate chemical and dynamical evolution self-consistently.
Using a chemo-dynamical evolution code, Kawata (2001b, hereafter K01b) were able to reproduce the colourmagnitude relation (CMR) of observed elliptical galaxies. K01b showed that SNe feedback affects the evolution of lower mass systems more strongly, and induces a stronger galactic wind in lower mass systems. Such low mass systems have lower metallicities and bluer colours than higher mass systems, as predicted by analytic models (e.g. Arimoto & Yoshii 1987;Gibson 1997) and more phenomenological numerical simulations (e.g. Larson 1974;Carlberg 1984). Moreover, K01b found that galactic winds are triggered mainly by SNe Ia rather than SNe II, although K01b ignored the lifetime of progenitors of SNe II and SNe Ia (i.e. K01b assumes the instantaneous recycling approximation for SNe II, and all SNe Ia occur simultaneously 1.5 Gyr after stars were born).
We have developed a new chemo-dynamical evolution code, called GCD+, which adopts a more sophisticated chemical evolution model than that discussed in K01b. We have relaxed the instantaneous recycling approximation for SNe II, and adopt a more sophisticated SNe Ia model (Kobayashi, Tsujimoto, & Nomoto 2000, hereafter KTN00). Our code takes into account the metallicity dependent lifetime of stars (Kodama 1997), and metal enrichment from intermediate mass stars (van den Hoek & Groenewegen 1997). The purpose of this paper is to introduce details of GCD+. In addition, as an application of GCD+, we carry out simulations of elliptical galaxy formation similar to the ones of K01b. Then we re-examine the CMR, the Kormendy relation, and the [Mg/Fe]-magnitude relation, which were studied in K01b. Since these observed relations are well-established, and are expected to reflect the chemo-dynamical evolution of elliptical galaxies, this is an important application for GCD+. Next, taking advantage of the sophisticated SNe II and SNe Ia models in GCD+, we clarify the role of SNe Ia, compared to SNe II, in the evolution of elliptical galaxies.
The outline of this paper is as follows. Details of GCD+ are described in Section 2. In Section 3, we use GCD+ to simulate elliptical galaxy formation. We show the galaxy formation model used in this paper briefly in Section 3.1. In Sections 3.2-3.6, results of numerical simulations are presented, and are compared to the observed scaling relations. Discussion about the study of elliptical galaxy formation is given in Section 3.7. Finally, we present our conclusions from this paper in Section 4.

THE NEW CHEMO-DYNAMICAL EVOLUTION CODE
The new code, GCD+ is an update version of the code used in Kawata (1999) and K01b. The code is essentially based on TreeSPH (Hernquist & Katz 1989;Katz, Weinberg, & Hernquist 1996), which combines the tree algorithm (Barnes & Hutt 1986) for the computation of the gravitational forces with the smoothed particle hydrodynamics (SPH : Lucy 1977;Gingold & Mnoraghan 1977) approach to numerical hydrodynamics. The dynamics of the dark matter and stars is calculated by the N-body scheme, and the gas component is modeled using the SPH. It is fully Lagrangian, three-dimensional, and highly adaptive in space and time owing to individual smoothing lengths and individual time steps. Moreover, it includes self-consistently almost all the important physical processes in galaxy formation, such as self-gravity, hydrodynamics, radiative cooling, star formation, SNe feedback, and metal enrichment. The code is vectorized and parallelized. The parallelization is done using the MPI library, so that the code can run on any type of computer, including supercomputers and PC clusters. We have revised the modeling of cooling, star formation, and SNe feedback drastically from the code used in K01b. Here, we will only deal with the new aspects of our N-body/SPH code, which has already been described in Kawata (1999).

Gas Cooling
To model gas cooling, we use the cooling curves computed by MAPPINGS III 1 written by R.S. Sutherland. Using MAP-PINGS III, we have computed the cooling function of the ionization equilibrium gas with various metallicities. Fig. 1 shows all the cooling curves which we use. With linear interpolation of these cooling curves, the code calculates the cooling rate of Λ(T, [Fe/H]) as a function of the temperature and metallicity for each gas particle. We assume that the cooling rate of the gas with metallicity of lower than [Fe/H] = −3 is the same as the one with [Fe/H] = −3, and the cooling rate of the gas with [Fe/H] > 0.5 is equivalent to that at [Fe/H] = 0.5 to avoid extrapolation. The cooling rate for the gas with the solar metallicity is larger than that for the gas of primordial composition by more than an order of magnitude, thus cooling by metals should not be ignored in numerical simulations of galaxy formation (Käellander & Hultman 1998;Kay et al. 2000). The temperature for each particle is derived by Ti = Piµmp/(kBρi), where Pi and ρi, are the pressure and density for i-th particle, and µ, kB, and mp are the mean molecular weight, Boltzmann's constant, and the proton mass, respectively. For simplicity, we fix µ = 0.6, regardless of the metallicity. We set the lower limit of the temperature to be T lim = 100 K.

Star Formation and Initial Mass Function
We model star formation using a method similar to that suggested by Katz (1992) and Katz, Weinberg, & Hernquist (1996). The criteria for star formation and the star formation rate are the same as those in K01b. We use the following three criteria for star formation: 1) the gas density is greater than a critical density, ρcrit = 2 × 10 −25 g cm −3 , i.e. nH ∼ 0.1cm −3 ; 2) the gas velocity field is convergent, ∇ · vi < 0; and 3) the Jeans unstable condition, h/cs > tg, is satisfied, here h, cs, and tg = 3π/16Gρg are the SPH smoothing length, the sound speed, and the dynamical time of the gas respectively (Kawata 1999). When a gas particle is eligible to form stars, its star-formation rate (SFR) is where c * is a dimensionless SFR parameter and tg is the dynamical time, which is longer than the cooling timescale in the region eligible to form stars. This formula corresponds to the Schmidt law that SFR is proportional to ρ 1.5 g . We fix c * = 0.5, following K01b.
We assume that stars, which are represented by a star particle, are distributed according to the Salpeter (1955) Figure 2. Chemical yields as a function of age of a star particle with mass of 1 M ⊙ . The upper (lower) panel shows the total ejected oxygen (iron) mass. The thick (thin) line indicates the history of a star particle with the metallicity of Z = 0.02 (10 −4 ). initial mass function (IMF). The IMF by number, Φ(m), in each mass interval dm is defined as 2 where x = 1.35 is the Salpeter index and the coefficient A is determined by the normalization in the mass range M l ≤ m ≤ Mu. We set M l = 0.2 M ⊙ and Mu = 60 M ⊙ in this paper.

Feedback
GCD+ takes account of energy feedback and metal enrichment by SNe, and metal enrichment from intermediate mass stars. We consider here both SNe II and SNe Ia. The code calculates the event rates of SNe II and SNe Ia, and the yields of SNe II, SNe Ia, and intermediate mass stars for each star particle at every time step, considering the IMF and stellar lifetimes. Here, we assume the same stellar lifetimes as the ones used in Kodama (1997) and Kodama & Arimoto (1997). The simulation follows the evolution of the abundances of several chemical elements ( 1 H, 4 He, 12 C, 14 N, 16 O, 20 Ne, 24 Mg, 28 Si, 56 Fe, and Z, where Z means the total metallicity). We refer to Table 2 of  for the solar abundances.

Type II Supernovae
We assume that each massive star (≥ 8 M ⊙ ) explodes as a Type II supernova. To calculate the ejected mass of gas and metals by SNe II, we use the stellar yields derived by .  provide yields for a grid of stellar masses between 11 and 40 M ⊙ , and metallicities between 0 and 1 Z ⊙ . With linear interpolation of these grids, we obtain yields as a function of the stellar mass and metallicity.  considered three different models for m ≥ 30 M ⊙ , their A, B, and C models. These differ in the amount of energy imparted by the piston in their models at explosion initiation. Following Timmes, Woosley, & Weaver (1995), we use the 'B' model in this mass regime. For m ≥ 40 M ⊙ , we assume the same abundance ratios for the yields and the same mass fraction for the remnants as those for 40 M ⊙ star. For Z ≥ Z ⊙ , we use the following simple scaling, e.g. the yield of carbon from a star with the mass, m, and metallicity, Z, where mej,C is the total mass of ejected carbon, including newly synthesized and initial carbon; mej = mej,H + mej,He + mej,Z indicates the total ejected mass; Z C,⊙ is the solar abundance of carbon. The yields from stars with initial masses between 8 and 11 M ⊙ are still unclear (e.g. Hashimoto, Iwamoto, & Nomoto 1993). Hence, we simply interpolate linearly between the yields for the lowest mass, i.e. ∼ 11 M ⊙ , star in  and those for the highest mass, i.e. 8 M ⊙ star in van den Hoek & Groenewegen (1997), as mentioned in the next section.  Chieffi, et al. 2001;Marigo et al. 2001).

Type Ia Supernovae
We adopt the SNe Ia model proposed by KTN00, which is based on the SNe Ia progenitor scenario proposed by Hachisu, Kato, & Nomoto (1999). This model is different from the conventional one proposed by Greggio & Renzini (1983). Details of this model are described in KTN00. Here, we briefly explain this model and implementation in our code. Following KTN00, SNe Ia are assumed to occur at binary systems which consist of primary and companion stars with appropriate masses and metallicities ( M ⊙ . KTN00 designate the binary system whose companion has the mass within the former (latter) mass range "RG+WD (MS+WD) system". Finally, the total number of SNe Ia is obtained by the following equation as a function of age, t, of a star particle with the mass of ms M ⊙ , where mt is the mass of star whose lifetime is equal to t, and the mass function of the companion stars is assumed to be Φ d (m) ∝ m −1.35 by number, following KTN00. The term before the square bracket indicates the number of C+O WDs, i.e. primary stars. The first term within the square bracket determines how much fraction of C+O WDs evolves into SNe Ia from the MS+WD systems, and the second term determines how much fraction of C+O WDs evolves into SNe Ia from the RG+WD systems. Following KTN00, we set bMS = 0.05 and bRG = 0.02. The nucleosynthesis prescriptions for SNe Ia are taken from the W7 model of Iwamoto et al. (1999). Fig. 2 shows the total amount of oxygen and iron ejected from a star particle with the mass of 1 M ⊙ as a function of its age. Initially, metals are ejected only by SNe II. This continues until the 8 M ⊙ star dies (∼ 0.04 Gyr in the case of Z = 0.02). Oxygen is produced mainly by SNe II. After SNe II cease, the continuous ejection of oxygen is mainly due to the contribution from intermediate mass stars. On the other hand, a significant fraction of iron is produced by SNe Ia. Stars with Z = 0.02 lead to SNe Ia depending on the lifetime of their companion stars. SNe Ia occur in MS+WD (RG+WD) systems in the range of age of 0.7 − 1.5 (> 2.8) Gyr. However, a star particle with Z = 10 −4 does not lead to SNe Ia, because the above SNe Ia model restricts the metallicity range for progenitors of SNe Ia to [Fe/H] ≥ −1.1. For simplicity, we assume that the metallicity range for SNe Ia is log Z/Z ⊙ ≥ −1.1, instead of [Fe/H] ≥ −1.1. In the implementation, the code possesses a look-up table of the yields of all the chemical elements, remnant masses, number of SNe as a function of the age and metallicity, and calculates those values for each star particle at every time step.

Energy Feedback
One of the most difficult and most critical processes to model in galaxy formation simulations is the way in which the energy feedback from SNe affects the surrounding gas. Unfortunately, there is no clear understanding of how this should be modeled. We adopt a simple model proposed by Navarro & White (1993) and used in K01b. This model assumes that the energy produced by SNe affects only the temperature and velocity field of the surrounding gas, and its effect is implemented by increasing the thermal (E th ) and kinetic (E kin ) energy of the gas neighbours of each star particle by an amount corresponding to the energy released by SNe. We assume that a parameter fv = E kin /(E kin + E th ) defines the fraction of the available energy to perturb the gas velocity field, and the rest of the energy of the SNe contributes to the increase in the thermal energy of the gas (see K01b for details). It is known that kinetic feedback affects the history of star formation more strongly than thermal feedback, which quickly dissipates due to efficient radiative cooling where the gas density is high enough to form stars. The parameter, fv, controls the magnitude of the effect of SNe (K01b; Navarro & White 1993). We assume that each SN yields the energy of ǫSN ×10 51 ergs, and then E kin = fvǫSN 10 51 ergs. Because an initial SN energy has not been established quantitatively yet, we consider the available SN energy to be a free parameter. Finally, in our code, there are two parameters, fv and ǫSN , to control the magnitude of the effect of SNe.

Implementation
Based on the above feedback processes, the code calculates the amount of the mass, energy, and heavy elements released from each star particle within each time step. Our code adopts the simple feedback scheme suggested by Katz (1992) (see also Lia, Portinari, & Carraro 2002, for a proposed alternative scheme). The mass, energy, and heavy elements are smoothed over the neighbouring gas particles using the SPH smoothing algorithm. For example, when the i-th star particle ejects the mass of MSN,i , the increment of the mass of the j-th neighbour gas particle is given by where and W (x) is an SPH kernel, Here, rij is the distance between the i-th and j-th particles, and hi is the smoothing length for the i-th particle. Note that the above equations employ hi, instead of hij = (hi + hj)/2 which is adopted in other SPH calculations, such as density and pressure gradient (see Kawata 1999). This is because feedback is a one-way process, which does not have to be symmetrised. In addition, if all neighbour particles have much smaller smoothing length than the i-th particle, hij is possible to become smaller than rij , and then the i-th particle loses the place to feedback. The use of hi ensures to avoid this problem. We revise the update algorithm for the smoothing length in Kawata (1999) who uses the algorithm suggested by Hernquist & Katz (1989). The new algorithm is based on that suggested by Thacker et al. (2000). First, we count neighbour particles for the i-th particle with hi using the following smoothed kernel, where r is the distance of the neighbour particle from the i-th particle, and W (x) is the same spline kernel as equation (7) 3 . The number of neighbours, N nb,i , is no longer integer, but real number. This revision helps alleviate the discontinuity in the number of neighbours. Next, to avoid a sudden change in the smoothing length, we modify equation (14) of Kawata (1999) which determines the evolution of the smoothing length. Setting s = (Ns/N n nb,i ) 1/3 , equation (14) of Kawata (1999) is expressed as, where Ns = 40 is the desired number of neighbours, h n i means the smoothing length of the i-th particle at n step. Equation (14) of Kawata (1999) corresponds to the case of a = 0.5. Here we change a as a function of s, In addition, Thacker et al. (2000) use a modified Courant condition depending on the relative velocities to neighbour particles. Because this condition requires a significant amount of additional computational cost, once star formation and kinetic feedback are involved, we abandon this condition. As shown in Fig. 1 of Thacker et al. (2000), the above two algorithms lead to a dramatical improvement in keeping the number of neighbour particles constant, and the modified Courant condition adds a marginal effect.
We follow the evolution of the smoothing length of star particles until the end of simulation, because we relax the instantaneous recycling for SNe II and consider the feedback and chemical enrichment from SNe Ia and intermediate mass stars. Even in the above new algorithm, involving kinetic feedback enables gas particles to escape from the position of star particles more quickly than the speed of the increase in the smoothing length of the star particles. In serious cases, it happens that there is no neighbour particle, and the mass conservation is broken due to the missing of the mass deposition from stars to gas. To solve this problem, we introduce the iteration process using equation (9), once the number of neighbour particles is less than ten. Although this iteration requires an extra computational cost, it ensures complete mass conservation.

APPLICATION TO ELLIPTICAL GALAXY FORMATION
To see how the new code works, in this section, we apply GCD+ to simulations of elliptical galaxy formation in the semi-cosmological model described in K01b. K01b examined the mass-dependent properties of elliptical galaxies, such as the CMR, the Kormendy relation, and the [Mg/Fe]magnitude relation. Our new code relaxes some of the simple assumptions in K01b, and becomes a more sophisticated code. Hence, the re-examination of those properties is an interesting application as well as an important test for GCD+.

Elliptical Galaxy Formation Model
Following K01b, we consider an isolated sphere, as a seed galaxy, on which small-scale density fluctuations corresponding to a cold dark matter (CDM) power spectrum are superimposed. Here, we use Bertschinger's software COS-MICS (Bertschinger 1995) to generate initial density fluctuations. To incorporate the effects of fluctuations with longer wavelengths, the density of the sphere has been enhanced and a rigid rotation corresponding to a spin parameter, λ, has been added. The initial conditions of this model are determined by the following four parameters: λ, Mtot, σ8,in, and zc. The spin parameter is defined by where J is the total angular momentum of the system, E is the total energy, and Mtot is the total mass of the sphere, which is composed of dark matter and gas; σ8,in is the rms mass fluctuation in a sphere of radius 8 h −1 Mpc, which normalizes the amplitude of the CDM power spectrum; zc is the expected collapse redshift. If the top-hat density perturbation has an amplitude of δi at the initial redshift, zi, we ob-tain zc = 0.36δi(1+zi)−1 approximately (e.g. Padmanabhan 1993). Thus, when zc is given, δi at zi is determined. Kawata (1999) found that the seed galaxy which has a slow rotation corresponding to λ = 0.02 and small-scale density fluctuations evolves into an elliptical-like system. Thus, we employ λ = 0.02. This spin parameter is close to the minimum value possible in a CDM universe, according to the results of Nbody simulations (Barnes & Efstathiou 1987;Warren et al. 1992). We fix σ8,in = 0.5 and zc = 3.5. Our simulations assume a flat universe (Ω = 1) with a baryon fraction of Ω b = 0.1 and a Hubble constant of H0 = 50 km s −1 Mpc −1 . We carry out each simulation using 9171 particles for gas and dark matter, respectively. We simulate the evolution of each model from zi = 40 to z = 0.
The morphological evolution of all the models which are shown in this paper are similar to the evolution seen in Fig. 1 of Kawata (1999). A nearly spherical stellar system is formed at z = 0 in all the models. Following K01b, we focus on the chemical and photometric properties for the simulation end-products at z = 0. In our simulations the stellar particles each carry their own age and metallicity "tag", which, when combined with population synthesis, enables us to derive the photometric properties of the simulated stellar systems. The photometric properties are derived by the same analysis as K01b. Details of this analysis are described in Section 3.2 of Kawata (2001a) and Section 4 of K01b.
In this analysis, the spectral energy distribution (SED) of each stellar particle is assumed to be that of a single stellar population (SSP) that means a coeval and chemically homogeneous assembly of stars. Since the observational data with which our results should be compared provide the luminosity distribution projected to a plane, we have to derive the projected distribution of SED from the three dimensional distribution of stellar particles. Finally, we obtain the x-y projected images as shown in Figure 5 of Kawata (2001a), when z-axis is set to be the initial rotation axis. We confirmed that the results do not depend on the direction of the projection, because the end-products are nearly spherical. The flux of each stellar particle is smoothed using a Gaussian filter with the filter scale of 1/4 of the softening length of the stellar particle. These images provide similar information to the imaging data obtained in actual observations. Thus, we can obtain various photometric properties from these images in the same way as in the analysis of observational imaging data. We use the images similar to the one displayed in Figure 5 of Kawata (2001a), but employing a 1001 × 1001 pixel mesh to span the squared region with 100 kpc on a side.
We adopt the data of SSPs of Kodama & Arimoto 97 model (Kodama 1997;Kodama & Arimoto 1997) The global photometric properties, such as the total luminosity, the colours, and the effective radius, are obtained from the projected image data. Following K01b, the total magnitude and the effective radius are derived by fitting the surface brightness profile to the r 1/n law (eq. [11] of K01b). We set the center to the position of a pixel which has the maximum V band luminosity.
These procedures enable us to compare the properties for simulation end-products with those for observed elliptical galaxies directly and quantitatively. In the following sections, we begin by showing the best model which can reproduce the CMR of observed elliptical galaxies, and examine what is responsible for the CMR based on the detailed analysis of the simulation results. Next, we discuss the physical sizes of the simulation end-products, comparing them with the observed Kormendy relation. Finally, we examine the correlation between the total magnitude and the abundance ratio of Mg to Fe, i.e. the [Mg/Fe]-magnitude relation, which is sensitive to the star formation history.

The Best Model: the Colour-Magnitude Relation
As explained in Section 2, due to the lack of knowledge of the physics of star formation and SNe feedback, our code has a number of free parameters, such as the mass range  of the IMF (M l and Mu in Section 2.2) and the feedback parameters of fv and ǫSN in Section 2.3.4. We fix these parameters so that the observed CMR is reproduced. To compare with the CMR, we simulate the evolution of seed galaxies with three different masses of Mtot = 4 × 10 12 , 8 × 10 11 , and 2 × 10 11 M ⊙ (hereafter, models HM, MM, and LM, respectively). The mass and spatial resolutions of each model are summarized in Table 1.
The mass range of the IMF controls yields per mass of new born star particles. Higher Mu or M l leads to larger yields, and makes the metallicity of the end-products higher. We chose (Mu, M l ) = (60, 0.2) M ⊙ , which produces enough yields for the highest mass model (HM) to reproduce the observed red colour of galaxies with similar luminosities. As shown in K01b, we also find that the stronger feedback causes the steeper slope in the CMR. The magnitude of the feedback effect is controlled by the parameters ǫSN and fv. We fix the parameter set of (ǫSN , fv) = (0.1, 0.002) such that the observed CMR is reproduced, when the above IMF is adopted. According to high resolution 1D simulations of an SN remnant in Thornton et al. (1998), 90 % of the initial SN energy is lost in radiation during its early expansion phase; this is not resolved in our simulations. The available SN energy of 10 50 ergs which we adopt is consistent with the result of Thornton et al. (1998), when a canonical initial SN energy of 10 51 ergs is assumed. Thornton et al. (1998) also suggested that in the last stage of their simulation about ∼ 8.5 × 10 49 erg, i.e. 85 % of the available energy, is found in the form of kinetic energy. This corresponds to fv = 0.85 in our definition, and is much larger than that we assumed here. However, the resolution limit of our simulation (a few kpc) is much larger than the physical scale of their simulation (∼0.15 kpc). Unfortunately, little is know about the fraction of the kinetic energy which is still available to affect the interstellar medium on such large scales. Therefore, fv suggested by Thornton et al. (1998) is considered to be an upper limit, which does not invalidate the lower value which we adopted here. Although this parameter set is not a unique solution, it is one which best explains the observed CMR. Therefore, the set of models HM, MM, and LM with the above parameters is called the "best model" in this paper. Fig. 3 shows the CMR of the best model which reproduces the observed CMR for galaxies in the Coma Cluster very well. The data for galaxies in the Coma Cluster are the observed CMR of Bower, Lucey, & Ellis (1992a). Since there is no difference between S0s and ellipticals in the scaling relations which we now discuss, we do not distinguish between S0s and ellipticals. Bower, Lucey, & Ellis (1992a) supply U − V and V − K colours which refer to an aperture size of 11 arcsec, and the V band total magnitude de-   Here, the ejected gas is defined as the sum of all the gas particles whose galactocentric radius is greater than 20 kpc.
rived from a combination of their data and the literature. Throughout this paper we adopt the distance modulus of the Coma cluster of m − M = 34.7 mag; the Virgo distance modulus of m − M = 31.01 (Graham et al. 1999) and the relative distance modulus of the Coma with respect to the Virgo of m − M = 3.69 (Bower, Lucey, & Ellis 1992b). This gives the luminosity distance of 87.1 Mpc for the Coma. We assume that the angular diameter distance equals the luminosity distance, because the redshift of the Coma cluster (z ∼ 0.023) is nearly zero cosmologically. Then the aperture size of 11 arcsec at the distance of the Coma Cluster corresponds to ∼ 5 kpc. The colours for simulation end-products are derived in the same aperture size as that used in the observations. K01b claimed that the aperture effect should not be ignored, when discussing the CMR observed in an aperture of a small size. Observed elliptical galaxies (especially in luminous galaxies) have colour gradients which make the colour at the centre redder than that in the outer regions (e.g. Peletier et al. 1990). The end-products in our simulations for modes HM and MM also have significant colour gradients as seen in Table 2. Since colours within a fixed aperture yield colours in a more central region for larger galaxies, it is possible that the colours of large galaxies become redder than those of small galaxies, even if the mean colour of the whole galaxy is the same between the large and small galaxies. Fig. 4 shows the CMR for the mean colour within the aperture of 99 kpc, which covers almost the entire galaxy in all the models which we have examined. As expected, the larger galaxy's colours are changed more dramatically by the aperture, which makes the slope between models HM and MM shallower than that between models MM and LM. The aperture effect contributes significantly to the slope of the CMR in Fig. 3 in the mass range between models HM and MM (see also Kodama et al. 1998).
As seen in K01b, Table 2 shows that the colour and metallicity gradients become shallower, and Sersic law index n (see eq. [11] of K01b) in the surface brightness profile becomes smaller, with decreasing mass of the system. It may appear strange that the mean colour within the 99 kpc aperture is bluer than that within the 5 kpc aperture in model LM, although the colour gradient for model LM is flat. This is because Table 2 shows the gradients within the radius where the B-band magnitude is brighter than µB = 24.5 mag arcsec −2 , following K01b, and the gradient at larger radii becomes negative.
The sequences of colours and line strengths among elliptical galaxies can almost equally well be attributed to either age difference or metallicity differences (Worthey 1994). In our simulation this degeneracy can be broken completely, because we can directly analyse the metallicity and age in the simulation end-products. Fig. 5 shows the metallicity and age for each model. Here, the metallicities and ages are luminosity weighted values. Although the ages are all very old irrespective of the models and the apertures, the metallicities exhibit similar behavior to the colours as a function of luminosity. Thus, we conclude that the slope of the CMR Figure 8. The density-temperature distribution of the gas particles with higher (left), intermediate (middle), and lower (right) metallicity at z = 3.27 for models HM (upper panels) and LM (lower panels). The solid curves separate the region where the cooling times are shorter (upper region) and longer than the age of the universe at z = 3.27. of the best model in Fig. 3 is caused solely by the changes in metallicity. Fig. 6 shows the histories of SNe II and SNe Ia event rates normalized by the total mass of the system. In our simulations, the histories of SNe II events roughly trace those of star formation, because the lifetime of SNe II progenitors is shorter than the size of the bins of the histograms. We can clearly see that star formation ceases abruptly around t = 1.7 Gyr (z = 1.7), regardless of the mass of the system. We confirmed that gas particles begin to blow out from the stellar system, i.e. the galactic wind occurs when star formation stops. Subsequently, all the gas particles overcome the binding energy of the dark matter and stars, and escape from the system, as shown in Fig. 9 of K01b. Fig. 7 shows the mass fraction of the ejected gas as a function of the mass of the system. The ejected gas mass in Fig. 7 is defined as the total mass of all the gas particles whose galactocentric radius is greater than 20 kpc at z = 0. Comparing between models MM and LM, the lower mass system ejects a greater fraction of gas. Since the ejected gas cannot contribute to further metal enrichment, metal enrichment is more strongly suppressed in the lower mass system, leading to its lower metallicity in Fig. 5. On the other hand, between models HM and MM, the lower mass system ejects a slightly smaller fraction of the gas. In this luminosity range, the aperture effect appears to be a dominant factor in explaining the slope of the mass-metallicity relation. However, there is a significant slope between models HM and MM in mean metallicity within the large aperture (Fig. 5). We find that this is because there is another factor which contributes to the mass-metallicity relation. Fig. 8 plots density against temperature for the gas particles at z = 3.27 for models HM and LM. At z = 3.27 the system is almost relaxed, the gas particles which stay in the lower right region in the panel are hot gas which has not yet cooled since they virialised. On the other hand, the gas particles which stay in the upper left (log ρgas ≥ −28 and log T ≤ 4) region represent the cold gas. The gas particles in the lower left region are escaping from the system. Comparing between models, the high mass system has more hot gas and less cold gas than the low mass system. The solid lines in the panels indicate where the cooling time equals the age of the universe at z = 3.27. Cooling is efficient in the upper region from this line, and inefficient in the lower region from this line. In the high mass system, the hot gas particles fall into where cooling is inefficient. Therefore, a large fraction of the gas cannot cool or infall into the central region where stars are forming, thus the hot gas is not affected by metal enrichment from stars (Fig. 8). Since the amount of gas is small in the central region, relative to the amount of metals returned from stars, the metallicity of the gas in the central region increases rapidly. In other words, a low infall rate of the gas induced by inefficient cooling leads to efficient chemical enrichment. This is a similar mechanism to the infall model which was proposed to solve the G-dwarf problem of the closed box chemical evolution model for the solar vicinity (e.g. Pagel 1997). On the other hand, in the low mass system the hot gas stays close to the line in Fig. 8. A large fraction of this gas which cools efficiently infalls into the central region, and suffers from metal enrichment. Since there is a large amount of gas in the star forming region, the metallicity of the gas in the central region slowly increases. In other words, a high infall rate of gas, caused by efficient cooling, leads to inefficient chemical enrichment. Consequently, the high mass system makes metal rich stars more quickly than the low mass system (Fig. 9), and the mean metallicity of the high mass system becomes higher than that of the low mass system. Hence, we conclude that not only the galactic wind, which ejects a larger fraction of gas in the lower mass system, but also the mass dependence of the gas infall rate, contributes to the slope in the mass-metallicity relation (Fig. 5) and the CMR (Fig. 3).

The Role of SNe Ia
To assess the relative importance of SNe Ia to SNe II during the evolution of elliptical galaxies, we have run an additional model which is similar to model LM, except that SNe Ia are not included. Hereafter, this model is called "noIa-LM". Figs. 3-7 also show results of model noIa-LM. Comparison between models LM and noIa-LM demonstrates that SNe Ia are crucial sources for the galactic wind. For example, Fig. 6 shows that model noIa-LM does not have a clear cessation of star formation. Consequently, model noIa-LM ejects a smaller fraction of gas from the system (Fig. 7), and has a higher metallicity and a redder colour than model LM (Figs. 3 and 5). In the U − V and MV diagram, model noIa-LM appears to follow the CMR of the best model. This is due to the younger age of model noIa-LM as seen in the lower panel of Fig. 5. Model noIa-LM adopts the same feedback parameters as those of model LM. In addition, model noIa-LM has a larger number of SNe than model LM, because more stars have formed. Thus, this result indicates that the galactic wind is caused mainly by SNe Ia. It is also notable that the event rate of SNe Ia exceeds that of SNe II, when the galactic wind occurs in models HM, MM, and LM. We therefore conclude that SNe Ia play a crucial role in driving and maintaining the galactic wind, and thus in contributing to the evolution of elliptical galaxies. This conclusion is consistent with K01b's suggestion, and our sophisticated treatment of SNe II and SNe Ia reinforces this idea.

Instantaneous recycling SNe II and burst SNe
Ia model K01b adopted instantaneous recycling for SNe II and a simple burst SNe Ia model, in which all the SNe Ia occur simultaneously 1.5 Gyr after star formation. To compare the feedback effect of this simple description of SNe in K01b with that of the more realistic SNe model used in this paper, we now examine another model which employs instantaneous recycling SNe II and 1.5 Gyr delay burst SNe Ia. The total number of SNe II is given by the integral of number of stars with the mass higher than 8 M ⊙ . The total number of SNe Ia expected in each star particle is assumed to be given by the sum of the MS+WD and RG+WD binary systems, which corresponds to the integral of equation (4) from t = 0 to t = ∞. We also assume that SNe Ia occur only in star particles with log Z/Z ⊙ ≥ −1.1, following the KTN00 model. Using this SNe model, we carry out the simulation with the same total mass as model LM. Hereafter, this model is called "iIIbIa-LM". Fig. 6 shows that the galactic wind ceases star formation in this model, although the amount of ejected gas is much smaller than that for model LM (Fig. 7). As a result, model iIIbIa-LM has a higher metallicity and a redder colour than model LM (Figs. 3 and 5). In both models LM and iIIbIa-LM, the galactic wind occurs just after SNe Ia ignite. However, SNe Ia in model iIIbIa-LM begin to occur later than those in model LM, due to the simple 1.5 Gyr delay, and model iIIbIa-LM has a longer duration of star formation than model LM. As a result, in model iIIbIa-LM the chemical enrichment progresses for a longer time, which leads to the higher metallicity and redder colour. Thus, the final properties of elliptical galaxies are sensitive to the time delay of SNe Ia, which is the main trigger of the galactic wind. It is also worth noting that the peak rate of star formation in model iIIbIa-LM is significantly higher than that in model LM (Fig. 6). This means that the suppression of star formation by continuous SNe II feedback is stronger than that from instantaneous recycling SNe II. These results demonstrate that we should not ignore the lifetime of both SNe II and SNe Ia progenitors, in studying the dynamical evolution of elliptical galaxies.

The Kormendy Relation
We now examine the physical size of the simulation endproducts for all the models presented above. Fig. 10 shows a comparison of the Kormendy relations for the simulation end-products and the Coma Cluster galaxies both in the V and K bands. We refer to the data of the Coma Cluster galaxies of Pahre (1999). In deriving the absolute magnitude and effective radius in the kpc from the data set in Pahre (1999), we assume the same distance modulus as mentioned in Section 3.2. The best model shows a tendency for higher mass galaxies to have larger effective radii, which is qualitatively consistent with the observed Kormendy relation. However, the effective radii are systematically larger and the slope is slightly shallower than the observational data. In comparison among the low mass systems, model LM shows a significantly larger effective radius than models noIa-LM and iIIbIa-LM. Sections 3.3 and 3.4 show that the effect of SNe on star formation in the models decreases as: LM, iIIbIa-LM, and noIa-LM, and the effective radius decreases in the same order. In other words, stronger feedback causes an expansion of the system. Due to this expansion effect, the best model fails to explain the observed Kormendy relation, which is the same problem as that found in K01b (see also Chiosi & Carraro 2002).

The [Mg/Fe]-Magnitude Relation
The abundance ratio of [Mg/Fe] of elliptical galaxies has been well-studied by several groups ( and SNe Ia, respectively, and because SNe Ia have a longer delay than SNe II after the formation of stars, this correlation provides a strong constraint on the star formation history of elliptical galaxies (Matteucci, Ponzone, & Gibson 1998;Thomas, Greggio, & Bender 1999;Romano et al. 2002;Tantalo & Chiosi 2002). In Fig. 11, we compare [E/Fe] for the simulation end-products of all the models with those derived by Trager et al. (2000a). We use the absolute B band magnitudes in Trager et al. (2000b). In both the data of Trager et al. (2000a)  Model HM appears to have larger fraction of those stars than model LM. In addition, model LM stops star formation before enrichment by SNe II whose yields have [Mg/Fe]∼ 0.1 becomes dominant. These two effects lead to the slope of the best model seen in Fig. 11. The reason why [Mg/Fe] of model noIa-LM is smaller than that of model iIIbIa-LM is a little more complicated -it is due to the enrichment from intermediate mass stars. In theory, with the exception of carbon, nitrogen, and oxygen, the abundance pattern of the yields for intermediate mass stars is the same as the abundance pattern when stars are born, i.e. the initial abundance pattern. However, the abundance pattern for intermediate mass stars is set to be the solar abundance pattern in our code, due to the reason mentioned in Section 2.3.2. As time progresses, the yields for intermediate mass stars becomes more important, because more stars with high metallicities finish the main-sequence phase. Therefore, model noIa-LM which has the longest duration of star formation makes a significant amount of stars which suffer from the yields with the "artificial" solar abundance pattern. These stars are seen in Fig. 12 as stars with high metallicity ([Fe/H]> 0.2) and slightly lower [Mg/Fe] (around 0.05). As a result, model noIa-LM has lower Figure 11. Comparison of the simulation end-products with observed early-type galaxies (Trager et al. 2000a)

Discussion
We have applied our new chemo-dynamical evolution code, GCD+, to simulations of elliptical galaxy formation, and confirmed the results of K01b for the CMR, the Kormendy relation, and the [Mg/Fe]-magnitude relation. Tables 2 and  3 summarise the photometric properties and stellar population for all the models presented. The most fundamental difference between our new code and the one used in K01b is its taking into account of the lifetime of stars and the adoption of the new SNe Ia model proposed by KTN00. We have obtained a more realistic history of SNe Ia rate than K01b (Fig. 6). The history of the SNe Ia rate follows the SNe Ia progenitor model of KTN00. In the "best model" where star formation is stopped by galactic winds at high redshift, the histories of the SNe Ia rate have two peaks. These two peaks are caused by SNe Ia from the MS+WD and RG+WD systems (Section 2.3.3). Their evolution is similar to that shown in Fig. 5 of KTN00. After the second peak, this model predicts the continuous explosion of SNe Ia from the RG+WD system until z = 0. Capellaro, Evans, & Turatto (1999) provides the observed SNe rates at nearby galaxies. This SNe II rate for elliptical galaxies gives an upper limit, and is consistent with the best model where there is no SNe II. Points with error bars in Fig. 6 show their observed SNe Ia rates which are transfered from units of SNu to units used in Fig. 6 based on the B band luminosity for each model. In the best model SNe Ia rates are in good agreement with the observed rates. This result which comes from a self-consistent chemo-dynamical evolution provides strong support for the KTN00 model which is based on pure chemical evolution studies.
Our sophisticated SNe Ia model has clarified that the feedback effect of SNe Ia plays a crucial role in the evolution of elliptical galaxies rather than that of SNe II. The lifetimes of progenitor stars of SNe II are much shorter than those of SNe Ia. As shown in Fig. 2, SNe II occur immediately (between 4 and 40 Myr for stars with Z=0.02) after star formation events. On the other hand, there is a significant delay (> 0.7 Gyr for stars with Z=0.02) between star formation and SNe Ia. Therefore, SNe II occur in dense gaseous environments still actively forming stars, whereas SNe Ia act on tenuous gas left after star formation. This difference makes SNe Ia the major trigger of galactic winds. Comparison between models noIa-LM and LM shows this clearly in Section 3.3. Also model iIIbIa-LM demonstrates that the simplified model for SNe II and SNe Ia underestimates significantly the SNe effect on the dynamical evolution of elliptical galaxies.
The difference in SNe II and Ia model between K01b and this paper is one of the reasons why our new model requires only 10 50 ergs and fv = 0.002 as an SN feedback energy to reproduce the observed CMR, whereas K01b requires 4 × 10 51 ergs and fv = 0.9. This difference is caused by the combination of some physical processes, such as the IMF, SNe II and SNe Ia feedback, and the cooling function, which are updated in the new code. Comparing the results, our new code requires a more reasonable amount of the SN feedback energy than that in K01b, to reproduce the observed CMR. However, we might still ignore some physical processes which contribute to the CMR, and be able to further reduce the feedback energy. For example, the UV background radiation is a possible candidate as such physical process, as suggested by K01b. A recent semi-analytical study demonstrates that the UV background radiation has a similar effect to the SNe feedback on the CMR (Nagashima & Gouda 2001).
We find that mass dependences of both the galactic wind and the gas infall rate are responsible for the slope of the CMR. The conventional galactic wind scenario considers only the former mechanism. Our numerical simulation shows that the higher mass system has a higher virial temperature and lower efficiency of radiative cooling. Therefore, the gas infall rate decreases with increasing mass of the system. This lower gas infall rate leads to the increase in the fraction of the high metallicity stars (Fig 9). As a result, the mean metallicity of the higher mass system becomes higher than that of the lower mass system. This is a novel mechanism contributing to the CMR. Also in the [Mg/Fe]-magnitude relation, we reach the same conclusion as K01b qualitatively. However, the value of [Mg/Fe] is systematically different from that in K01b. This is due to the difference in the stellar yields of SNe II. We use the yields calculated by , whereas K01b adopts those by Notomo et al. (1997). Although  and Notomo et al. (1997) are two representative computations for SNe II nucleosynthesis, it is known that the amount of iron yield derived in  is significantly larger than that in Notomo et al. (1997), and therefore the yields of  lead to smaller [Mg/Fe] (Timmes, Woosley, & Weaver 1995;Gibson 1997;Gibson, Lowenstein, & Mushtzky 1997;Thomas, Greggio, & Bender 1998). The theoretical stellar yields still have a large uncertainty, and to discuss the zero-points of [Mg/Fe] for simulation endproducts we need a more reliable yield model. Nevertheless, we can discuss the observed slope of the [Mg/Fe]-magnitude relation, and our model so far cannot explain this.
To explain the observed slope, it is a plausible scenario that low mass systems have a longer duration of star formation and are enriched by SNe Ia more than high mass systems (Faber, Worthey, & González 1992;Matteucci 1994;Matteucci, Ponzone, & Gibson 1998). K01b suggests that the UV background radiation is a possible candidate to realise the scenario, because it suppresses cooling and star formation more strongly in lower mass systems (Efstathiou 1992), and is expected to prolong the duration of star formation. Since we have developed a sophisticated chemo-dynamical evolution code, it is important to introduce the UV background radiation into GCD+, and examine the effect of the UV background radiation on galaxy evolution. In a forthcoming paper, we hope to throw a light on this issue.

CONCLUSION
We have developed a new chemo-dynamical evolution code called GCD+. The code includes both SNe II and SNe Ia modeling taking into account the lifetime of progenitor stars, and the chemical enrichment from intermediate mass stars. In this paper we have applied the code to simulations of elliptical galaxy formation described as the evolution of an isolated seed galaxy. We have shown that GCD+ is a useful and unique tool which enables us to compare simulation results with the observational data directly and quantitatively, and gives us new insights into the effects of dynamics on the chemical evolution of galaxies. Encouraged by the success of this study, we plan to apply GCD+ to high-resolution cosmological simulations in the near future. In the real universe, even isolated galaxies at z = 0 are expected to experience interactions with the another galaxies in the past. In fact, Gnedin (1998) suggests that interactions between galaxies are a dominant mechanism to transport metals from the galaxies to intergalactic medium rather than SNe feedback. Therefore, using simulations covering a cosmological scale, it will be valuable to examine what physical process, or combination of processes, induces the tight observed CMR. In addition, cosmological simulations using GCD+ can provide detailed metallicity distributions not only in the stellar component of individual galaxies, but also in the intergalactic medium, such as Lyman-α clouds and the hot gas in clusters of galaxies. As mentioned in Section 2.3, the different chemical elements have progenitors with different masses, and different mass stars have different lifetimes. It is highly probable that the abundance ratios are fairly sensitive to the metal transport mechanism and its epoch. New observational facilities will offer a wealth of information about the distribution of metallicity, as well as the abundance ratio of each element, such as [N/O], [Mg/Fe] and [Si/Fe], with unprecedented accuracy. Studies comparing new observations with our future cosmological simulation results will show clearly which metal transport mechanism is dominant, and thus contribute to understanding the formation and evolution of galaxies as well as clusters of galaxies from a new perspective.