Formation and evolution of galactic disks with a multiphase numerical model

1 - Department of Physics, Waseda University, Oh-kubo, Shinjuku-ku, Tokyo 169-8555, Japan
2 - Observatoire de Paris, DEMIRM, 61 Av. de l'Observatoire, 75014 Paris, France

Abstract
The formation and evolution of galactic disks are complex phenomena, where gas and star dynamics are coupled through star formation and the related feedback. The physical processes are so numerous and intricate that numerical models focus, in general, on one or a few of them only. We propose here a numerical model with particular attention to the multiphase nature of the interstellar medium; we consider a warm gas phase ( $\ge$ 10⁴ K), treated as a continuous fluid by an SPH algorithm, and a cold gas phase (down to 10 K), fragmented in clouds, treated by a low-dissipation sticky particle component. The two gas phases do not have the same dynamics, nor the same spatial distribution. In addition to gravity, they are coupled through mass exchanges due to heating/cooling processes, and supernovae feedback. Stars form out of the cold phase, and re-inject mass to the warm phase through SN explosions and stellar winds. The baryons are embedded in a live cold dark matter component.
Baryonic disks, initially composed of pure gas, encounter violent instabilities, and a rapid phase of star formation that slows down exponentially. Stars form in big clumps that accumulate in the center to build a bulge. Exponential metallicity gradients are obtained. External infall of gas should be included to maintain a star formation rate in the disk comparable to what is observed in present disk galaxies.

1 Introduction

The formation of a galactic disk involves such complex physical processes that only numerical computations can give some insight in their relative importance. The essential force is gravity, which is nowadays easily controlled through N-body codes, based either on grid calculations (such as FFT), or on hierarchical grouping of particles (such as the tree-code). Both algorithms optimize the computing time, growing roughly as $N \ln(N)$ . The main improvement in this domain is obtained by going towards higher and higher spatial resolution, with an ever-growing number of particles. The gravitational force is softened at small scales, to reduce unphysical two-body relaxation, resulting in a resolution equal or slightly smaller than the inter-particle distance (Romeo 1998; Dehnen 2001; Knebe et al. 2001).

It has been recognized for a long time that the dissipative component is also an essential feature in galaxy evolution, and the interstellar medium has been widely introduced in galaxy simulations, either as a continuous fluid (Eulerian grid codes, van Albada & Roberts 1981; Lagrangian SPH codes, Hernquist & Katz 1989), or through sticky particle algorithms (Combes & Gerin 1985). However, contrary to the pure gravity case, the result of simulations now depends not only on the spatial resolution and computer power, but on the physical assumptions made on the nature of the ISM and related processes.

Star formation and feedback, coupling the stars and gas by non-gravitational processes are essential (Katz 1992; Mihos & Hernquist 1994). Since the detailed processes involved in star formation are not yet well known, this introduces further uncertainties and liberties in the modeling. The most widely adopted recipe to control star formation is based on a local Schmidt law (i.e. the star formation rate is locally proportional to some power of the volume density), sometimes associated with a threshold for star formation (e.g. Friedli & Benz 1995). However, the local Schmidt law is not actually observed in galaxies, and the justification comes from an observed global empirical Schmidt law, averaged over the whole galaxy (e.g. Kennicutt 1998). Other star formation recipes are based on Jeans instability (e.g. Steinmetz & Müller 1994; Gerritsen & Icke 1997), or cloud-cloud collisions (Noguchi & Ishibashi 1986).

Due to the limited number of particles, each star "particle'' represents in fact a stellar cluster of the order of 10⁵ to $10^6~M_\odot$ . The conversion of a fraction of a gas particle into stars requires the creation of a large number of star-particles (e.g. Katz 1992), or the consideration of "hybrid'' particles, that are transiently containing both gas and stars, and have the same dynamics for a while (Mihos & Hernquist 1994). Alternatively, some "starlet'' particles are transiently created, and merged with the nearest neighbors (Jungwiert et al. 2001). Whatever the choice, some approximations are made, related to the limited resolution in time, mass and space, respectively of the order of 10⁶ yr, 10 $^6~M_\odot$ and 0.3 kpc for a typical giant spiral galaxy.

The same limitations occur when mass loss and feedback are considered. Either only stellar heating is considered (Gerritsen & Icke 1997), or the supernovae mechanical energy only (Mihos & Hernquist 1994), or more sophisticated models are used including several types of supernovae, planetary nebulae, stellar winds, evaporation and condensation (e.g. Theis et al. 1992; Samland et al. 1997; Berczik 1999).

Stellar nucleosynthesis can be followed to describe the detailed metal enrichment history of the galaxy (e.g. Lia et al. 2001). Models reveal that the frequently used instantaneous recycling is a too simple approximation (see also Jungwiert et al. 2001).

Of first importance is the thermal evolution of the gas, since thermal instabilities are responsible for cloud condensation, feeding the cold gas phase, and thus star formation (Hultman & Pharasyn 1999). Cooling and heating processes are taken into account depending on the temperature, density and metallicity of the gas (e.g. Dalgarno & McCray 1972; Sutherland & Dopita 1993). However, the time-scales involved in the thermal processes can be much smaller than the dynamical time-scales, and the spatial scales of the corresponding processes are much below the spatial resolution of the simulations. Phenomenological recipes are therefore used to convey at large-scales the resulting effects of the processes occurring below the resolution (e.g. Thacker & Couchman 2000): the heating energy is smoothed over the resolution scale, and the effective time-scales slowed down.

Many simulation models consider only one gas phase, treated as a continuous fluid at the virial temperature corresponding to a galaxy potential ( $\sim$ 10⁴ K). A simple approximation is a strictly isothermal gas, since the cooling is very efficient above 10⁴ K (Barnes & Hernquist 1991; Mihos & Hernquist 1994). Some allow gas to cool down to low temperatures ( $\sim$ 10 K), and to spread into several temperature phases (e.g. Gerritsen & Icke 1997; Yepes et al. 1997; Thacker & Couchman 2001). However, their dynamics is still that of a one-phase medium, given the insufficient density contrast that can be dealt with the numerical simulations on large-scales. Significant contrasts can be achieved only in simulations of small volumes, that approach more realistically the multiphase nature of the interstellar medium (e.g. Rosen & Bregman 1995; Wada & Norman 1999, 2001).

