format: text

text:
# GCD+: A New Chemodynamical Approach to Modeling Supernovae and Chemical Enrichment in Elliptical Galaxies

## 1 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$`\stackrel{`}{\mathrm{e}}`$, 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 & 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 self-consistently (e.g. Carraro, Lia, & Chiosi, 1998; Mosconi et al., 2001; Chiosi & Carraro, 2002). Three dimensional chemo-dynamical 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, 2001). 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, 1999, 2000, 2001) 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 colour-magnitude 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.23.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.

## 2 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).

### 2.1 Gas Cooling

To model gas cooling, we use the cooling curves computed by MAPPINGS III111MAPPINGS III is the successor of MAPPINGS II described in Sutherland & Dopita (1993), and available at http://www.mso.anu.edu.au/ralph/map.html written by R.S. Sutherland. Using MAPPINGS 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 $`\mathrm{\Lambda }(T,[\mathrm{Fe}/\mathrm{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 $`[\mathrm{Fe}/\mathrm{H}]=3`$ is the same as the one with $`[\mathrm{Fe}/\mathrm{H}]=3`$, and the cooling rate of the gas with $`[\mathrm{Fe}/\mathrm{H}]>0.5`$ is equivalent to that at $`[\mathrm{Fe}/\mathrm{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 $`T_i=P_i\mu m_p/(k_\mathrm{B}\rho _i)`$, where $`P_i`$ and $`\rho _i`$, are the pressure and density for $`i`$-th particle, and $`\mu `$, $`k_\mathrm{B}`$, and $`m_p`$ are the mean molecular weight, Boltzmann’s constant, and the proton mass, respectively. For simplicity, we fix $`\mu =0.6`$, regardless of the metallicity. We set the lower limit of the temperature to be $`T_{\mathrm{lim}}=100`$ K.

### 2.2 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, $`\rho _{\mathrm{crit}}=2\times 10^{25}\mathrm{g}\mathrm{cm}^3`$, i.e. $`n_\mathrm{H}0.1\mathrm{cm}^3`$; 2) the gas velocity field is convergent, $`𝒗_i<0`$; and 3) the Jeans unstable condition, $`h/c_s>t_\mathrm{g}`$, is satisfied, here $`h`$, $`c_s`$, and $`t_\mathrm{g}=\sqrt{3\pi /16G\rho _\mathrm{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

$$\frac{d\rho _{}}{dt}=\frac{d\rho _\mathrm{g}}{dt}=\frac{c_{}\rho _\mathrm{g}}{t_\mathrm{g}},$$
(1)
where $`c_{}`$ is a dimensionless SFR parameter and $`t_\mathrm{g}`$ 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 $`\rho _\mathrm{g}^{1.5}`$. 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) initial mass function (IMF). The IMF by number, $`\mathrm{\Phi }(m)`$, in each mass interval $`dm`$ is defined as222Equation (2) corresponds to an IMF by mass $`m^x`$

$$\mathrm{\Phi }(m)dm=Am^{(1+x)}dm,$$
(2)
where $`x=1.35`$ is the Salpeter index and the coefficient $`A`$ is determined by the normalization in the mass range $`M_\mathrm{l}mM_\mathrm{u}`$. We set $`M_\mathrm{l}=0.2\mathrm{M}_{\mathrm{}}`$ and $`M_\mathrm{u}=60\mathrm{M}_{\mathrm{}}`$ in this paper.

### 2.3 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 (1H, 4He,12C, 14N, 16O, 20Ne, 24Mg, 28Si, 56Fe, and Z, where Z means the total metallicity). We refer to Table 2 of Woosley & Weaver (1995) for the solar abundances.

#### 2.3.1 Type II Supernovae

We assume that each massive star ($`8\mathrm{M}_{\mathrm{}}`$) 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 Woosley & Weaver (1995). Woosley & Weaver (1995) 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. Woosley & Weaver (1995) considered three different models for $`m30`$ 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 $`m40`$ 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 $`ZZ_{\mathrm{}}`$, we use the following simple scaling, e.g. the yield of carbon from a star with the mass, $`m`$, and metallicity, $`Z`$,

$`m_{ej,C}(Z,m)`$ $`=`$ $`m_{ej,C}(Z_{\mathrm{}},m)+m_{ej}(Z_{\mathrm{}},m)Z_{C,\mathrm{}}`$ (3)
$`\times (Z/Z_{\mathrm{}}1),`$
where $`m_{ej,C}`$ is the total mass of ejected carbon, including newly synthesized and initial carbon; $`m_{ej}=m_{ej,H}+m_{ej,He}+m_{ej,Z}`$ indicates the total ejected mass; $`Z_{C,\mathrm{}}`$ 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 Woosley & Weaver (1995) and those for the highest mass, i.e. 8 M star in van den Hoek & Groenewegen (1997), as mentioned in the next section.

#### 2.3.2 Intermediate Mass Stars

The yields and remnant masses for intermediate mass stars ($`8`$ M) are taken from van den Hoek & Groenewegen (1997). van den Hoek & Groenewegen (1997) present the yields of 1H, 4He, 12C, 13C, 14N, and 16O for stars with initial masses between 0.8 and 8 M and initial metallicities of $`Z=0.0010.04`$. We do not use 13C yields. van den Hoek & Groenewegen (1997) provide the newly formed and ejected metals unlike Woosley & Weaver (1995) who present the total ejected metals. Therefore, we calculate the total ejected metals using Table 1 of van den Hoek & Groenewegen (1997) for H, He, C, N, and O and the initial abundances for the other elements, i.e. Ne, Mg, Si, and Fe, which are simply scaled to the solar abundance set. As for stars with $`Z>0.04`$, we use the same scaling as equation (3) based on $`Z=0.04`$ yields. For $`Z<0.001`$, we simply assume the same yields as $`Z=0.001`$ yields. Whilst this simple assumption may overestimate the ejected metals, since the metals are ejected from intermediate mass stars later than those from SNe II, which eject much larger amounts of metals, this simplification should not affect the final metallicity. Nevertheless, the knowledge of the yields for low and zero metallicity stars is important in studies of the chemical composition of extremely metal poor objects observed at both low and high redshifts, although it is not yet well-established (e.g. Fujimoto, Ikeda, & Iben, 2000; Chieffi, et al., 2001; Marigo et al., 2001).

#### 2.3.3 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 ($`[\mathrm{Fe}/\mathrm{H}]1.1`$). The primary stars have the main-sequence mass in the range of $`m_{p,l}=3`$ M and $`m_{p,u}=8`$ M, and evolve into the C+O white dwarf (WD). The mass range of companion stars are restricted to between $`m_{d,RG,l}=0.9`$ M and $`m_{d,RG,u}=1.5`$ M and between $`m_{d,MS,l}=1.8`$ M and $`m_{d,MS,u}=2.6`$ 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 $`m_s`$ M,