In the present paper, we implement most of the above processes to investigate the formation and evolution of disk galaxies (gravity, stellar and gas dynamics, star formation, feedback and massloss heating and cooling, metal enrichment...), and we focus on the large-scale multi-phase dynamics of the gas. The "warm gas'' component is modeled by a continuous fluid with an SPH code, and contains a large range of temperatures, up to hot gas re-injected by supernovae, though the bulk of the gas is found between 10⁴ and $2\times 10^4$ K. The "cold gas'' component corresponds to the cloudy and fragmented, essentially molecular medium between 10 K and 100 K. It is modeled as a separate phase, via individual cloud-particles, subject to inelastic collisions. Since at large galactic scales, the density contrast necessary to form such a fragmented structure (12 orders of magnitude) is far beyond the present computational power, the cloud-particle approximation appears well suited. For the cloudy component, the pressure forces or viscosity forces are negligible. The dissipation occurring during cloud collisions at AU scales, much below our spatial resolution, is modeled phenomenologically.

$\begin{figure} \par\includegraphics[width=9cm,clip]{ms2301f1.eps}\end{figure}$	Figure 1: Schematic presentation of the model. The four phases are represented with the main physical processes responsible for mass and energy exchanges.
Open with DEXTER

Two general classes of models have been proposed for the formation of galaxies: either big galaxies can form through the monolithic collapse of a very massive baryonic system (e.g. Eggen et al. 1962), or dwarf galaxies form first, and act as the building blocks of larger systems, that subsequently form through recursive mergers, in the hierarchical scenario (e.g. Kauffmann et al. 1999). The first scenario assumes that the gas experiences violent 3D star formation, so quickly that it has no time to settle into a disk. This has been proposed to form large elliptical galaxies or bulges. However, it is today recognized that many (if not all) elliptical galaxies can form through mergers of spiral galaxies, or smaller entities, as first suggested by Toomre & Toomre (1972), and developed further by Schweizer (1986), while bulges can form similarly by mergers, and also through secular evolution (e.g. Combes 2000). Although the most successful cosmological scheme today, based on inflation, promotes primordial density fluctuations that are scale-invariant, which together with the existence of cold dark matter, results in hierarchical formation of structures, the collapse of baryons inside dark matter halos is still an unsolved problem, with many uncertainties and free parameters. The rate at which stars are formed in a given system, and whether the stellar populations observed today can be attributed to a sudden event like the monolithic collapse, or a more progressive one like the gradual and recursive merging, is a matter of debate. It is likely that both point of views may correspond to some observed galaxy examples today, since the major merger of two gas rich giant spirals have some similarities with a monolithic event. More subtle diagnoses should be searched for, and large number statistics should be used to determine the dominant processes. In this work, we focus on the formation and evolution of a single galaxy. To achieved the required resolution we use the monolithic collapse scenario.

In Sect. 2 we detailed the assumptions of the multiphase model, and present the code and numerical methods in Sect. 3. In Sects. 4, 5 and 6 we present and discuss the results from our simulations. We give our conclusions in Sect. 7.

2 The multiphase model

We will describe in detail the different aspects of the model. Figure 1 presents a schematic view for quick reference.

2.1 Description of the four components

2.1.1 Dark matter

2.1.2 Stars

2.1.3 Warm gas

$\begin{displaymath}P=(\gamma -1) \rho u \qquad \left( \gamma = {5 \over 3} \right)\cdot \end{displaymath}$

(1)

$\begin{displaymath}\rho {{\rm d}u \over {\rm d}t}+ \tilde{P} \nabla \hbox{\vec v}= n_{\rm H}\Gamma_{\hbox{\tiny UV}}+ n_{\rm H}^2 \Lambda(u), \end{displaymath}$

(2)

2.1.4 Cold gas

A number of numerical models stop the cooling of SPH particle at 10 000 K (e.g., Weil et al. 1998; Thacker & Couchman 2000). This can only be a first approximation.

Indeed, a large amount of the gas in galaxies is observed in cold (10-100 K) dense molecular clouds. These objects are essential to the global dynamics of the galaxies since they are the location of star formation. Their physical properties and dynamics differ strongly from those of the surrounding diffuse warm gas. Indeed, each molecular cloud is in local virial equilibrium and the overall filling factor of the cold gas is very small. Consequently, it cannot be accurately described as a fluid (which supposes a rather continuous medium). Moreover, current treatment of SPH uses the ideal gas state equation to relate pressure, density and internal energy of particles which, at present days numerical resolution, have the mass of a giant molecular cloud. This ideal gas description is hardly adequate for cold, self-gravitating, GMC-like particles, which should rather obey Larson laws (Larson 1981). For these reasons, we feel that a specific non-SPH treatment should be given to these star-forming gas particles.

Some attempts have been made in this direction. Let us first mention the works by Wada & Norman (1999) and Wada (2001) who include proper cooling down to 10 K in Eulerian simulations of small regions of the galaxy ( $\sim$ 100 pc). Yepes et al. (1997), using a PPM scheme, and Hultman & Pharasyn (1999), using a SPH simulation, both identify a phase of cold, star-forming gas. However, both still treat the dynamics of the gas as a single fluid. Relying on the fact that the largest fraction of the gas in the ISM is in the form of cold clumpy fragments (molecular clouds), Noguchi (1999) adopts a sticky particle scheme to describe the collisional gas. He does not account for the warm diffuse gas. Andersen & Burkert (2000) develop a more elaborate model to treat the cold gas clouds as a collisional fluid. They take warm gas into account only through a constant external pressure exerted on cold clouds. In this work we attempt to fully model a cold and a warm gas phase with different dynamics.

While we describe the warm gas with the usual SPH technique, the cold gas particles experience gravitational forces only, and are subject to inelastic collisions following the sticky particle scheme (Levinson & Roberts 1981). Within one time step, collisions occur with a probability $P \sim \rho_{\hbox{\tiny loc}}\delta v$ , where $\rho_{\hbox{\tiny loc}}$ is the local cold gas density, and $\delta v$ is the local velocity dispersion. The only "thermal" evolution for cold gas particles is the possibility to be evaporated back to the warm gas phase through supernovae feedback (see Sect. 2.2.3)

2.2 Physical processes

2.2.1 Cooling and heating

The evolution of the thermal energy of warm gas particles follows Eq. (2). The sources of variation of the thermal energy are the following: adiabatic expansion, heating from viscosity, heating from an external radiation field (the $\Gamma_{\hbox{\tiny UV}}$ term) and radiative cooling (the $\Lambda$ term).

We use the tabulated values of $\Lambda$ in Sutherland & Dopita (1993). We have not included the effect of metallicity on the cooling: we use a constant solar metallicity in the cooling function. Indeed, most of the gas lies in the 10 000 K-20 000 K range, where metallicity has little effect. In the typical conditions of galactic dynamics, the cooling time can become much shorter than the dynamical time in resolved dense regions. Reducing the dynamical time step accordingly is not appealing since it would drastically increase de CPU requirements. Moreover, it would not be consistent with the time scale resolution of gravitational effects. Consequently we choose to damp the thermal fluctuations (the same solution is applied by Weil et al. 1998). Whenever the variation of the internal energy of a particle within one time step is larger than one fourth of its current value, it is damped to one fourth. In addition, adaptative small steps are used to follow the very fast variations of $\Lambda$ as a function of u in the 10 000 K-20 000 K range (this is similar to the integral scheme use by Thomas & Couchman 1992). This results in a stable thermal behavior for the gas.

A constant uniform heating $\Gamma_{\hbox{\tiny UV}} \sim 10^$ erg s^-1 is applied to the gas to model the background UV radiation field. The variations of $\Gamma$ only shifts the mass balance between warm and cold phase as other parameter can, such as the minimal temperature of the warm gas. It is however important to include a non-zero $\Gamma$ to help sustain the warm phase. In our study we keep $\Gamma$ constant. Gerritsen & Icke (1997) propose a detailed model for the heating.

The second source of heating comes from supernovae feedback. This is the only way in our model for cold gas to return directly to the warm gas state. It will be described in Sect. 2.2.3.

Let us explained now how the transition from warm gas to cold gas is handled. A parameter u₀ is set which defines the thermal energy at the transition. It is a free parameter, which value corresponds to a temperature of 11 000 K. A simple method would be to transform into cold gas any warm gas particle cooling down to u₀. However, owing to the different density dependence of cooling and heating terms in Eq. (2), this choice induces a sharp density threshold for turning warm into cold gas. When the average warm gas density falls noticeably below the threshold the transfer to cold gas is stopped and the cold gas phase is progressively depleted by star formation. We obtained better equilibrium between the phases with a smoother transition. Accordingly, and for the purpose of transition to cold gas only, the thermal energy of each SPH particle of warm gas is considered as a Gaussian distribution around a central value u with a dispersion of $\sigma_u$ , rather than a Dirac distribution $\delta(u)$ . This makes sense from a physical point of view since the scale resolution of the simulation is about 100 pc, leaving out thermal fluctuations below this scale, which can be crudely taken into account using the Gaussian distribution. We use a $\sigma_u$ corresponding to $\sigma_{\rm T} \sim 1 500~\rm K$ . The tail of the Gaussian distribution extending below u₀ can be considered as cold gas. Consequently we can define the fraction $f_{\rm c}$ of cold gas in a warm particle as

$\begin{displaymath}f_{\rm c}= 0.5\left(1-\hbox{erf}\left({u-u_0 \over \sigma_u}\right) \right)\cdot \end{displaymath}$

2.2.2 Star formation

There are two steps involved in modeling star formation. First, select under which conditions stars can be formed from the gas, then specify at which rate the formation occurs. Katz (1992) proposes a set of criteria: stars can form only in Jeans unstable regions and within a convergent flow of the gas. In these regions they forms at a rate following a Schmidt law. Navarro & White (1993) add the condition that the star forming regions are cooling rapidly. They reduce the cooling time criterion and the Jeans instability criterion to a density threshold. This is coherent with results derived from the observations (Martin & Kennicutt 2001). Other authors argue that the limited resolution of the simulations, where a gas particle has usually the mass of a whole GMC, makes the use of dynamical criteria unreliable (e.g Mihos & Hernquist 1994; Hultman & Pharasyn 1999). We agree with these considerations and use a simple Schmidt law to describe star formation:

$\begin{displaymath}{ {\rm d} \rho_{\hbox{\tiny star}} \over {\rm d}t} = C \rho_{\hbox{\tiny gas}}^n~. \end{displaymath}$

(3)

2.2.3 Supernovae feedback

A lower estimate for the amount of energy released in supernovae for a Salpeter IMF is 10⁴⁸ erg per solar mass of formed stars. During the first 2 Gyr of its life, a galaxy converts a large fraction of its gas content into stars. The total energy released by supernovae during this process is many times larger than the total kinetic energy of the galaxy. It appears important to take this phenomenon into account. However since most of the feedback is in the form of thermal energy in the dense gas environment of young stars, it will be quickly radiated away. A small fraction is also fed back as kinetic energy in expanding gas shells. This phenomenon is strong enough to eject gas from dwarf galaxies. For a detailed analysis of the different algorithms see Thacker & Couchman (2000). In our model, we include both thermal and kinetic feedback. Kinetic energy is imparted to any neighboring gas particles in the form of a velocity kick opposite to the forming star direction. Thermal energy is directly deposited on a small number of neighboring cold gas particles ( $\sim$ 3), which are converted into warm gas particle and brought to a temperature of 20 000 K to 100 000 K depending on the simulation. The thermal feedback is largely affected by the damping of the thermal evolution of SPH particles described in Sect. 2.2.1. It actually has an effect similar to imposing an adiabatic period of a few Myr as described in Thacker & Couchman (2000). We will show in Sect. 6.1 that the thermal feedback associated to the damping has a definite effect on the mass balance between the phases.