$`N_{\mathrm{SNeIa}}(t)`$ $`=`$ $`m_s{\displaystyle _{m_{p,l}}^{m_{p,u}}}m^{(1+x)}𝑑m\left\{{\displaystyle _{M_l}^{M_u}}m^x𝑑m\right\}^1`$ (4)
$`\times [b_{\mathrm{MS}}{\displaystyle \frac{_{\mathrm{max}(m_{d,MS,l},m_t)}^{m_{d,MS,u}}\mathrm{\Phi }_d(m)𝑑m}{_{m_{d,MS,l}}^{m_{d,MS,u}}\mathrm{\Phi }_d(m)𝑑m}}`$
$`+b_{\mathrm{RG}}{\displaystyle \frac{_{\mathrm{max}(m_{d,RG,l},m_t)}^{m_{d,RG,u}}\mathrm{\Phi }_d(m)𝑑m}{_{m_{d,RG,l}}^{m_{d,RG,u}}\mathrm{\Phi }_d(m)𝑑m}}],`$
where $`m_t`$ is the mass of star whose lifetime is equal to $`t`$, and the mass function of the companion stars is assumed to be $`\mathrm{\Phi }_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 $`b_{\mathrm{MS}}=0.05`$ and $`b_{\mathrm{RG}}=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.71.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 $`[\mathrm{Fe}/\mathrm{H}]1.1`$. For simplicity, we assume that the metallicity range for SNe Ia is $`\mathrm{log}Z/Z_{\mathrm{}}1.1`$, instead of $`[\mathrm{Fe}/\mathrm{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.

#### 2.3.4 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_{\mathrm{th}}`$) and kinetic ($`E_{\mathrm{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 $`f_v=E_{\mathrm{kin}}/(E_{\mathrm{kin}}+E_{\mathrm{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, $`f_v`$, controls the magnitude of the effect of SNe (K01b; Navarro & White, 1993). We assume that each SN yields the energy of $`ϵ_{SN}\times 10^{51}`$ ergs, and then $`E_{\mathrm{kin}}=f_vϵ_{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, $`f_v`$ and $`ϵ_{SN}`$, to control the magnitude of the effect of SNe.

#### 2.3.5 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 $`M_{\mathrm{SN},i}`$, the increment of the mass of the $`j`$-th neighbour gas particle is given by

$$\mathrm{\Delta }M_{\mathrm{SN},j}=\frac{m_j}{\rho _{\mathrm{g},i}}M_{\mathrm{SN},i}W(r_{ij}/h_i),$$
(5)
where

$$\rho _{\mathrm{g},i}=\rho _\mathrm{g}(𝒙_i)=\underset{ji}{}m_jW(r_{ij}/h_i),$$
(6)
and $`W(x)`$ is an SPH kernel,

$$W(x)=\frac{8}{\pi h^3}\{\begin{array}{cc}16x^2+6x^3& \mathrm{if}0x1/2,\\ 2(1x)^3& \mathrm{if}1/2x1,\\ 0& \mathrm{otherwise}.\end{array}$$
(7)
Here, $`r_{ij}`$ is the distance between the $`i`$-th and $`j`$-th particles, and $`h_i`$ is the smoothing length for the $`i`$-th particle. Note that the above equations employ $`h_i`$, instead of $`h_{ij}=(h_i+h_j)/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, $`h_{ij}`$ is possible to become smaller than $`r_{ij}`$, and then the $`i`$-th particle loses the place to feedback. The use of $`h_i`$ 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 $`h_i`$ using the following smoothed kernel,

$$W_{nn}(r/h_i)=\{\begin{array}{cc}1& \mathrm{if}0r/h_i3/4,\\ \frac{\pi h_i^3}{8}W(4(r/h_i3/4))& \mathrm{if}3/4r/h_i1,\end{array}$$
(8)
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) 333Note that our smoothing definition is different from that used in Thacker et al. (2000), in other words our smoothing length is twice as large as theirs.. The number of neighbours, $`N_{\mathrm{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=(N_s/N_{\mathrm{nb},i}^n)^{1/3}`$, equation (14) of Kawata (1999) is expressed as,

$$h_i^{n+1}=h_i^n(1a+as),$$
(9)
where $`N_s=40`$ is the desired number of neighbours, $`h_i^n`$ 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`$,

$$a=\{\begin{array}{cc}0.2(1+s^2)& \mathrm{if}s<1,\\ 0.2(1+1/s^3)& \mathrm{if}s1.\end{array}$$
(10)
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.

## 3  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+.

### 3.1 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 COSMICS (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, $`\lambda `$, has been added. The initial conditions of this model are determined by the following four parameters: $`\lambda `$, $`M_{\mathrm{tot}}`$, $`\sigma _{8,\mathrm{in}}`$, and $`z_\mathrm{c}`$. The spin parameter is defined by

$$\lambda \frac{J|E|^{1/2}}{GM_{\mathrm{tot}}^{5/2}},$$
(11)
where $`J`$ is the total angular momentum of the system, $`E`$ is the total energy, and $`M_{\mathrm{tot}}`$ is the total mass of the sphere, which is composed of dark matter and gas; $`\sigma _{8,\mathrm{in}}`$ is the rms mass fluctuation in a sphere of radius $`8h^1`$ Mpc, which normalizes the amplitude of the CDM power spectrum; $`z_\mathrm{c}`$ is the expected collapse redshift. If the top-hat density perturbation has an amplitude of $`\delta _i`$ at the initial redshift, $`z_i`$, we obtain $`z_\mathrm{c}=0.36\delta _i(1+z_i)1`$ approximately (e.g. Padmanabhan, 1993). Thus, when $`z_\mathrm{c}`$ is given, $`\delta _i`$ at $`z_i`$ is determined. Kawata (1999) found that the seed galaxy which has a slow rotation corresponding to $`\lambda =0.02`$ and small-scale density fluctuations evolves into an elliptical-like system. Thus, we employ $`\lambda =0.02`$. This spin parameter is close to the minimum value possible in a CDM universe, according to the results of N-body simulations (Barnes & Efstathiou, 1987; Warren et al., 1992). We fix $`\sigma _{8,\mathrm{in}}=0.5`$ and $`z_c=3.5`$. Our simulations assume a flat universe ($`\mathrm{\Omega }=1`$) with a baryon fraction of $`\mathrm{\Omega }_\mathrm{b}=0.1`$ and a Hubble constant of $`H_0=50\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. We carry out each simulation using 9171 particles for gas and dark matter, respectively. We simulate the evolution of each model from $`z_i=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 $`\times `$ 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). Kodama & Arimoto 97 model supplies the database of SSPs with two types of IMF: ($`x`$, $`M_\mathrm{u}`$, $`M_\mathrm{l}`$) = (1.35, 60, 0.1) and ($`x`$, $`M_\mathrm{u}`$, $`M_\mathrm{l}`$) = (1.1, 60, 0.1) in the definition of Section 2.2. We adopt the data of SSPs with the IMF of ($`x`$, $`M_\mathrm{u}`$, $`M_\mathrm{l}`$) = (1.35, 60, 0.1), while we use the IMF of ($`x`$, $`M_\mathrm{u}`$, $`M_\mathrm{l}`$) = (1.35, 60, 0.2) in the numerical simulations. As mentioned in Section 4 of K01b, this inconsistency does not affect the final photometric properties.

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.  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.

### 3.2 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_\mathrm{l}`$ and $`M_\mathrm{u}`$ in Section 2.2) and the feedback parameters of $`f_v`$ 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 $`M_{\mathrm{tot}}=4\times 10^{12}`$, $`8\times 10^{11}`$, and $`2\times 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 $`M_\mathrm{u}`$ or $`M_\mathrm{l}`$ leads to larger yields, and makes the metallicity of the end-products higher. We chose $`(M_\mathrm{u},M_\mathrm{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 $`f_v`$. We fix the parameter set of $`(ϵ_{SN},f_v)=(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\times 10^{49}`$ erg, i.e. 85 % of the available energy, is found in the form of kinetic energy. This corresponds to $`f_v=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, $`f_v`$ 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 $`UV`$ and $`VK`$ colours which refer to an aperture size of 11 arcsec, and the $`V`$ band total magnitude derived from a combination of their data and the literature. Throughout this paper we adopt the distance modulus of the Coma cluster of $`mM=34.7`$ mag; the Virgo distance modulus of $`mM=31.01`$ (Graham et al., 1999) and the relative distance modulus of the Coma with respect to the Virgo of $`mM=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 ($`z0.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.  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 $`\mu _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 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 ($`\mathrm{log}\rho _{\mathrm{gas}}28`$ and $`\mathrm{log}T4`$) 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).

### 3.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. 37 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 $`UV`$ and $`M_V`$ 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.

### 3.4 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=\mathrm{}`$. We also assume that SNe Ia occur only in star particles with $`\mathrm{log}Z/Z_{\mathrm{}}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.

### 3.5 The Kormendy Relation

We now examine the physical size of the simulation end-products 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).

### 3.6  The \[Mg/Fe\]–Magnitude Relation

The abundance ratio of \[Mg/Fe\] of elliptical galaxies has been well-studied by several groups (Worthey, Faber, & Gonzalez, 1992; Jørgensen, 1999; Kuntschner, 2000; Trager et al., 2000b; Terlevich & Forbes, 2002). Strictly speaking, \[Mg/Fe\] is not “observed”, but is derived by the comparison of the theoretical model with the observed absorption line features. We compare our results with the observational results in terms of \[Mg/Fe\] rather than the absorption line features, because of convenience, although we plan to analyse the absorption line features from simulation results and compare them with the observational data more directly in future work. The observational results indicate that \[Mg/Fe\] correlates with the luminosity, i.e. the \[Mg/Fe\]–magnitude relation. Because Mg and Fe are mostly produced by SNe II 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) and the simulation results, \[E/Fe\] is measured in the $`r_{e,B}/2`$ aperture. Because Trager et al. (2000a, b) adopt \[E/Fe\] instead of \[Mg/Fe\] (see also Tantalo et al., 1998a), we derived \[E/Fe\] as well as \[Mg/Fe\] 444Note that Fig. 12 of K01b compares \[Mg/Fe\] of the simulation end-products with \[E/Fe\] of Trager et al. (2000a).. The definition of \[E/Fe\] is as follows:

$$[\mathrm{E}/\mathrm{Fe}]=([\mathrm{Z}/\mathrm{H}][\mathrm{Fe}/\mathrm{H}])/A,$$
(12)
where $`A=0.929`$ for model 4 of Trager et al. (2000a) which we refer. In model 4 of Trager et al. (2000a), \[E/Fe\] indicates the mean enhancement of “E” group, which contains C, O, Ne, Na, Mg, Si, and S, with respective to the Fe-peak elements. The quantity of \[E/Fe\] should be similar to \[Mg/Fe\], because most of the E group elements are nucleosynthetically linked. Table 3 summarises luminosity weighted \[E/Fe\] and \[Mg/Fe\] for all the model, and shows that the difference between \[E/Fe\] and \[Mg/Fe\] is less than 0.04. Therefore, we consider that it is safe to use the \[E/Fe\]-magnitude relation in order to examine the \[Mg/Fe\]-magnitude relation. The Trager et al. (2000a) data show a clear tendency that \[E/Fe\] increases with the galactic luminosity. However, the best model is incapable of reproducing this tendency and has almost constant \[E/Fe\] regardless of the luminosity. Models noIa-LM and iIIbIa-LM have similar \[E/Fe\] to the best model. As mentioned above, these \[E/Fe\] abundance ratios reflect the star formation history shown in Fig. 6. In all the models except noIaLM, the galactic wind causes the cessation of star formation, and star formation stops before the metal enrichment by SNe Ia progresses. Since \[E/Fe\] is determined only by SNe II yields, these models have similar \[E/Fe\] to that for model noIa-LM. Finally, it is revealed that the best model to reproduce the observed CMR leads to a short period of star formation and constant \[E/Fe\], irrespective of the mass of the system, which is inconsistent with the \[E/Fe\]-magnitude relation derived from observations. Because the results of \[E/Fe\] are almost identical to those of \[Mg/Fe\], we focus on \[Mg/Fe\] in the following discussion.

A more detailed look shows that the best model has a slight slope which is opposite to the \[Mg/Fe\]-magnitude relation derived from observational data, and model iIIbIa-LM (noIa-LM) has a little lower \[Mg/Fe\] than model LM (iIIbIa-LM). These trends also reflect their star formation histories. Fig. 12 shows \[Mg/Fe\] vs \[Fe/H\] of all star particles for all the models. In all the models, there is large scatter in \[Mg/Fe\] at \[Fe/H\]$`<1`$, and the mean \[Mg/Fe\] approaches 0.1 at \[Fe/H\]$`>1`$ . Model iIIbIa-LM has the least scatter, compared to the other models, because of the assumption of instantaneous recycling for SNe II. To interpret Fig. 12, the mean \[Mg/Fe\] for the total SNe II yields is useful, which is calculated by

$`[\mathrm{Mg}/\mathrm{Fe}]_{\mathrm{II}}(Z)`$ $`=`$ $`\mathrm{log}{\displaystyle \frac{_8^{60}M_{\mathrm{Mg},\mathrm{II}}(m,Z)\mathrm{\Phi }(m)𝑑m}{_8^{60}M_{\mathrm{Fe},\mathrm{II}}(m,Z)\mathrm{\Phi }(m)𝑑m}}`$ (13)
$`\mathrm{log}(Mg/Fe)_{\mathrm{}}.`$
This value equals \[Mg/Fe\] of the SNe II yields in model iIIbIa-LM, and depends on the metallicity, such as $`[\mathrm{Mg}/\mathrm{Fe}]_{\mathrm{II}}(10^4Z_{\mathrm{}})0.3`$, $`[\mathrm{Mg}/\mathrm{Fe}]_{\mathrm{II}}(10^2Z_{\mathrm{}})0.08`$, and $`[\mathrm{Mg}/\mathrm{Fe}]_{\mathrm{II}}(Z_{\mathrm{}})0.1`$. Therefore, the \[Mg/Fe\] versus \[Fe/H\] distribution in model iIIbIa-LM can be easily explained as follows. Some low metallicity stars with $`0.3>`$\[Mg/Fe\]$`>0.1`$ formed from the gas enriched by the extremely low metallicity ($`Z10^4Z_{\mathrm{}}`$) stars. The other stars (\[Mg/Fe\]$`0.1`$) are enriched mainly by stars with $`Z>10^2`$ 555The stars with \[Mg/Fe\]$`=`$0.18 and \[Fe/H\]$`<1`$ are enriched only by zero metallicity stars, because of $`[\mathrm{Mg}/\mathrm{Fe}]_{\mathrm{II}}(0)0.18`$.. The other models follow the same trend as model iIIbIa-LM, although they have much larger scatter in \[Mg/Fe\]. As seen in Fig. 3 of Gibson (1997), \[Mg/Fe\] of SNe II yields by Woosley & Weaver (1995) also depends on the progenitor mass (e.g. for solar metallicity stars, \[Mg/Fe\]$`=`$1.2 at 40 M and \[Mg/Fe\]$`=`$0.59 at 11 M). Although the difference in the life-times of SNe II progenitors is small, this time difference causes the inhomogeneous enrichment for the interstellar medium, and makes the scatter seen in Fig. 12. Some fraction of stars with \[Mg/Fe\]$`<0`$ and \[Fe/H\]$`>1`$ are enriched by SNe Ia whose yields have \[Mg/Fe\]$`=1.5`$. 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 mean \[Mg/Fe\] than models LM and iIIbIa-LM which have a shorter duration of star formation. Although inclusion of SNe Ia makes this effect negligible, this is caused by the necessarily simple recipe for the yields for intermediate mass stars, which should be improved in the near future.

### 3.7 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 $`f_v=0.002`$ as an SN feedback energy to reproduce the observed CMR, whereas K01b requires $`4\times 10^{51}`$ ergs and $`f_v=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 Woosley & Weaver (1995), whereas K01b adopts those by Notomo et al. (1997). Although Woosley & Weaver (1995) and Notomo et al. (1997) are two representative computations for SNe II nucleosynthesis, it is known that the amount of iron yield derived in Woosley & Weaver (1995) is significantly larger than that in Notomo et al. (1997), and therefore the yields of Woosley & Weaver (1995) 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 end-products 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.

## 4 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-$`\alpha `$ 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.

## Acknowledgments

We thank Michael A. Beasley for his helpful advice during the completion of this manuscript, and the referee for their thorough reading of the manuscript. We are grateful to Nobuo Arimoto and Tadayuki Kodama for kindly providing the tables of their SSPs data. We are also grateful to Edmund Bertschinger for generously providing the COSMICS program. We acknowledge the Yukawa Institute Computer Facility, the Astronomical Data Analysis Center of the National Astronomical Observatory, Japan, the Australian and Victorian Partnership for Advanced Computing, where the numerical computations for this paper were performed. This work is supported in part by the Australian Research Council through the Large Research Grant Program (A0010517), Japan Society for the Promotion of Science through the Grants-in-Aid for Scientific Research (No. 14540221), and Swinburne University through the Research Development Grants Scheme.


text:
# Low Temperature Relaxation of the Phase In an Inhomogeneous Bose Gas

## I The effective action for the difference of phases.

Below the approach to the problem of the phase difference relaxation without using the perturbation theory in the tunneling amplitude is developed. For this purpose the effective action for the difference between phases of the Bose field at the right-hand and the left-hand boundaries of the potential barrier is derived. This effective action is obtained by the integration over the bulk components of the Bose field with the fixed values of this field at the potential barrier boundaries. The obtained effective action contains the relaxation part which is proportional to the first power of the frequency of the phase field. For the small frequencies the relaxation part is found to be much larger than the usual kinetic part of the phase dynamics which is proportional to the second power of the frequency. In order to describe the phase difference relaxation process the consideration is developed in the framework of the Schwinger-Keldysh technique ,  which is very convenient for the description of quantum kinetic processes.

The generating functional of the system can be written in the form $`Z=D\psi D\overline{\psi }e\mathrm{exp}\left\{i\left(S_U+S_j\right)\right\}`$. The action $`S_U`$ of the inhomogeneous non-equilibrium Bose gas is given by

$`S_U={\displaystyle 𝑑td^3r\left\{\overline{\psi }\left[i_t+\mu \left(x\right)U\left(x\right)\right]\psi \frac{1}{2}\left(\stackrel{}{}\overline{\psi }\right)\left(\stackrel{}{}\psi \right)\frac{\lambda }{2}\left(\overline{\psi }\psi \right)^2\right\}}`$

and the term $`S_j`$ reads as $`S_j=𝑑td^3r\left(\overline{\psi }j+\overline{j}\psi \right)`$, where $`j`$ and $`\overline{j}`$ are the infinitely small sources. The integration over the time t in the action is realized over the time reversing contour and is denoted as $`𝑑t`$. The Planck constant $`\mathrm{}`$ and mass m of a Bose particle are set equal to unity $`\mathrm{}=m=1`$. The chemical potential $`\mu \left(x\right)`$ is supposed to have the constant value $`\mu _1=\mu \frac{1}{2}\mathrm{\Delta }\mu `$ for the coordinate x$`>`$d and $`\mu _2=\mu +\frac{1}{2}\mathrm{\Delta }\mu `$ for x$`<`$-d where the chemical potential $`\mu `$ obeys the equality $`\mu =\lambda n`$ and the chemical potential difference obeys the equality $`\mathrm{\Delta }\mu =\lambda \mathrm{\Delta }n`$. The small value of $`\mathrm{\Delta }n`$ results in the small value of the chemical potential difference $`\mathrm{\Delta }\mu <<\mu `$.

Due to the large magnitude of potential $`U_0`$ in the region x$`[d,d]`$ the interaction between Bose particles can be neglected and the integral for Z over the fields $`\psi ,\overline{\psi }`$ takes the Gaussian form in this region of the x-coordinate. In this connection the integral for Z over the fields $`\psi (x,\stackrel{}{r}_{};t),\overline{\psi }(x,\stackrel{}{r}_{};t)`$ for x$`[d,d]`$ with the fixed magnitudes of these fields at the boundary of the region of the nonzero external potential can be calculated. Moreover, the characteristic scale of the $`\psi `$-field variation in the time and in the $`\stackrel{}{r}_{}`$ space is supposed to be much larger than $`1/U_0`$ and $`1/\sqrt{U_0}`$, respectively. It is these fields that describe the slow relaxation process possessing the characteristic frequency proportional to the small tunneling coefficient. Due to this assumption the tunneling amplitude of the $`\psi `$ fields has the local character for the time variable t and space variables $`\stackrel{}{r}_{}`$. The magnitudes of the $`\psi `$-fields at the boundary of the region of the nonzero external potential are denoted as $`\psi (d,\stackrel{}{r}_{},t)=\psi _{s1}(\stackrel{}{r}_{},t)`$; $`\psi (d,\stackrel{}{r}_{},t)=\psi _{s2}(\stackrel{}{r}_{},t)`$. The calculation of the integral for Z over the $`\psi `$-fields in the region x$`[d,d]`$ gives the action in the form of the sum of four terms $`S=S_{vol}+S_{surf}^{\left(\rho \right)}+S_{surf}^{\left(tunn\right)}+S_j`$. The part of the action $`S_{vol}`$ is

$$S_{vol}=𝑑t\underset{U\left(x\right)=0}{}𝑑xd^2r_{}\left\{\overline{\psi }\left(i_t+\mu \right)\psi \frac{1}{2}\left(\stackrel{}{}\overline{\psi }\right)\left(\stackrel{}{}\psi \right)\frac{\lambda }{2}\left(\overline{\psi }\psi \right)^2\right\}$$
(1)
where the region of the integration over the x-coordinate is the sum of the regions $`(L,d)(d,L)`$. The parts of the action $`S_{surf}^{\left(\rho \right)}`$ and $`S_{surf}^{\left(tunn\right)}`$ can be represented in the form

$$S_{surf}^{\left(\rho \right)}=\frac{\varkappa }{2}\mathrm{coth}\left(2\varkappa d\right)𝑑td^2r_{}\left\{\psi _{s1}(\stackrel{}{r}_{},t)^2+\psi _{s2}(\stackrel{}{r}_{},t)^2\right\}$$
(2)
$$S_{surf}^{\left(tunn\right)}=\frac{\varkappa }{2\mathrm{sinh}\left(2\varkappa d\right)}𝑑td^2r_{}\left\{\overline{\psi }_{s1}(\stackrel{}{r}_{},t)\psi _{s2}(\stackrel{}{r}_{},t)e^{i\mathrm{\Delta }\mu t}+\overline{\psi }_{s2}(\stackrel{}{r}_{},t)\psi _{s1}(\stackrel{}{r}_{},t)e^{i\mathrm{\Delta }\mu t}\right\}$$
(3)
where the magnitude $`\varkappa `$ is $`\varkappa =\sqrt{2U_0}`$, the index 1 denotes the right-hand side of the bulk and the index 2 denotes the left-hand side of the bulk. The term $`S_j`$ can be written as $`S_j=𝑑t\underset{\alpha =1,2}{}\left(\overline{\psi }_{s\alpha }j+\overline{j}\psi _{s\alpha }\right)`$ where $`j`$ and $`\overline{j}`$ are the infinitely small sources. Later on we assume that the width d of the potential barrier is sufficiency large so that the inequality $`\varkappa d>>1`$ takes place and, thus, $`\mathrm{sinh}\left(2\varkappa d\right)\mathrm{cosh}\left(2\varkappa d\right)\frac{1}{2}e^{2\varkappa d}>>1`$. The term $`S_{surf}^{\left(tunn\right)}`$ describes the tunneling between the right-hand and left-hand sides of the bulk. From the form of this tunneling part of the action it can easily be seen that the tunneling amplitude for the condensate particles and for the low energy non-condensate particles has the same nonzero magnitude in contrast with the work .

Our goal is the calculation of the functional integral for Z over the $`\psi `$-fields in the bulk, i.e., in the region x$`(L,d)(d,L)`$ with the fixed values at the potential barrier boundary, and, thus, obtain the effective action for the fields at the boundary $`\psi _{s1}`$, $`\psi _{s2}`$.

Below the onedimensional character of the system and its excitations is supposed. In this supposition the considering fields depend on the x coordinate only. Later on, it is convenient to introduce the dimensionless space and time coordinates via the following change of these variables $`x\xi x;`$ $`tt/\mu `$. Thus, the system of units which we use later measures the quantities of the dimension of length in units of $`\xi `$ and the quantities of the dimension of energy in units of $`\mu `$. Moreover, we translate the x coordinate for the right-hand side of the bulk on the value -d and for the left-hand side on +d, so as the right-hand and left-hand boundaries of the potential barrier take the zero coordinates x=0.

It is convenient to represent the fields $`\psi `$ in the modulus-phase form $`\psi =Re^{i\phi }`$; $`\overline{\psi }=Re^{i\phi }`$. Below the modulus $`R`$ of the field $`\psi `$ is represented as the sum of the saddle-point configuration $`F\left(x\right)`$ and the fluctuations of the modulus field $`\rho `$ which characterize the density fluctuations

$$R(x,t)=F\left(x\right)+\rho (x,t)$$
(4)
The saddle-point configuration $`F\left(x\right)`$ obeys the equation

$$\left(1+\frac{1}{2}_x^2F^2\left(x\right)\right)F\left(x\right)=\varkappa F\left(x\right)\delta \left(x\right)$$
(5)
and has the form $`F\left(x\right)=\mathrm{tanh}\left(x+X_0\right)`$, and the constant $`X_0`$ is the small value $`X_0=1/\varkappa `$. The fluctuations $`\rho `$ can be represented as the sum of two summands

$$\rho (x,t)=\delta \rho (x,t)+\frac{_xF\left(x\right)}{_xF\left(0\right)}\rho _s\left(t\right)$$
(6)
where the fluctuations $`\delta \rho (x,t)`$ obey the zero boundary conditions $`\delta \rho (0,t)=0`$. Note, that the second summand of (6) is the zero mode of the Hamiltonian $`\widehat{H}_1=\frac{1}{2}_x^21+3F^2`$ which describe the dynamics of the fluctuations $`\rho (x,t)`$. To describe the slow process of the tunneling through the potential barrier we are interesting for the fields $`\delta \rho `$ and $`\phi `$ which are slow varying values at the space scale $`\xi `$ and the time scale $`1/\mu `$. For the derivation of the slow field effective action we substitute the field $`\psi `$ into the action (1-3) in the modulus-phase representation where the field $`R(x,t)`$ has the form (4) and where $`\rho `$ is determined by (6). At the same time, in certain terms of the action due to the slowness of $`\delta \rho `$ and $`\phi `$ the function like $`F^2\left(x\right)`$ can be replaced by the function $`1\delta \left(x\right)`$, where $`\delta \left(x\right)`$ is the $`\delta `$-function. Integrating over the fields $`\delta \rho (x,t)`$ and $`\rho _s\left(t\right)`$ in the generation functional we obtain the effective action for the slow phase field

$`S^{\left(slow\right)}\left[\phi \right]=K{\displaystyle 𝑑t\underset{U\left(x\right)=0}{}𝑑x\left\{\frac{1}{2}\left(_t\phi \right)^2\frac{1}{2}\left(_x\phi \right)^2+\left[\frac{1}{2}\left(_t\phi _s\right)^2+\gamma \mathrm{cos}\left(\mathrm{\Delta }\phi _s\left(t\right)+\mathrm{\Delta }\mu t\right)\right]\delta \left(x\right)\right\}}`$

where $`\gamma =\frac{2}{\varkappa }\mathrm{exp}\left(2\varkappa d\right)`$, the value $`\mathrm{\Delta }\phi _s\left(t\right)`$ is the difference of phases at the right-hand and the left-hand boundaries of the potential barrier $`\mathrm{\Delta }\phi _s\left(t\right)=\phi _{s2}\left(t\right)\phi _{s1}\left(t\right)`$, the constant $`K`$ is $`K=n\xi S_{}>>1`$ and $`\delta \left(x\right)`$ is the Dirac $`\delta `$-function. The value $`S_{}`$ is the square of the bulk section which is supposed to obey the inequality $`S_{}\xi ^2`$ to ensure the 1D character of the system. The dimensionless constant $`K=n\xi S_{}`$ has the large value $`K>>1`$ due to the large value of the parameter $`n\xi ^3>>1`$. The last inequality is ensured by the application of the gas approximation.

The last step to the derivation of the effective action for the phase difference $`\mathrm{\Delta }\phi _s\left(t\right)`$ is the integration in Z over the bulk components of the phase $`\phi (x,t)`$ with the fixing boundary values $`\phi _s\left(t\right)`$. For this purpose the field $`\phi (x,t)`$ is represented in the form of the Fourier expansion $`\phi (x,t)=\underset{k}{}\phi _k\left(t\right)e^{ikx}`$ with the Fourier components $`\phi _k`$ obeying the condition $`\phi _s\left(t\right)=\underset{k}{}\phi _k\left(t\right)`$ where the coordinate of the right-hand boundary of the potential barrier we choose as x=0. The integration over $`\phi _k`$ gives

$$S_{eff}\left[\phi _s\right]=𝑑t𝑑t^{}\left\{\frac{1}{2}\left(_t\phi _s\right)^2\delta \left(tt^{}\right)+\frac{1}{2}\phi _s\left(t\right)\left[<\widehat{D}>\right]_{t,t^{}}^1\phi _s\left(t^{}\right)+\gamma \mathrm{cos}\left(\mathrm{\Delta }\phi _s\left(t\right)\mathrm{\Delta }\mu t\right)\delta \left(tt^{}\right)\right\}$$
(7)
where we denote $`<\widehat{D}>=\frac{dk}{2\pi }\widehat{D}_k`$. In the frequency-momentum representation in the ”triangular” form ,  at the temperature T the phonon propagator $`\widehat{D}(\omega ,k)`$ can be written

$$\widehat{D}(\omega ,k)=\left(\begin{array}{cc}0& D^A(\omega ,k)\\ D^R(\omega ,k)& D^K(\omega ,k)\end{array}\right)$$
(8)
where $`D^{R,A}(\omega ,k)=\left[\left(\omega \pm i\delta \right)^2k^2\right]^1`$ and $`D^K(\omega ,k)=\mathrm{coth}\left(\frac{\omega }{2T}\right)\left(D^RD^A\right)`$. The simple calculations of $`<\widehat{D}>_\omega `$ give

$$<\widehat{D}>_\omega =\frac{1}{2}\left(\begin{array}{cc}0& <D^A>_\omega \\ <D^R>_\omega & <D^K>_\omega \end{array}\right)\text{ }$$
(9)
where $`<D^R>_\omega =\frac{1}{i\left(\omega +i\delta \right)}`$; $`<D^A>_\omega =\frac{1}{i\left(\omega i\delta \right)}`$; $`<D^K>_\omega =2\mathrm{coth}\left(\frac{\omega }{2T}\right)\frac{1}{i\omega }`$.

## II The effective action for small frequencies and the Fokker-Planck equation.

Introducing the half-sum and the difference of the fields instead of the fields $`\phi _{s1}`$ and $`\phi _{s2}`$ we obtain the effective action for the difference of phases at the boundary of the potential barrier. In the ”triangle” representation it can be written as

$$S\left[\phi \right]=2K𝑑t\left\{\phi \left[\frac{1}{2}\widehat{\omega }^2+i\widehat{\omega }\right]\mathrm{\Phi }+\mathrm{\Phi }\left[\frac{1}{2}\widehat{\omega }^2i\widehat{\omega }\right]\phi +i\phi \left[2\widehat{\omega }\mathrm{coth}\left(\frac{\widehat{\omega }}{2T}\right)\right]\phi V[\mathrm{\Phi },\phi ]\right\}$$
(10)
where $`\widehat{\omega }=i_t`$, the function $`V[\mathrm{\Phi },\phi ]`$ is $`V[\mathrm{\Phi },\phi ]=\gamma \mathrm{sin}\left(\mathrm{\Phi }+\mathrm{\Delta }\mu t\right)\mathrm{sin}\left(\phi \right)`$. In this representation the fields $`\mathrm{\Phi }`$ and $`\phi `$ are expressed via the fields belonging to the upper and the lower time contour branches as $`\mathrm{\Phi }=\frac{1}{2}\left(\phi _++\phi _{}\right)`$ and $`\phi =\frac{1}{2}\left(\phi _+\phi _{}\right)`$. Note that as usually in the ”triangle” representation of the Keldysh-Schwinger technique the field $`\mathrm{\Phi }`$ describes the kinetics of the phase difference and the field $`\phi `$ describes the quantum noise of this value.

The fields $`\mathrm{\Phi }`$ and $`\phi `$ can be fragmented into fast and slow components. The slow fields are the fields having the frequencies $`\omega `$ obeying the inequality $`\omega <\omega _0`$ and the fast fields are the fields with the frequencies $`\omega >\omega _0`$ where $`\gamma <<\omega _0<<1`$. The part of the action for the fast fields is obtained by the expansion of the potential $`V`$ up to the second order of the fast components $`\delta \mathrm{\Phi }`$, $`\delta \phi `$ and has the form

$$S^{\left(fast\right)}[\delta \phi ,\delta \mathrm{\Phi }]=2K𝑑t\left\{\begin{array}{c}\delta \phi \left[\frac{1}{2}\widehat{\omega }^2+i\widehat{\omega }\right]\delta \mathrm{\Phi }+\delta \mathrm{\Phi }\left[\frac{1}{2}\widehat{\omega }^2i\widehat{\omega }\right]\delta \phi +i\delta \phi \left[2\widehat{\omega }\mathrm{coth}\left(\frac{\widehat{\omega }}{2T}\right)\right]\delta \phi \\ V_{\phi \mathrm{\Phi }}^{\prime \prime }\delta \phi \delta \mathrm{\Phi }\frac{1}{2}V_{\mathrm{\Phi }\mathrm{\Phi }}^{\prime \prime }\left(\delta \mathrm{\Phi }\delta \mathrm{\Phi }+\delta \phi \delta \phi \right)\end{array}\right\}$$
(11)
$`V_{\phi \mathrm{\Phi }}^{\prime \prime }=\gamma \mathrm{cos}\left(\mathrm{\Phi }+\mathrm{\Delta }\mu t\right)\mathrm{cos}\left(\phi \right)\text{}V_{\mathrm{\Phi }\mathrm{\Phi }}^{\prime \prime }=V_{\phi \phi }^{\prime \prime }=\gamma \mathrm{sin}\left(\mathrm{\Phi }+\mathrm{\Delta }\mu t\right)\mathrm{sin}\left(\phi \right)`$

where $`\delta \phi ,\delta \mathrm{\Phi }`$ are the fast components of fields and $`\phi ,\mathrm{\Phi }`$ are the slow components. The part of the action for the slow fields can be written as

$$S^{\left(slow\right)}\left[\phi \right]=2K𝑑t\left\{2\phi _t\mathrm{\Phi }V[\mathrm{\Phi },\phi ]+i\phi \left[2\widehat{\omega }\mathrm{coth}\left(\frac{\widehat{\omega }}{2T}\right)\right]\phi \right\}$$
(12)
Due to the large value of the constant $`K>>1`$, so as the value $`\frac{1}{K}\mathrm{ln}\frac{1}{\gamma }`$ obeys the inequality $`\frac{1}{K}\mathrm{ln}\frac{1}{\gamma }<<1`$, the action $`S^{\left(fast\right)}[\delta \phi ,\delta \mathrm{\Phi }]`$ can be approximated as

$$S^{\left(fast\right)}[\delta \phi ,\delta \mathrm{\Phi }]=2K𝑑t\left\{\begin{array}{c}\delta \phi \left[\frac{1}{2}\widehat{\omega }^2+i\widehat{\omega }\right]\delta \mathrm{\Phi }+\delta \mathrm{\Phi }\left[\frac{1}{2}\widehat{\omega }^2i\widehat{\omega }\right]\delta \phi +i\delta \phi \left[2\widehat{\omega }\mathrm{coth}\left(\frac{\widehat{\omega }}{2T}\right)\right]\delta \phi \\ V_{\phi \mathrm{\Phi }}^{\prime \prime }<\delta \phi \delta \mathrm{\Phi }>_0\frac{1}{2}V_{\mathrm{\Phi }\mathrm{\Phi }}^{\prime \prime }<\delta \mathrm{\Phi }\delta \mathrm{\Phi }>_0\frac{1}{2}V_{\phi \phi }^{\prime \prime }<\delta \phi \delta \phi >_0\end{array}\right\}$$
(13)
where

$`<`$ $`\delta \phi _t\delta \phi _t>_0=0\text{}<\delta \mathrm{\Phi }_t\delta \mathrm{\Phi }_t>_0={\displaystyle \frac{i}{2K}}{\displaystyle \underset{\omega >\omega _0}{}}{\displaystyle \frac{d\omega }{2\pi }}D^{\left(K\right)}\left(\omega \right)={\displaystyle \frac{1}{4\pi K}}\mathrm{ln}\left({\displaystyle \frac{1}{\omega _0}}\right)\text{;}`$ ()
$`<`$ $`\delta \phi _t\delta \mathrm{\Phi }_t>_0=i{\displaystyle \underset{\omega >\omega _0}{}}{\displaystyle \frac{d\omega }{2\pi }}D_0^{\left(R\right)}\left(\omega \right)={\displaystyle \frac{i}{8K}}`$ (14)
The correlator $`\widehat{D}_0`$ is defined by the equality $`\widehat{D}_0=\left(\frac{\omega ^2}{2}\widehat{\sigma }_x+\frac{1}{2}<\widehat{D}>^1\right)^1`$, where $`\widehat{\sigma }_x`$ is the $`\widehat{\sigma }_x`$\- Pauli matrix. Thus, taking into account the renormalizations resulting from the fast fluctuations the effective action for the slow fields can be written as

$$S^{\left(slow\right)}[\phi ,\mathrm{\Phi }]=2K𝑑t\left\{2\phi \left(_t\mathrm{\Phi }\right)V[\mathrm{\Phi },\phi ]+i\phi \left[2\widehat{\omega }\mathrm{coth}\left(\frac{\widehat{\omega }}{2T}\right)\right]\phi \frac{i}{8K}V_{\phi \mathrm{\Phi }}^{\prime \prime }\right\}$$
(15)
In this expression the renormalizations proportional to the correlators $`<\delta \mathrm{\Phi }_t\delta \mathrm{\Phi }_t>_0`$ give the small renormalization of the constant $`\gamma `$ in the potential $`V`$. They can be neglected due to the smallness of the value $`\frac{1}{K}\mathrm{ln}\left(\frac{1}{\omega _0}\right)<<1`$. On the other side the renormalizations proportional to the correlators $`<\delta \phi _t\delta \mathrm{\Phi }_t>_0`$ are taken into account because they can become comparable with the relaxation term $`\phi \left[2\omega \mathrm{coth}\left(\frac{\omega }{2T}\right)\right]\phi `$ in the case of the low temperatures and small frequencies $`T\frac{1}{K}`$, $`\omega \frac{1}{K}`$.

It is convenient to shift the field $`\mathrm{\Phi }\mathrm{\Phi }\mathrm{\Delta }\mu t`$. As a result of this shift the action (16) takes the form

$`S^{\left(slow\right)}[\phi ,\mathrm{\Phi }]=2K{\displaystyle 𝑑t\left\{2\phi \left(_t\mathrm{\Phi }\right)V^{\left(0\right)}[\mathrm{\Phi },\phi ]+i\phi \left[2\widehat{\omega }\mathrm{coth}\left(\frac{\widehat{\omega }}{2T}\right)\right]\phi \frac{i}{8K}V_{\phi \mathrm{\Phi }}^{\left(0\right)\prime \prime }\right\}}`$

where $`V^{\left(0\right)}[\mathrm{\Phi },\phi ]=\gamma \mathrm{sin}\mathrm{\Phi }\mathrm{sin}\left(\phi \right)v_0\phi `$, at that $`v_0=2\mathrm{\Delta }\mu `$. Taking into account the smallness of the quantum noise $`\phi `$ due to the large value of the constant $`K>>1`$ the action can be expanded in the powers of the field $`\phi `$. This expansion gives

$$S_{eff}^{\left(slow\right)}[\phi ,\mathrm{\Phi }]=2K𝑑t\left\{2\phi \left(_t\mathrm{\Phi }\right)+i\phi \left[2\widehat{\omega }\mathrm{coth}\left(\frac{\widehat{\omega }}{2T}\right)+\frac{1}{16K}\gamma \mathrm{cos}\mathrm{\Phi }\right]\phi A\left(\mathrm{\Phi }\right)\phi i\frac{1}{8K}A^{}\left(\mathrm{\Phi }\right)\right\}$$
(15)
where

$$A\left(\mathrm{\Phi }\right)=\gamma \mathrm{sin}\mathrm{\Phi }v_0\text{}A^{}\left(\mathrm{\Phi }\right)=\gamma \mathrm{cos}\mathrm{\Phi }\text{ }$$
(16)
The first correction to the renormalized amplitude of the potential $`V[\mathrm{\Phi },\phi ]`$, i.e., to the constant $`\gamma `$ with respect to the value $`1/K`$ has the form $`\gamma ^{}=\gamma \left(1\frac{1}{4\pi K}\mathrm{ln}\left(\frac{K}{\gamma }\right)\right)`$. Integrating over $`\phi `$ in the generation functional Z we obtain

$`S_{eff}^{\left(slow\right)}\left[\mathrm{\Phi }\right]={\displaystyle \frac{1}{2}}Sp\mathrm{ln}\left(\widehat{\mathrm{\Gamma }}\right)+2K{\displaystyle 𝑑t𝑑t^{}\left\{\frac{i}{4}\left[2\left(_t\mathrm{\Phi }\right)A_t\right]\widehat{\mathrm{\Gamma }}_{t,t^{}}\left[2\left(_t^{}\mathrm{\Phi }\right)A_t^{}\right]i\frac{1}{8K}A_t^{}\delta \left(tt^{}\right)\right\}}`$

where the correlator $`\widehat{\mathrm{\Gamma }}_{t,t^{}}`$ is

$$\widehat{\mathrm{\Gamma }}=\left[2\widehat{\omega }\mathrm{coth}\left(\frac{\widehat{\omega }}{2T}\right)+\frac{1}{16K}\gamma ^{}\mathrm{cos}\mathrm{\Phi }\right]^1\left[2\widehat{\omega }\mathrm{coth}\left(\frac{\widehat{\omega }}{2T}\right)+\frac{1}{16K}\gamma ^{}\right]^1$$
(17)
In the approximation of the local in time action we have

$$S_{eff}^{\left(slow\right)}\left[\mathrm{\Phi }\right]=i2K𝑑t\left\{\frac{1}{4}\left[2\left(_t\mathrm{\Phi }_t\right)+A\left(\mathrm{\Phi }_t\right)\right]^2\widehat{\mathrm{\Gamma }}_0\frac{1}{8K}A^{}\left(\mathrm{\Phi }_t\right)\right\}$$
(18)
$$\widehat{\mathrm{\Gamma }}_0=\frac{1}{4T+\frac{1}{16K}\gamma ^{}}=\frac{1}{4T^{}}$$
(19)
Finally, the action can be represented in the form

$$S_{eff}^{\left(slow\right)}\left[\mathrm{\Phi }\right]=i2K𝑑t\left\{\frac{M}{2}\left(_t\mathrm{\Phi }\right)^2+\frac{M}{8}A^2\frac{1}{8K}A^{}\right\}$$
(20)
where $`M=2\widehat{\mathrm{\Gamma }}_0=\frac{1}{2T^{}}`$. Note that the generation potential Z with the action $`S_{eff}^{\left(slow\right)}\left[\mathrm{\Phi }\right]`$ (22) has the form corresponding to the quantum mechanics with the ”imaginary time” in which the role of the ”imaginary time” plays the usual real time. This action corresponds to the action which is obtained for the classical random process describing by the Langevin equation with the white noise  and having the form of the supersymmentric quantum mechanics. The Schredinger equation and the Hamiltonian corresponding to the action $`S_{eff}^{\left(slow\right)}\left[\mathrm{\Phi }\right]`$ are

$$\frac{1}{2K}_t\mathrm{\Psi }=\widehat{H}_{FP}\mathrm{\Psi }$$
(21)
$$\widehat{H}_{FP}=\frac{1}{\left(2K\right)^2}\frac{d^2}{2Md^2\mathrm{\Phi }}+\frac{M}{8}A^2\frac{1}{8K}A^{}$$
(22)
This Hamiltonian $`\widehat{H}_{FP}`$ can be represented in the factorized form

$`\widehat{H}_{FP}=\frac{1}{2M\left(2K\right)^2}\left(\frac{d}{d\mathrm{\Phi }}MKA\right)\left(\frac{d}{d\mathrm{\Phi }}+MKA\right)`$. The transformation $`\mathrm{\Psi }=\mathrm{exp}\left\{MKV\left(\mathrm{\Phi }\right)\right\}P\left(\mathrm{\Phi }\right)`$, where $`V\left(\mathrm{\Phi }\right)=\gamma ^{}\left(1\mathrm{cos}\mathrm{\Phi }\right)v_0\mathrm{\Phi }`$ and $`A=\frac{d}{d\mathrm{\Phi }}\left(V\left(\mathrm{\Phi }\right)\right)`$ transfer the Hamiltonian (18) to the form of the Fokker-Planck operator $`\widehat{L}_{FP}=\frac{1}{2M\left(2K\right)^2}\frac{d}{d\mathrm{\Phi }}\left(\frac{d}{d\mathrm{\Phi }}+2MKA\right)`$. The stationary solution of the corresponding Fokker-Planck equation in the case of the zero velocity $`v_0=0`$ has the form $`P\left(\mathrm{\Phi }\right)=C_0\mathrm{exp}\left(\frac{KV\left(\mathrm{\Phi }\right)}{T^{}}\right)`$, where $`V\left(\mathrm{\Phi }\right)=\gamma ^{}\left(1\mathrm{cos}\mathrm{\Phi }\right)`$ and the constant $`C_0`$ is $`C_0=\left(\frac{K\gamma ^{}}{2\pi T^{}}\right)^{\frac{1}{2}}`$. In the case of the nonzero velocity $`v_0`$ the stationary solution of the Fokker-Planck equation is determined by the equation

$`\left({\displaystyle \frac{d}{d\mathrm{\Phi }}}+2MKA\right)P\left(\mathrm{\Phi }\right)=j`$

and can be found taking into account the condition of the periodicity of the probability distribution $`P\left(0\right)=P\left(\mathrm{\Phi }\right)`$ and the normalization condition $`\underset{0}{\overset{2\pi }{}}𝑑\mathrm{\Phi }P\left(\mathrm{\Phi }\right)=1`$ , . The constant j is found from the normalization condition for $`P\left(\mathrm{\Phi }\right)`$ and for the case of the small $`\mathrm{\Delta }n`$ we obtain $`j=\frac{v_0}{\gamma ^{}}\mathrm{exp}\left\{2\left(\frac{K\gamma ^{}}{T^{}}\right)\right\}`$. The probability of the phase slip process is determined by the value $`\left(_t<\mathrm{\Phi }>\right)`$ $`j`$. For the zero temperature and small velocities $`v_0`$, so as $`v_0=\mathrm{\Delta }\mu <<T^{}`$, the value $`T^{}`$ can be considered as the effective temperature and for the large time t, when the action can be considered in the local in time form, determines the diffusion coefficient. In the opposite case $`v_0=\mathrm{\Delta }\mu >>T^{}`$ the diffusion coefficient is defined by the value $`\widehat{\mathrm{\Gamma }}_{t,t}=\mathrm{ln}\left(\frac{1}{v_0}\right)`$ and this value should be substituted instead of the value $`\left(\gamma ^{}/T^{}\right)`$. This can be seen from the expression (20) for $`\widehat{\mathrm{\Gamma }}`$ with taking into account that in the case of the large value of $`\mathrm{\Delta }\mu `$ compared with $`T^{}`$ the frequencies $`\omega `$ should be cut of for the small values by the chemical potential difference $`\mathrm{\Delta }\mu `$. Note, that the result for the probability of the phase slip $`\left(1/\tau _{ps}\right)v_0\mathrm{exp}\left\{2K\mathrm{ln}\left(\frac{1}{v_0}\right)\right\}`$ coincides with the result of the works \[15-17\].

The Fokker-Planck equation describing the relaxation process of the difference of phases has been obtained both for high temperatures and for low temperatures. For high temperatures $`T\gamma `$ the obtained Fokker-Planck equation is equal to the classical Langevin equation with the white noise with the correlator proportional to 1/T . For the low temperatures $`T<<\gamma `$ and the small density difference $`\mathrm{\Delta }n`$, so as the inequality $`v_0<<1/K`$ takes place, the relaxation kinetics of the difference of phases is described by the Fokker-Planck equation with the effective temperature $`T^{}`$ and the diffusion coefficient determined by both the usual temperature and the quantum fluctuations of the Bose gas. For the small temperatures and the large value of the density difference $`\mathrm{\Delta }n`$, so as $`v_0>>\gamma /K`$, the diffusion coefficient is determined by the value $`v_0`$. In the last case the result is analogues to the result of the works \[15-17\].

The author thanks S. Burmistrov, A. Kozlov and Yu. Kagan for helpful discussions. This work was supported by the Russian Foundation for Basic Research, by the Netherlands Organization for Scientific Research (NWO) and by INTAS -2001-2344.