2.2.4 Stellar mass-loss

Jungwiert et al. (2001) show that within 10 Gyr, supernovae explosions and stellar winds can return over 40% of the stellar mass to the ISM for a star population following Scalo's IMF (Scalo 1998). Of these 40%, less than 10% are expelled in supernovae explosions. Consequently, the instantaneous recycling approximation (Tinsley 1980) used in most numerical simulations including star formation and metallicity enrichment (e.g. Hultman & Pharasin 1999) can only be a first step. Indeed, stars keep releasing gas through stellar winds long after they have separated from the gas clouds where they were born, and refill the gaseous disk proportionally to the smooth star density field, not only in locations where stars are formed. Rosen & Bregman (1995) improved on the instantaneous recycling approximation by using a constant mass loss rate for the stars. Jungwiert et al. (2001) and Lia et al. (2001) propose much more detailed models. They derive time dependent mass loss rates for various IMF, and include them in simulations. In our model we use a simple version of Jungwiert et al. mass loss rate:

$\begin{displaymath}{{\rm d} M_{\hbox{\tiny gas}} \over {\rm d}t} = M_{\hbox{\tiny star}} {c \over (t-t_{\hbox{\tiny birth}}+T_0)}\cdot \end{displaymath}$

(4)

2.2.5 Metal enrichment

The production of heavy elements in stars and subsequent ejection in SN explosions or stellar winds lead to the modification of the chemical composition and cooling properties of the ISM. Since our model includes both stellar formation and stellar mass loss, it is natural to follow the chemical evolution of the ISM. However, as mentioned in Sect. 2.2.1, we will not take the modifications of the cooling properties into account.

We use the yield y, defined as the mass of metal produced per unit mass of stellar hydrogen (Tinsley 1980). If a star was initially formed from gas with a metallicity $Z_{\rm s}$ , the gas produced by this star through the mass loss process has a metallicity:

The main difference of our model with the instantaneous recycling approximation (Tinsley 1980) is that it works within the continuous mass-loss scheme, so enriched gas is not only released in star forming regions but throughout the stellar disk.

3 Numerical methods

3.1 Tree-SPH

We use a N-body representation of the system. Gravitation and hydrodynamics are the main factors driving the evolution. To compute the corresponding forces acting on each particle, we use the now familiar tree-SPH method (Hernquist & Katz 1989).

Gravitation forces are computed using the so-called tree algorithm (Barnes & Hut 1986). This particular implementation was used in previous works, for example Semelin & Combes (2000). Following this algorithm, the contribution by a group of distant particles to the gravitational force acting a given particle can be computed using the multipole moments of the group mass distribution. It is actually computed if the angular size of the group, seen from the particle, is smaller than an opening criterion $\theta$ . Otherwise, the group is divided into subcomponents. Throughout this work we use a value of $\theta = 0.8$ . The multipole expansions are carried out to quadrupole moments. With this setup, the typical error on the gravitational force is smaller than $1\%$ (Hernquist 1987).

Hydrodynamic forces are computed through the SPH method (Gingold & Monaghan 1977; Lucy 1977). For a complete description of the method, see, for example, Monaghan (1992). However, implementations of the SPH technique come in many different flavors (see Thacker et al. 2000 for a comparative study). Let us summarize our choices for this implementation. We use the ideal gas equation of state:

To measure the quality of our implementation, we run a simulation that is becoming a standard validation test for this type of codes: the collapse of an initially static, isothermal sphere of self-gravitating gas. This was first studied by Evrard (1988) and has been reproduced by many authors since then (e.g., Hernquist & Katz 1989; Thacker et al. 2000; Springel et al. 2001; ...). Let us consider a sphere of radius R, mass M, density profile

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{ms2301f2.eps}\end{figure}$	Figure 2: Evolution of total, thermal, potential and kinetic energies of an initially isothermal sphere of gas with density profile $\rho (r) \sim {1 \over r}$ . We use dimensionless units, and G =1. The thermal and potential extrema correspond to the instant of maximal collapse.
Open with DEXTER

3.2 Implementation of mass exchanges

As presented in Sect. 2, our model takes into account four different types of particles: dark matter, stars, warm gas and cold gas. As a result of cooling, warm gas can be converted to cold gas. In dense regions, cold gas is then converted to star following the star formation rate fixed by a Schmidt law. In turn, stars expel warm gas following the mass-loss prescription (Eq. (4)). And finally, cold gas can be transformed into warm gas as a result of supernovae feedback. Except for the vaporization of cold gas to warm gas, which receive a simple treatment, we use a common numerical procedure to realize these exchanges.

A number of procedures have been used in the literature, applied mostly to star formation. In one way or another, all procedures use a variable to keep track of the fractional mass of stars already produced in a given gas particle (this can be generalized to any two phases). We can then separate procedures in two types. In the first type, pioneered by Katz (1992), the gas particle spawns small star particles each time its fractional star-mass reaches a chosen low threshold value. In this procedure the dynamics of collisionless star particles is differentiated early-on from the dynamics of the gas. The drawback is that it leads to the creation of many star particles with small masses. This increases the CPU cost and leads to two-body relaxation problems: the young stellar disk, composed of low-mass particles will be sensitive to the dynamical heating by heavier particles. We will show in Sect. 5 that this is a sensitive issue, especially with heavy dark-matter particles. Jungwiert et al. (2001) have applied this procedure to the mass-loss of the stars. They have taken care of some of the drawbacks by including schemes to merge small-mass particles of the same type.

In the second type of procedure, the entire gas particle is converted to star when the fractional star-mass reaches a high threshold, 80% for example. This was used by Mihos & Hernquist (1994). They also choose to subtract the fractional star-mass from SPH density evaluations to avoid having a 70% star particle behave as a full gas particle. Other authors choose not to, e.g. Thacker & Couchman (2000). An extension of this procedure is to use a probabilistic approach to convert the gas particle into a star particle, with a probability function depending on the fractional star-mass (see, e.g., Lia et al. 2001). The main drawback is that large amounts of star-mass may have to follow gas-particle trajectories for a long time.

$\begin{figure} \par\includegraphics[width=10cm,clip]{ms2301f3.eps}\end{figure}$	Figure 3: Reference simulation. Configuration of the three baryonic phases after 200 Myr (top) and after 1.2 Gyr (bottom).
Open with DEXTER

Our scheme, similar to Steinmetz & Muller (1994), combines the advantages of keeping the number of particles constant (no small mass particles are produced) and keeping the fractional mass of stars trapped in the gas particles very small. We simply use an SPH-like gathering procedure. When the fraction of star-mass in a gas particle is higher than 5%, we examine the star-mass content of its SPH neighbors. If there is enough star-mass in the neighboring particles to make one full star particle, the mass is gathered in the initial particle, which is then converted into a star. Using 15 SPH neighbors, the average star content of a gas particle is smaller than 10%. This procedure is also applied to convert warm gas into cold gas, and stars into warm gas. Consequently we have to keep track of a "smoothing length'' for all three types of particles.

4 Formation of a stellar disk: Reference case

4.1 Initial conditions

$\begin{displaymath}\rho_{\rm dm}(r)={ 3 M_{\rm dm} \over 4 \pi r_{\rm dm}^3 } \left( 1 + {r^2 \over r_{\rm dm}^2}\right)^\cdot \end{displaymath}$

The baryonic matter (warm gas at t=0) is distributed in a rotationally supported thin disk defined by a Miyamoto-Nagai density profile (Miyamoto & Nagai 1975):

$\begin{displaymath}\rho_{\rm bm}(r,z)= ( {z_{\rm bm}^2 M_{\rm bm} \over 4 \pi })... ...r{z})^2 \over [ r^2 + (r_{\rm bm}+\bar{z})^2]^\bar{z}^3 } \end{displaymath}$

$\begin{displaymath}\sigma_r \sim {3.36 G \Sigma \over \kappa } \qquad \sigma_{\t... ...er 2 \Omega} \qquad \sigma_z \sim \sqrt{G z_{\rm bm} \Sigma}, \end{displaymath}$

Initially, all baryonic particles are SPH gas particle with a temperature of 15 000 K, and metallicity Z=0.

For the reference simulation we use the following parameters. A Schmidt law with exponent 1.5 is used for star formation. The number of baryonic particle is 50 000, the number of dark matter particle is 10 000. This makes dark matter particles 15 times heavier than baryonic particles. The gravitational smoothing is 30 pc for baryonic matter and 300 pc for dark matter. The minimal SPH smoothing is accordingly set to 30 pc. The time step is 1 Myr and the integration is carried out over 2 Gyr.

4.2 Results

Face-on views of the various baryonic phases are presented in Fig. 3. The snapshots are taken at 0.2 Gyr and 1.2 Gyr. They show the formation of the stellar disk and the depletion of the gas phases. One interesting point in the stellar disk formation is that stellar clumps form in the outer, most unstable regions of the disk during the first few hundred Myr. They subsequently merge in the central part of the disk. We have encountered this phenomenon for a wide choice of parameters. Only with very high initial values of Toomre criterion ( $Q \sim 8$ ) throughout the disk does it disappear. It is physically unlikely that such a Q=8 disk is formed before any star formation takes place, since dissipation tends to decrease Toomre criterion in the gaseous disk. Noguchi (1999) also forms clumps for a cold collisional disk of gas and suggests that it can be a mechanism for the formation of a bulge.

Another morphological aspect worth mentioning is the behavior of the spiral structure. It is strong in the three phases during the first few 100 Myr. It then fades, especially in the stellar disk. Possible causes for this fading are the absence of external perturbations to drive the spiral structure, and the lack of resolution leading to unphysical heating of various baryonic components (see also Sect. 5)

$\begin{figure} \par\includegraphics[width=7cm,clip]{ms2301f4.eps}\end{figure}$	Figure 4: Evolution of the composition of the baryonic matter. The mass fraction of each of the three components is plotted against time.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=7cm,clip]{ms2301f5.eps} \end{figure}$	Figure 5: Star formation rate in the reference simulation, as a function of time. The star formation law is a Schmidt law with exponent 1.5.
Open with DEXTER

Figure 4 shows the time-evolution of the mass of each baryonic phase. A large part of the warm gas is almost instantaneously turned into cold gas, as thermal equilibrium is reached within a few time steps. Then mass exchanges proceed smoothly. Within 500 Myr, the warm gas mass reaches a quasi-constant value, as the balance between the gain from stellar mass loss and the loss from cooling into cold gas is reached. This is a self-regulated equilibrium: if the warm gas mass is depleted through cooling, its density falls and the cooling becomes less efficient. If it increases through stellar mass loss, the density rises, the cooling becomes more efficient, increasing the rate of transfer to the cold gas phase. On long time scales (several Gyr), the warm gas mass would eventually decrease as stars grow older and release less warm gas. The cold gas phase is continuously depleted although the process is slow after 2 Gyr. We will examine in Sect. 6 the influence of various parameters on the mass evolution of the phases. Figure 5 shows the evolution of the star formation rate. The rate is high during the first Gyr, during which 75% of the gas is turned into stars. It decreases steadily with time. This is not fully into accord with the observations, which suggest more constant star formation rate. The missing factor here is probably the accretion of gas, which would replenish the gas disk and sustain star formation rate.

$\begin{figure} \par\includegraphics[width=7cm,clip]{ms2301f6.eps}\end{figure}$	Figure 6: Surface density profiles of the various baryonic components of the disc as functions of the radius. The density scale is logarithmic. The profile are roughly exponential with 10 kpc typical radius for the warm gas, and 2 kpc for the stars.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=7cm,clip]{ms2301f7.eps}\end{figure}$	Figure 7: Rotation curve of the disc (circular velocity as a function of radius) at the beginning and the end of the reference simulation. The curve shows little evolution and is rather flat in the outer part of the disk. The dark matter contribution is large, producing a stable disk.
Open with DEXTER

The surface density profiles of each baryonic phase at 2 Gyr are plotted in Fig. 6. The stellar disk exhibits an exponential profile with a 2 kpc typical radius. The warm gas disk profile is also reasonably fitted by an exponential, with a rather larger typical radius of $\sim$ 10 kpc. The cold gas disk falls in between. Finally, Fig. 7 shows that the circular velocity profile exhibits little evolution during the first 2 Gyr of the simulation. The profile is rather flat at r > 5 kpc, with a typical velocity of $\sim$ 250 km s^-1.

4.3 Modification of the parameters

Table 1: Description of the various models used in the simulations.
Differences with model 1

Model 1 Reference model (see Sect. 4)

Model 2 Unique gravitational softening (30 pc)

Model 3 Analytical dark matter halo

Model 4 High resolution simulation: 500 000 particles cm

Model 5 Schmidt law with exponent 1.

Model 6 No stellar massloss

Model 7 No supernovae feedback

Model 8 Purely thermal supernovae feedback

Model 9 High thermal supernovae feedback

Model 10 Lighter dark matter halo

**Table 1:** Description of the various models used in the simulations.
	Differences with model 1
Model 1	Reference model (see Sect. 4)
Model 2	Unique gravitational softening (30 pc)
Model 3	Analytical dark matter halo
Model 4	High resolution simulation: 500 000 particles cm
Model 5	Schmidt law with exponent 1.
Model 6	No stellar massloss
Model 7	No supernovae feedback
Model 8	Purely thermal supernovae feedback
Model 9	High thermal supernovae feedback
Model 10	Lighter dark matter halo

$\begin{figure} \par\includegraphics[width=10cm,clip]{ms2301f8.eps}\end{figure}$	Figure 8: Edge-on views of the stellar component of the disk of different models, 1.8 Gyr after identical initial configurations. Models exhibit different degree of thickening.
Open with DEXTER

The relative complexity of the model that we propose makes it impossible to study the action of all the parameters independently. Some parameters have been tuned to produce physically acceptable results, and remain unchanged in our study. This is the case, for example, for the rate of dissipation through inelastic collisions in the cold gas: it is tuned so that cold gas loses a few percent of its velocity dispersion every 100 Myr. Another example is the temperature of the phase transition between warm and cold gas, set at 11 000 K. Changing the value would tilt the mass-equilibrium between the phases. We have tuned it in connection with other parameters such as the heating rate $\Gamma$ to produce similar masses for the two components after 2 Gyr in Model 1.

The modifications we have chosen to study define 10 models described in Fig. 8. Model 1 is the reference model described in Sect. 4. In model 2, heavy dark matter particles have the same gravitational softening as light baryonic particles. In model 3, the dark matter halo is not realized with particles, but using the static analytic profile given in Sect. 4.2. Model 4 includes 400 000 baryonic particles and 100 000 dark matter particles with a unique 15 pc gravitational softening, the time step is reduced to 0.5 Myr. In model 5, a Schmidt law with exponent 1 instead of 1.5 is used to compute the star formation rate. In model 6, the stellar mass loss process is simply shut off, and an instantaneous recycling approximation is used for metal enrichment. In model 7, supernovae feedback is switched off. In model 8, the kinetic part of the supernovae feedback is switched off. In model 9, the supernovae feedback is purely thermal again, but gas particles are heated up to 100 000 K instead of 20 000 K. Last, in model 10, the mass of the dark matter halo is three times smaller than in the reference model.

5 Disk thickness: The numerical heating

Numerical simulations of galactic disks usually focus on the baryonic matter behavior which can be easily compared to observations. It is then tempting for a given CPU cost to increase the scale resolution of baryonic matter at the expenses of the dark matter, using many light baryonic matter particles and few heavy dark matter particles. This is the case in most of the simulations in this work, where the mass ratio between dark matter and baryonic matter particles is between 5 and 15. This raises the question of how to handle two-body relaxation within each component, and more importantly between the 2 components. The two-body relaxation time-scale is usually controlled by the value of the gravitational smoothing length. There is actually an optimal value which limits both two-body relaxation and the bias in the force evaluation (see, e.g. Romeo 1997; Athanassoula 2000). How to deal with several particles species with different masses is not as well established. The main lead is to use variable softening parameters. Dehnen (2001) presents an analytic study of the effect of using variable smoothing length on the bias in the force computation. Binney & Knebe (2001) present a numerical study including two species of particles with different masses in a cosmological context, and find that using different softening values can reduce mass segregation, an effect arising from two-body relaxation between the different species.

Our early simulations showed a thickening of the stellar disk. It was quickly realized that the heavy dark matter particles were responsible for the heating of the baryon disk. Being organized in a spherical halo, they have a large velocity dispersion perpendicular to the disk which is transmitted to the baryonic particles. We have studied the importance of this effect under various simulation settings.

The thickening of the stellar disk was compared for models 1, 2, 3 and 4. Figure 9 shows the aspect of the disk seen edge on after 1.8 Gyr. Model 2, with a constant softening length of 30 pc produces the thickest disk. Model 1, with softening 30 pc for the baryonic matter and 300 pc for the dark matter, produces a thiner disk, but still quite thicker than models 3 and 4. Model 3, with an analytic halo component removes all heating of the disk due to the dark matter particles, thus providing a reference point. Interestingly model 4, including 500 000 particles also produces a very thin disk. The relaxation time has increased enough with the number of particles to produce little thickening within 2 Gyr.

Figure 10 gives a quantitative evaluation of the phenomenon. We plot the average vertical velocity dispersion of the baryonic matter within a radius of 7 kpc as a function of time for the four models described above. Initially the matter is mainly in the form of dissipative gas, which explains the initial decrease of the velocity dispersion. Conclusions pertaining to the heating by dark matter can be drawn from the latter times (>0.5 Gyr), when the disk has settled into a quiet evolution. The conclusions are the same as drawn from Fig. 9. The only new element is that, in model 4, the heating is not completely suppressed, which was not obvious in Fig. 9.

Our main conclusion is that, while different smoothing length helps reduce the transfer of vertical momentum between dark matter and baryonic matter, a greater number of particle is the best solution.

6 Sensitivity of physical quantities to parameter variations

6.1 Mass equilibrium between the phases

As explained in Sect. 4.4 and shown in Fig. 4, the parameters of the various physical processes have been chosen in the reference simulation to produce a mass equilibrium between the phases that is only slowly evolving after 2 Gyr. The gas fraction of the baryonic mass is then about 21% and slowly decreasing, which is compatible with the value of $\sim$ 10% observed in present day, $\sim$ 10 Gyr old galaxies.

$\begin{figure} \par\includegraphics[width=7cm,clip]{ms2301f9.eps}\end{figure}$	Figure 9: The vertical velocity dispersion of the stellar disc is averaged in the inner part of the disc (r < 7 kpc). Its evolution is plotted for different models, tracking the thickening of the disk.
Open with DEXTER

Most of the physical processes included in the model can affect the mass balance between the phases. We have studied the effect of varying some of the parameters involved in the processes modeling. The result are summarized in Fig. 11. The mass fraction of each baryonic component is given at 0.5 Gyr, which marks the end of the initial most unstable period., and at 2. Gyr, the end of the simulations.

Model 5 shows the effect of changing the exponent from 1.5 to 1. in the Schmidt law. The constant C in Eq. (3) is adjusted to yield similar overall star production; model 5 produces less stars in the initial period of high gas density, and more in the later period of low gas density. The final fraction of cold gas is smaller in model 5 than in the reference model, more of it having been turned into stars. The warm gas fraction is barely affected.

Model 6 is a run where stellar mass loss has been switched off. It is interesting to notice that, although stars eject warm gas particle, it is the cold phase that is mostly depleted when mass loss is switched off. We interpret this as follows. The warm gas equilibrium mass which appears in Fig. 4, is associated with a critical density fixed by the cooling process. If the gas falls below this density, transfer to the cold gas phase almost stops. The injection of mass from stellar mass loss does barely increases the density above the critical density since the cooling increases in proportion to the density. Its main effect is actually to maintain a mass-flow from stars to warm gas to cold gas. If this flow dies up at the source, by switching off stellar massloss, the warm gas settles at its critical density and stops producing cold gas. Cold gas is then steadily depleted through star formation.

Models 7, 8 and 9 study the effect of supernovae feedback. From the result of model 7, where SN feedback is absent, we assert that SN feedback does not affect the average star production, but that it has a strong influence on the balance between cold and warm gas. We can guess that the evaporation of cold gas by thermal input from SN, associated with the damping of thermal evolution is a more important factor in this regard than the kinetic feedback. This is checked in model 8, where SN feedback is purely thermal. The fraction of each baryonic component is very similar to the reference model. Model 9 shows that using a high thermal feedback does not change the balance much more.

The largest variation in this study is from 5.1% to 17.6% of warm gas at 2 Gyr between models 6 and 7. The largest variation with respect to the reference model is from 17.6% to 10.2% of warm gas at 2 Gyr. These are limited variations which point to the stability of the model as far as the composition of the baryonic matter is concerned.

Table 2: Composition of the baryonic matter at 500 Myr and 2 Gyr in different models. See Fig. 8 for properties of the models.
Model At 0.5 Gyr At 2. Gyr

Stars Warm gas Cold Gas Stars Warm gas Cold gas

1 58.7% 12.5% 28.8% 78.9% 10.9% 10.2%

5 53.7% 13.3% 33% 84.8% 10.3% 4.9%

6 65.4% 11.8% 22.8% 86.2% 8.7% 5.1%

7 58% 7.6% 34.4% 77.5% 4.9% 17.6%

8 61.2% 11.3% 27.5% 79.5% 10.4% 10.1%

9 61.4% 13.2% 25.4% 78.5% 11.2% 10.3%

**Table 2:** Composition of the baryonic matter at 500 Myr and 2 Gyr in different models. See Fig. 8 for properties of the models.
Model	At 0.5 Gyr	At 2. Gyr
	Stars	Warm gas	Cold Gas	Stars	Warm gas	Cold gas
1	58.7%	12.5%	28.8%	78.9%	10.9%	10.2%
5	53.7%	13.3%	33%	84.8%	10.3%	4.9%
6	65.4%	11.8%	22.8%	86.2%	8.7%	5.1%
7	58%	7.6%	34.4%	77.5%	4.9%	17.6%
8	61.2%	11.3%	27.5%	79.5%	10.4%	10.1%
9	61.4%	13.2%	25.4%	78.5%	11.2%	10.3%

6.2 Metallicity profiles

$\begin{figure} \par\includegraphics[width=7cm,clip]{ms2301f10.eps}\end{figure}$	Figure 10: Metallicity (log scale) as a function of radius for different models (see Fig. 8). All profiles are roughly exponential.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=7cm,clip]{ms2301f11.eps}\end{figure}$	Figure 11: Relative contribution of the dark matter halo to the disc circular velocity as a function of radius (see main body for exact definition). In model 10 (light dark matter halo), the inner part of the disc is dominated by baryonic matter.
Open with DEXTER

We establish metallicity profiles at 2 Gyr to test the reliability of the overall star formation process. Let us recall that metal enrichment is included in the simulation only as a diagnosis, without any influence on the cooling of the gas (see Sect. 2.2.1 for the justification). The metallicity profiles of the gas component for models 1, 5 and 6 are presented in Fig. 12. All three roughly fit an exponential profile. As can be expected the metallicity gradient is steeper in models 1 and 6, with a Schmidt law of exponent 1.5, than in model 5, with a Schmidt law of exponent 1. In model 6, stellar mass loss is switched off, and metals are returned to the gas using the instantaneous recycling approximation. For every unit mass of star formed, 0.4 unit mass of gas is instantaneously enriched using a yield y=0.02. The difference with the reference simulation shows at small radii, where instantaneous recycling produces a steeper metallicity gradient. Indeed this is the region where star formation is the most active, and the instantaneous recycling, which does not allow enriched gas to be carried out of star forming regions, may overestimate the process of enrichment.

6.3 Bar formation in a light dark matter halo

$\begin{figure} \par\includegraphics[width=11.5cm,clip]{ms2301f12.ps}\end{figure}$	Figure 12: Face-on plots of the four type of matter at 900 Myr for the reference model (top) and for model 10 (bottom).
Open with DEXTER

$\begin{displaymath}{\hbox{\vec v}(\rho_{\hbox{\tiny DM}}) \over \hbox{\vec v}(\rho_{\hbox{\tiny tot}})}~~ , \end{displaymath}$

The analysis of the evolution of the baryonic matter composition doesn't show any qualitative difference between models 1 and 10. Figure 14 shows the configuration of the various components at t=900 Myr, for the reference model and for model 10. The existence of a bar and a m=2 spiral mode is obvious in model 10, especially for the warm gas component. Simulation plots show that the bar is actually present during several 100 Myr.

7 Conclusions

The first motivation for this work is the recognition that the cold gas present in the ISM has morphological and dynamical properties very different from those of the warm diffuse gas. Semelin & Combes (2000) and Semelin & Combes (2002) explore some of the properties of the cold gas and its effect on the galaxy formation process. Those were local studies, where the small scale fractal structure of the gas is fully taken into account. In the present work we try, using a simple model, to include cold gas physics in a global model of galaxy formation and evolution. Ideally, what we learn from local studies of cold gas physics could be included as "sub-resolution" processes in global simulations. This interplay will hopefully develop in the future.

This work presents a multiphase numerical model for galaxy formation and evolution. The model takes into account four different types of matter obeying different dynamics: dark matter and stars, warm gas and cold gas. Dark matter and star particles obey gravity only, computed in a tree algorithm. Warm gas particles follow self-gravitating hydrodynamics implemented by a tree-SPH algorithm, and cold gas obeys gravity and undergo inelastic collisions (sticky particle scheme). Warm and cold gas are treated as two different fluids, and this is one of the innovations in the present work. In galaxies, gas and stars exchange matter and energy through a number of processes. The model takes into account the following: heating and cooling of the gas, star formation, stellar mass loss, kinetic and thermal feedback from SN explosions, metallic enrichment. We use a continuous massloss model which goes beyond the usual instantaneous recycling approximation and has not been integrated in such a multiphase model before.

Section 4 follows the evolution over 2 Gyr of an initially pure gas disk within a dark matter halo. The resulting galaxy exhibits properties similar to those of observed Milky Way like galaxies. We present in particular the evolution of the composition of baryonic matter, the star formation history, the surface density profile of the various components and the initial and final rotation curves. These quantities (except for star formation history) are reasonably similar to observed values.

The formation of a stellar disk from a pure gas component gives rise to a large number of self-gravitating clumps; those heat the disk, merge, and are driven towards the center through dynamical friction, where they contribute to form a bulge. The cold gas phase is always more concentrated than the warm gas phase. Although the continuous stellar mass loss replenishes the disk in gas along its evolution, the star formation rate decreases exponentially with a short characteristic time-scale of the order of 1 Gyr, and therefore gas infall would be required to reproduce the observed star formation history in giant galaxies like the Milky Way.

Section 5 details the problem of the thickening of the stellar disk during the 2 Gyr of the simulation. This thickening is due to the two-body relaxation between star and dark matter particles. It is shown that using different, larger, smoothing length for the heavier dark matter particles reduces the transfer of velocity dispersion, but a better result is obtained by using an analytic halo or by increasing the number of particles to 500 000. We emphasize this point since many simulations of galaxy evolution implementing complex models (including ours) are still limited to a few 10⁴ particles for systematic studies.

Section 6 explores the effect of several modifications to the model. We first study how the evolution of the baryonic matter composition is modified under several type of circumstances, such as a different star formation law, different stellar massloss, or different SN feedback. The model is rather stable in this respect. The mass composition responds to these changes, but within reasonable limits. Constraining the final baryonic matter composition from observations allows to discriminate between the various possible processes. The metallicity profiles are also computed for several models and are found to be exponential, which validates the star formation process.

If we try to estimate which of the processes and parameters studied are the most significant for the overall evolution of the galaxy, it appears that taking into consideration a proper massloss scheme for the stars comes first. Indeed it affects noticeably the composition of baryonic matter, it affects the metallicity profiles, and, although it was not studied in detail in this work, it may affect the morphology of the galaxy by changing the sites of star formation. Next, comes the specific star formation law, which is of course central to reproduce the star formation history. And third, the existence of SN feedback which allow to sustain a warm gas phase. It is an interesting point that the specific strength of the feedback is not a sensitive parameter, due to the cooling properties of the gas. Accretion of intergalactic matter, which was not studied in this work, may prove a key ingredient to reproduce the star formation history of spiral galaxies. It will be included in a future work.

Finally a simulation was run with a lighter dark matter halo, that is a less stable initial configuration of the gas disk. As expected, the disc evolved to form a bar and an m=2 spiral structure.

Our model takes a large range of physical processes into account. For some of them, such as heating of the gas by the UV background, we use a very simple algorithm. There is room for improvement at this level. But more importantly, some of the physical processes, the star formation law in particular, are still not well understood. More simulation work is needed to validate or reject the various criterion in use. Finally, there is the important issue of modeling the cold gas dynamics. We used a simple sticky particle scheme. Other choices (such as Andersen & Burkert 2000) are worth investigating in association with the SPH warm gas dynamics. Another possibility is to use SPH dynamics for the cold gas, but with an equation of state different from the ideal gas. With the correct multi-phase description of the gas, it will be possible to study a wide range of phenomena such as galaxy mergers or accretion, where cooling flows play an important role.