A&A 449, 509-518 (2006)
DOI: 10.1051/0004-6361:20042190

Modelling galaxies with a 3d multi-phase ISM

S. Harfst^1,2 - Ch. Theis^1,3 - G. Hensler³

1 - Institut für Theoretische Physik und Astrophysik, Universität Kiel, 24098 Kiel, Germany
2 - Centre for Astrophysics and Supercomputing, Swinburne University, Hawthorn, Victoria 3122, Australia
3 - Institut für Astronomie, Universität Wien, Türkenschanzstr. 17, 1180 Wien, Austria

Received 15 October 2004 / Accepted 11 November 2005

Abstract
We present a new particle code for modelling the evolution of galaxies. The code is based on a multi-phase description for the interstellar medium (ISM). We include star formation (SF), stellar feedback by massive stars and planetary nebulae, phase transitions, and interactions between gas clouds and ambient diffuse gas, namely condensation, evaporation, drag, and energy dissipation. The last is realised by radiative cooling and inelastic cloud-cloud collisions. We present new schemes for SF and stellar feedback that include a consistent calculation of the star-formation efficiency (SFE) based on ISM properties, as well as a detailed redistribution of the feedback energy into the different ISM phases.

As a first test we show a model of the evolution of a present day Milky-Way-type galaxy. Though the model exhibits a quasi-stationary behaviour in global properties like mass fractions or surface densities, the evolution of the ISM is strongly variable locally depending on the local SF and stellar feedback. We start only with two distinct phases, but a three-phase ISM is formed soon and consists of cold molecular clouds, a warm gas disk, and a hot gaseous halo. Hot gas is also found in bubbles in the disk accompanied by type II supernovae explosions. The volume-filling factor of the hot gas in the disk is $\sim$$35\%$. The mass spectrum of the clouds follows a power-law with an index of $\alpha \approx -2$. The star-formation rate (SFR) is $\sim$ $1.6~\ensuremath{{M}_{\odot}} ~\ensuremath{~{\rm yr}^{-1}}$ on average, decreasing slowly with time due to gas consumption. In order to maintain a constant SFR, gas replenishment, e.g. by infall, of the order $1~\ensuremath{{M}_{\odot}} ~\ensuremath{~{\rm yr}^{-1}}$ is required. Our model is in fair agreement with Kennicutt's (1998, ApJ, 498, 541) SF law including the cut-off at $\sim$ $10~\ensuremath{{M}_{\odot}} ~\ensuremath{{\rm pc}^{-2}}$. Models with a constant SFE, i.e. no feedback on the SF, fail to reproduce Kennicutt's law.

We performed a parameter study varying the particle resolution, feedback energy, cloud radius, SF time scale, and metallicity. In most these cases the evolution of the model galaxy was not significantly different to our reference model. Increasing the feedback energy by a factor of 4-5 lowers the SF rate by $\sim$ $0.5~\ensuremath{{M}_{\odot}} ~\ensuremath{~{\rm yr}^{-1}}$, while decreasing the metallicity by a factor of $\sim$100 increases the mass fraction of the hot gas from about 10% to 30%.

Key words: methods: N-body simulations - ISM: evolution - Galaxy: evolution - galaxies: ISM - galaxies: kinematics and dynamics

1 Introduction

The treatment of the interstellar medium (ISM) is a crucial part in modelling galaxies due to the importance of dissipation, star formation (SF), and feedback. Single-phase TREE-SPH models suffer from over-cooling (Thacker et al. 2000; Pearce et al. 1999) because the multi-phase structure of the ISM is not resolved. As a result, baryonic matter loses angular momentum to the dark matter (DM) halos, and the disks formed in these simulations are too concentrated. Also, SF is overly efficient. On the other hand, a multi-phase ISM is a pre-requisite for the classical grid-based chemo-dynamical (CD) models (Burkert et al. 1992; Samland et al. 1997; Samland & Gerhard 2003; Theis et al. 1992), but they lack the geometrical flexibility of a (3d-)particle code.

An approach to including a multi-phase treatment of the ISM in a particle code has been presented by Semelin & Combes (2002), who combined smoothed particles hydrodynamics (SPH) with a sticky particle method (Noguchi 1988; Carlberg 1984; Combes & Gerin 1985; Palous et al. 1993; Theis & Hensler 1993). In their multi-phase ISM code, sticky particles are used to describe a cold, dense cloudy medium, while the diffuse ISM is modelled by SPH. Stars form from cold gas, and feedback processes include the heating of cold gas by supernovae and stellar mass loss (including metal enrichment). Cooling leads to a phase transition from warm to cold gas. The main advantage of this approach is the dynamical treatment of the cloud component: clouds are assumed to move on ballistic orbits and to dissipate kinetic energy in collisions.

The main purpose of this paper is to present a new extended TREE-SPH code with a multi-phase description of the ISM adopted from CD models. As in Semelin & Combes (2002), this is achieved by a combination of a sticky particle scheme and SPH. However, we use a different sticky particle scheme based on an individual treatment of single clouds. Advantages of this scheme are the build-up of a cloud mass spectrum, the individual handling of cloud collisions, and the localised treatment of phase transitions between clouds and the ambient diffuse gas. Furthermore, we apply a refined SF recipe suggested by Elmegreen & Efremov (1997, hereafter EE97). This scheme considers star formation within individual clouds including local feedback processes. As a result the star-formation efficiency is locally variable depending on the cloud mass and ambient ISM pressure.

A detailed description of the model including all the processes is given in Sect. 2. Details on the numerical treatment can be found in Sect. 3. The code is then applied to an isolated Milky-Way-type galaxy (Sect. 4), and the results (including a short parameter study) are discussed in Sect. 5. Conclusions and future applications of this code are finally presented in Sect. 6.

   
2 The model

	Figure 1: The network of processes connecting the different gaseous and stellar components. Clouds are the sites of SF. Stars return mass and inject energy into the multi-phase ISM. The two gas phases exchange mass and momentum by means of condensation or evaporation (C/E) and drag. The energy dissipation is due to radiative cooling and cloud-cloud collisions.
Open with DEXTER

We use particles to describe all galactic components. Each particle is assigned a special type representing either stars, clouds, (diffuse) gas, or DM. The different components are connected by a network of processes including SF, stellar feedback, cooling, cloud-cloud collisions, condensation, evaporation, and drag, as shown schematically in Fig. 1.

   
2.1 Stars

Each star particle represents a cluster of stars formed at the same time. The initial stellar mass distribution is given by the multi-power law IMF of Kroupa et al. (1993) with the lower and upper mass limits of  $0.1~\ensuremath{{M}_{\odot}}$ and  $100~\ensuremath{{M}_{\odot}}$, respectively.

The stellar life times (depending on stellar mass and metallicity) are approximated according to Raiteri et al. (1996) to give reasonable fits to the Padova evolutionary tracks (Bertelli et al. 1994; Alongi et al. 1993; Bressan et al. 1993) in the mass range from $0.6~\ensuremath{{M}_{\odot}}$ to $120~\ensuremath{{M}_{\odot}}$ and for metallicities Z = 7 $\times$ 10^-5-0.03.

We assume that stars with masses above $8~\ensuremath{{M}_{\odot}}$ end their lives as type II supernovae (SNII), and stars from $1{-}8~\ensuremath{{M}_{\odot}}$ end as planetary nebulae (PNe). The mass fractions of these stars are $9\%$ and $34\%$, respectively. The life times of the high mass stars range from  $3~\ensuremath{{\rm Myr}}$to  $39~\ensuremath{{\rm Myr}}$. The stellar feedback is described below.

   
2.2 Clouds

The cloud particles describe a cold, clumpy phase of the ISM, i.e. the molecular clouds, which are also the sites of SF. The radius  $\ensuremath{r_{\rm cld}}$ of an individual cloud is given by the following relation based on observations (Rivolo & Solomon 1988):

$$% r_{\rm cld} = \ensuremath{A_{\rm cld}}\cdot \sqrt{\ensurema... ...}\over 10^6~\ensuremath{{M}_{\odot}} }~\ensuremath{{\rm pc}} ,$$ (1)

where the constant $\ensuremath{A_{\rm cld}}$ is set to 50.

The giant molecular clouds represented by cloud particles consist of gas at different densities and temperatures. Since we do not resolve this structure, we do not define a temperature or pressure for individual clouds. However, we assume that the clouds are influenced by the pressure of the surrounding ISM (see Sect. 2.6). Radiative cooling within clouds is not taken into account. However, we allow for dissipation of kinetic energy in cloud-cloud collisions using the dissipation scheme of Theis & Hensler (1993). Two clouds undergo a complete inelastic collision if they physically collide within the next time step. Additionally, their orbital angular momentum must also be less than the break-up spin of the final cloud, and the cloud mass must not exceed a critical mass limit (here: $M_{\rm coll,lim} = \ensuremath{1 \times 10^{7}} ~\ensuremath{{M}_{\odot}}$). In case of a sticky collision, the two cloud particles are replaced by one particle keeping the centre-of-mass data of the former two clouds.

2.3 Diffuse gas

In a first approximation, the diffuse gas can be described as continuous self-gravitating fluid. We model this component using an SPH approach. The gas is assumed to obey the ideal equation of state with adiabatic index $\gamma=5/3$:

$$% \ensuremath{P_{\rm gas}} =(\gamma-1)\rho_{\rm gas} u,$$ (2)

where $\ensuremath{P_{\rm gas}}$ is the pressure, $\rho$ the density, and u the energy per unit mass of the gas. The thermal energy evolution is governed by cooling and heating processes. For the latter we only take the most dominant source, the SNII, into account. Dissipation is performed by the metal-dependent cooling function from Böhringer & Hensler (1989).

The diffuse gas can have temperatures from a few $10^3~\ensuremath{{\rm K}}$ up to several $10^6~\ensuremath{{\rm K}}$. In the following, we will also distinguish between warm and hot gas, meaning gas at temperatures below or above $10^4~\ensuremath{{\rm K}}$, respectively.

2.4 Dark matter

In our code DM can be modelled either with particles (live halo) or as a static potential. In this work we focus on the late evolution of an isolated galaxy where a static DM potential is a reasonable approximation. We use the potential from Kuijken & Dubinski (1995), which is a flattened, isothermal potential with a core and a finite extent.

2.5 Interactions between the ISM phases

The ISM is treated as two dynamically independent gas phases. This approach has been successfully applied in chemo-dynamical models (Samland et al. 1997; Theis et al. 1992). The two gas phases are linked by the processes of condensation/evaporation and a drag force. Following Cowie et al. (1981), the mass exchange rate by condensation and evaporation for an individual cloud particle is calculated from

$$% \begin{array}{ll} \dot{m}_{\rm cld} =& \left({\displaystyl... ...ht.\ensuremath{{\rm g}}\ensuremath{{\rm s}^{-1}} , \end{array}$$ (3)

where the parameter $\sigma_{0}$ is defined as

$$% \sigma_{0} = \left({\ensuremath{T_{\rm gas}}\over \ensurema... ..._{\rm gas}}\over \ensuremath{{\rm cm}^{-3}} }\right)^{-1}\cdot$$ (4)

Here, $\ensuremath{T_{\rm gas}}$ and $\ensuremath{n_{\rm gas}}$ are the local temperature and number density of the diffuse gas. The clouds are subject to a drag force  ${\vec F}_{\rm D}$ due to ram pressure, which is given by (Shu et al. 1972)

$$% {\vec F}_{\rm D} = -C_{\rm D}~ \pi r_{\rm cld}^2~ \rho_{\rm gas}~ v_{\rm rel}~ {\vec v}_{\rm rel},$$ (5)

where $\rho_{\rm gas}$ is the local gas density and  $v_{\rm rel}$the velocity of the cloud relative to the gas. The constant $C_{\rm D}$ is not well-determined, and it can be in the range 0.1-1. We used $C_{\rm D} = 0.1$.

   
2.6 Star formation

The treatment of SF is usually based on the Schmidt Law (Schmidt 1959)

$$% {\it SFR}/{\rm area} = c_n \cdot \Sigma_g^n,$$ (6)

with the gas surface density $\Sigma_{\rm g}$ and $n \approx 1.5$. The parameter c_n combines SF time scale and star-formation efficiency (SFE). A reasonable estimate for the time scale is given by the dynamical time. The SFE, on the other hand, cannot be determined from current theories of SF and, therefore, remains a free parameter. A constant SFE also cannot account for local feedback processes, thereby prohibiting any self-regulation of the SF. However, Köppen et al. (1995) have demonstrated the importance of negative stellar feedback, which leads to a Schmidt-type, strongly self-regulated SF.

Though the standard approach for SF works well for isolated galaxies, it has difficulties reproducing observations of interacting galaxies (Mihos et al. 1993,1992). Barnes (2004) suggested a second shock-induced SF mode, but this introduces another free parameter in the SF recipe.

	Figure 2: The SF scheme: (t0) inactive cloud (no SF); (t1) an embedded star cluster formed with a local SF efficiency; (t2) the cloud fragmented by SNII energy input.
Open with DEXTER

Different to such an approach, we apply a single SF law, but with a temporally and spatially variable SFE using a model suggested by EE97. They assign each molecular cloud an individual SFE controlled by the local ISM properties (mass of the molecular cloud, pressure of ambient diffuse gas) and the stellar feedback. According to EE97, the star formation is most efficient for massive clouds embedded in a high pressure ambient medium. A pre-requisite of using the EE97 SFE scheme is a multi-phase treatment of the ISM as proposed in this paper.

The main concept of our SF implementation is that stars form in clouds and that clouds are destroyed by stellar feedback, thereby self-regulating the SF (Fig. 2). Each cloud is assumed to form stars only after a time of inactivity  $\ensuremath{\tau_{\rm ia}}$ (t₀ in Fig. 2), because not all gas in clouds is dense or unstable enough for immediate SF. Once the SF process has started, an embedded star cluster is formed (t₁) instantaneously.

Following EE97 we use a variable SF efficiency $\epsilon$depending on local properties of the ISM. A typical SFE (cloud mass of $10^4{-}10^6~\ensuremath{{M}_{\odot}}$ and pressure typical for solar vicinity) is between 5% and 10% (see Fig. 4 in EE97). For details of the implementation see Appendix A.

The only free parameter in our SF recipe is $\ensuremath{\tau_{\rm ia}}$. It can be interpreted as a global SF time scale: assuming an SFE of $\epsilon \sim 0.1$, a total mass  $\ensuremath{M_{\rm cld}}$ in clouds of a few $10^9~\ensuremath{{M}_{\odot}}$ and  $\ensuremath{\tau_{\rm ia}}$ of a few $100~\ensuremath{{\rm Myr}}$, we get a global ${\it SFR} \sim \epsilon \ensuremath{M_{\rm cld}}\ensuremath{\tau_{\rm ia}} ^{-1}$ of $\sim$ $1~\ensuremath{{M}_{\odot}} ~\ensuremath{~{\rm yr}^{-1}}$.

2.7 Stellar feedback

Once the SFE has been determined, the energy released by SNII in each time step can be computed from the chosen IMF and stellar life times (Sect. 2.1). We assume $\ensuremath{\epsilon_{\rm SN}} =\ensuremath{2 \times 10^{50}} ~\ensuremath{{\rm erg}}$ for the energy released per SNII event.

This energy injection disrupts the cloud into smaller fragments (t₂ in Fig. 2) leaving a new stellar particle. The time for disruption $t_{\rm dis}$ ( =t₂ - t₁) is determined from the energy input: in our model the energy input of the SNII drives an expanding shell into the cloud, where the cloud mass accumulates. The expansion of the SN shell can be calculated from a self-similar solution (Brown et al. 1995; Theis et al. 1998). Assuming a spherical cloud with a 1/r-density profile in agreement with the mass-radius relation (Eq. (1)), we get

$$% r_{\rm sh}(t) = 0.918\left( {\dot{E} \over \rho_1} \right)^{0.25} \cdot t^{0.75},$$ (7)

$$% v_{\rm sh}(t) = 0.689\left( {\dot{E} \over \rho_1} \right)^{0.25} \cdot t^{-0.25},$$ (8)

where $\dot{E}$ is the energy injection rate and $\rho_1 = \ensuremath{m_{\rm cld}} /\ensuremath{r_{\rm cld}} ^2$.

It is assumed that the cloud is disrupted when the shell radius reaches the cloud radius. The fragments are then placed symmetrically on the shell with velocities corresponding to the expansion velocity at  $t_{\rm dis}$. The energy not consumed by cloud disruption, i.e. the energy of SNII exploding after  $t_{\rm dis}$, is returned as thermal feedback to the hot gas phase. Depending on cloud mass and SF efficiency, the fragmentation uses about 0.1-5% of the total SNII energy, of which $\sim$$20\%$ ends up as the kinetic energy of the fragments.

The mass returned by SNII is also added to the surrounding gas. It is determined from the IMF and stellar life times assuming that the mass ejected per SNII event $M_{\rm SNII,ej}$ is given by (Raiteri et al. 1996)

$$% M_{\rm SNII,ej}(m) = 0.77m^{1.06}~\ensuremath{{M}_{\odot}} ,$$ (9)

where m is the mass of a star. Along with the SNII feedback, we also allow for feedback by PNe. After each time step, the mass  $M_{\rm PN}$ returned to the ISM by PNe is calculated by assuming that the PN mass ejected by a star of mass m is given by (Weidemann 2000)

$$% M_{\rm PN,ej}(m) = 0.91m-0.45~\ensuremath{{M}_{\odot}} .$$ (10)

Gradual metal enrichment due to stellar mass release and type I SNe are not included so far. This is justified for the kind of models presented here (Sect. 4.1) as metal enrichment has little effect, as we show later (Sect. 5.3).

   
3 Numerical treatment

3.1 Stellar dynamics

In our code the gravitational force on a particle can be calculated as a combination of the self-gravity of the N-body system and an external force. External forces - if applied - usually mimic a static dark-matter halo. Self-gravity is calculated using a tree scheme. We have included two different tree schemes: the widely used one of Barnes & Hut (1986) and a new tree algorithm proposed by Dehnen (2002). An advantage of the Dehnen-TREE is that it conserves momentum by construction. Additionally, an error control is included and the CPU scaling is basically linear with the number of particles. In our implementation, the Dehnen-TREE is more than an order of magnitude faster than the Barnes&Hut-TREE. Hence, we mainly apply the Dehnen-TREE here. For easy comparison with other N-body solvers, we realise gravitational softening by a Plummer potential with a softening length $\ensuremath{{\epsilon}_{\rm s}} = 0.1~\ensuremath{{\rm kpc}}$.

3.2 Gas dynamics

We used the SPH formalism to compute the hydrodynamics of the diffuse gas and a sticky particle scheme for the clouds. Details of our SPH implementation are given in Appendix B. The dynamics of the clouds are described by the sticky particle scheme of Theis & Hensler (1993), according to whom clouds move on ballistic orbits unless they undergo inelastic collisions. The treatment of cloud-cloud collisions is described in Sect. 2.2.

3.3 Time integration

The time integration of the system is done with a leap-frog scheme. Individual time steps are used for each particle. The time step for an SPH particle p is chosen to satisfy a modified form of the Courant condition (for details see Appendix B). The time step of all other particles is based on dynamical criteria calculated by

$$% \Delta t_p \le \eta_{\rm ts} \cdot \min\left( {\ensuremath{... ..._p}, \sqrt{\ensuremath{{\epsilon}_{\rm s}}\over a_p } \right),$$ (11)

where a_p is the acceleration and $\ensuremath{{\epsilon}_{\rm s}}$ the softening length of particle p. The constant  $\eta_{\rm ts}$ is of the order of unity (we set it to 0.4).

The largest individual time step $\Delta t_p$ also defines the system time step  $\Delta t_{\rm sys}$, which is further limited not to exceed $\sim$ $1~\ensuremath{{\rm Myr}}$. Individual particle time steps are chosen to be a power-of-two subdivision of  $\Delta t_{\rm sys}$ (Aarseth 1985).

3.4 Mass exchange and new particles

Several processes included in our model lead to an exchange of mass between particles or the creation or removal of particles. These processes are introduced in such a way that matter and momentum are conserved.

The mass exchange by condensation and evaporation is calculated after each system time step  $\ensuremath{\Delta t_{\rm sys}}$. The transferred mass  $\dot{m}_{\rm cld}\ensuremath{\Delta t_{\rm sys}}$ is added to the cloud. The same mass is subtracted from the surrounding SPH particles taking the weighting from the kernel function into account. In addition, the velocity of particles gaining mass (i.e. clouds in the case of condensation and SPH particles for evaporation) is changed in order to conserve momentum. The mass exchange from stars to SPH particles (SNII) and clouds (PNe) is done in a similar way.

Whenever an SF process occurs, the star-forming cloud is removed from the simulation. Instead a new star particle of the mass  $\epsilon\ensuremath{m_{\rm cld}}$ is added at the same position with the same velocity. Furthermore, the clouds formed in expanding shells are created by adding four new cloud particles carrying the remaining mass. Their positions and velocities are chosen according to Eqs. (7) and (8) thereby conserving mass and momentum.

3.5 Mass limits for particles

Our model allows particles to be created, merged and changed in mass. A few constraints are used for the sake of numerical accuracy and particle numbers. In order to prevent the formation of too many low mass star particles, a lower mass limit $\ensuremath{M_{\rm SF,lim}} = \ensuremath{2.5 \times 10^{4}} ~\ensuremath{{M}_{\odot}}$ and a lower pressure limit $\ensuremath{P_{\rm SF,lim}} = 0.1~\ensuremath{P_{\odot}}$ for SF are introduced. The maximum mass for clouds is given by  $M_{\rm coll,lim}$(Sect. 2.2) but a lower mass limit for clouds  $M_{\rm cld,lim}$ is also set to

$$% M_{\rm cld,lim} = 0.1 \ensuremath{M_{\rm SF,lim}} .$$ (12)

If the mass of a cloud becomes smaller than  $M_{\rm cld,lim}$, the particle is removed from the system and the mass is distributed evenly among neighbouring cloud particles within  $300~\ensuremath{{\rm pc}}$.

The SPH particles can change their mass due to condensation and evaporation. An upper mass limit is introduced to avoid SPH particles getting too massive. Whenever an SPH particle reaches this limit, it is split into four new particles. The new equal mass particles have the same center of mass, velocity, angular momentum, and specific energy as the particle they replace. A lower mass limit is also used where the mass of an SPH particle is distributed to the neighbouring particles. The upper and lower mass limits are set to  $4~m_{\rm SPH}$ and  $0.1~m_{\rm SPH}$, respectively, where  $m_{\rm SPH}$ is the SPH particle mass at the beginning of the simulation.

3.6 Tests

We have tested our code thoroughly. Energy and angular momentum conservation is better than 1% for a pure N-body simulation following the evolution of a disk galaxy over many rotational periods. We tested our SPH implementation with the standard Evrard-test (Evrard 1988) and find that energy conservation is within 3%. In a full model, angular momentum is conserved within 3%.

   
4 Disk galaxy model

   
4.1 Initial conditions

A galaxy that is similar to the present-day Milky Way was set up using the galaxy-building package described by Kuijken & Dubinski (1995, hereafter KD95). The models of KD95 consist of three components: the baryonic disk and bulge and a DM halo. A distribution function for each component yields a unique density for a given potential. By solving the Boltzmann equation for the combined system, the galaxy can be set up in a self-consistent way. By construction these models are in virial equilibrium, and KD95 have also demonstrated that their models are stable over many rotation periods.

We choose model A of KD95 as a reference model (Table 1). The scale length of the disk is $R_{\rm d} = 4~\ensuremath{{\rm kpc}}$, and the total mass of the galaxy is 2.47 $\times$ $10^{11}~\ensuremath{{M}_{\odot}}$. The mass and radius of the disk are set to 0.34 $\times$ $10^{11}~\ensuremath{{M}_{\odot}}$ and $18.0~\ensuremath{{\rm kpc}}$. The corresponding parameters for the bulge are $\ensuremath{0.17 \times 10^{11}} ~\ensuremath{{M}_{\odot}}$ and $4.0~\ensuremath{{\rm kpc}}$ and for the DM halo $\ensuremath{1.94 \times 10^{11}} ~\ensuremath{{M}_{\odot}}$ and $87.9~\ensuremath{{\rm kpc}}$, respectively. The particle numbers for disk and bulge, $N_{\rm d} = 20~000$ and $N_{\rm b} = 10~000$, are chosen to have particles of nearly equal mass. The DM halo is added as a static potential. The disk is initially in rotational equilibrium with a small fraction of disordered motion superimposed. This results in a Toomre parameter $Q\equiv \kappa\sigma_R/(3.36G\Sigma)$ exceeding 1.9 all over the disk. Thus, according to Polyachenko et al. (1997), the disk is stable against all axisymmetric and non-axisymmetric perturbations.

Table 1: Properties of the initial disk galaxy model.

A cloudy gas phase is added by treating every fifth particle from the stellar disk as a cloud. The total mass in clouds is $\ensuremath{M_{\rm cld}}\approx 6.9$ $\times$ $10^9~{M}_{\odot}$, and each cloud is randomly assigned a time of inactivity  $\ensuremath{\tau_{\rm ia}}$ between 0 and  $200~\ensuremath{{\rm Myr}}$. Initially, all clouds have the same mass ($\sim$ $\ensuremath{1.7 \times 10^{6}} ~\ensuremath{{M}_{\odot}}$), but a mass spectrum is built up within the first Gyr of evolution.

Finally, a diffuse gas component with a total mass of ${M}_{\rm gas} = 1$ $\times$ $10^9~{M}_{\odot}$ is added. This is more difficult because no simple equilibrium state exists for this gas phase: it is not only affected by the gravitational potential but also by cooling and heating processes. The latter varies locally and temporally so that the diffuse gas can only reach a dynamical equilibrium, as is also revealed in high-resolution simulations of the ISM (de Avillez & Breitschwerdt 2004). Testing different initial setups, it turned out to be most efficient to start far from equilibrium (for the diffuse component only). Then, the diffuse gas phase relaxes on a few dynamical time scales to a quasi-equilibrium state. In detail, we started with a slowly rotating homogeneous gaseous halo with a radius of  $25~\ensuremath{{\rm kpc}}$ and a temperature of $\ensuremath{5.4 \times 10^{4}} ~\ensuremath{{\rm K}}$, which relaxed to a quasi-equilibrium within a few 100 Myr.

Since we do not follow the chemical evolution of the galaxy, we apply solar abundances where a metallicity is needed. In order to treat the feedback of "old'' star particles, each star particle is assigned an age. For the bulge stars the age is $10~\ensuremath{{\rm Gyr}}$ and for disk stars $0{-}10~\ensuremath{{\rm Gyr}}$.

   
4.2 Evolution of the reference model

Figure 3: The evolution of stellar surface density (at $1~\ensuremath{{\rm Gyr}}$ and $3~\ensuremath{{\rm Gyr}}$). The stellar disk becomes slightly thicker with time but is otherwise stable with weak transient spiral patterns. The surface density plots (Figs. 3 and 4) were computed using the NEMO Stellar Dynamics Toolbox (Teuben 1995). The particles were binned in a grid with a cell size of $40~\ensuremath{{\rm pc}}$. The resulting map was then smoothed with a Gaussian kernel with ${\it FWHM} = 1~\ensuremath{{\rm kpc}}$.

Open with DEXTER

	Figure 4: The evolution of cloud surface density (at  $1~\ensuremath{{\rm Gyr}}$ and  $3~\ensuremath{{\rm Gyr}}$). The cloud disk is stable, but the overall surface density slowly decreases. For details on how the surface density was computed see Fig. 3.
Open with DEXTER

Starting from these initial conditions, the galaxy model is evolved for  $3~\ensuremath{{\rm Gyr}}$ with all processes switched on. After an initial phase ( $500~\ensuremath{{\rm Myr}}$) dominated by the collapse of the gaseous halo, the galaxy reaches a quasi-equilibrium defined by the interplay of the gravitational potential, radiative cooling, heating by SNII, and condensation and evaporation. When the equilibrium is reached, the diffuse gas phase has formed a warm disk and a hot halo component. Their mass ratio of about 7 remains nearly constant for the rest of the simulation. The equilibrium state reached is almost independent of the intial setup of the gas halo; i.e. the following evolution is not influenced by our choice.

Snapshots of the evolution of the Milky Way model can be seen in Figs. 3 to 6. The stellar surface density plots in Fig. 3 demonstrate the stability of the disk during the $3~\ensuremath{{\rm Gyr}}$ evolution. Profiles of the surface density show that the radial mass distribution does not change except for the central $2~\ensuremath{{\rm kpc}}$ where the surface density is increased by up to a factor of 5, caused by the dissipation-induced infall of clouds that later form stars. Face-on views show weak transient spiral patterns, e.g. at $t = 2~\ensuremath{{\rm Gyr}}$. Edge-on views reveal a thickening of the disk. Throughout the disk the 95%-mass height  $z_{\rm 95,str}$increases by 47% from  $0.75~\ensuremath{{\rm kpc}}$ to about  $1.1~\ensuremath{{\rm kpc}}$. Only in the central  $2~\ensuremath{{\rm kpc}}$ can no thickening be seen.

Figure 5: The of evolution of gas density (at $1~\ensuremath{{\rm Gyr}}$ and $3~\ensuremath{{\rm Gyr}}$). In each panel the xy- ( bottom), xz-projection ( top), and density scale are shown. Starting from a homogenous distribution, the diffuse gas collapses and forms a thin gas disk within $1~\ensuremath{{\rm Gyr}}$. While the mass of the gaseous disk is nearly constant, its radius is growing slowly. Densities in the disk are of the order of  $1~\ensuremath{{\rm cm}^{-3}}$. The halo has densities from  $10^{-4}~\ensuremath{{\rm cm}^{-3}}$ to  $10^{-6}~\ensuremath{{\rm cm}^{-3}}$. Densities were computed on a 100 $\times$ 100-grid using the SPH formalism.

Open with DEXTER

Figure 6: The evolution of gas temperature (at  $1~\ensuremath{{\rm Gyr}}$ and $3~\ensuremath{{\rm Gyr}}$). In each panel the xy- ( bottom), xz-projection ( top), and temperature scale are shown. After the initial collapse, the diffuse gas forms a disk at about  $10^4~\ensuremath{{\rm K}}$. The central part of the disk, as well as small bubbles in the outer parts, are heated by SNII to $10^6~\ensuremath{{\rm K}}$. The halo gas it at $10^6{-}10^7~\ensuremath{{\rm K}}$.

Open with DEXTER

The evolution of the cloudy phase (Fig. 4) is influenced by dissipation and SF. The depletion of the cloudy medium by SF decreases its mean surface density. Due to inelastic cloud-cloud collisions, the clouds also show a less smooth distribution. As a consequence a mass flow towards the centre increases its central surface density. The radial profiles reveal that the surface density has decreased mostly around  $5~\ensuremath{{\rm kpc}}$, while it it almost unchanged beyond  $10~\ensuremath{{\rm kpc}}$. The cloud disk shows a thickening only outside of  $11~\ensuremath{{\rm kpc}}$ where the increase of  $z_{\rm 95,cld}$ is very similar to that of the stellar disk. Inside that radius no thickening can be seen. In fact, due to dissipation  $z_{\rm 95,cld}$ is even reduced inside $5~\ensuremath{{\rm kpc}}$ by a factor of 2-3.

The diffuse gas (Fig. 5) starts from a homogeneous distribution that collapses in less than  $300~\ensuremath{{\rm Myr}}$ to form a disk. Initially, the gaseous disk is only a few kpc in diameter but it slowly grows in size. After  $1~\ensuremath{{\rm Gyr}}$ it has reached a diameter of  $20~\ensuremath{{\rm kpc}}$, while at 2 gyr the diameter is $\sim$ $30~\ensuremath{{\rm kpc}}$, which is comparable to that of the stellar disk. At the same time the mass of the gaseous disk is nearly constant. It has particle densities of $0.1{-}1~\ensuremath{{\rm cm}^{-3}}$ and temperatures (Fig. 6) of less than  $10^4~\ensuremath{{\rm K}}$. Towards the centre the gas reaches higher temperatures of up to several $10^6~\ensuremath{{\rm K}}$. Some diffuse gas remains in the halo at densities of $\sim$ $10^{-4}~\ensuremath{{\rm cm}^{-3}}$ and less, while its temperature is in the range of $10^6{-}10^7~\ensuremath{{\rm K}}$. During the initial collapse, the diffuse gas gains a significant amount of angular momentum transferred from the clouds by the drag force. For this reason the gaseous disk rotates with a rotation velocity of $\sim$ $200~\ensuremath{{\rm km}~{\rm s}^{-1}}$. Transient spiral patterns are also observed in the gas disk. Stellar feedback induces the emergence of bubbles of hot gas. These bubbles have diameters of up to $1~\ensuremath{{\rm kpc}}$ and they cool within a few $10~\ensuremath{{\rm Myr}}$. The volume-filling factor of the hot gas ( $T>10^4~\ensuremath{{\rm K}}$) in the disk is 30-40%.

The velocity dispersion of the stars increases with time. At  $8~\ensuremath{{\rm kpc}}$the heating rate is about  $1.3~\ensuremath{{\rm km}~{\rm s}^{-1}} ~\ensuremath{{\rm Gyr}^{-1}}$, and it is nearly constant over the whole disk. In order to see how much heating is caused by unphysical two-body relaxation (which stems from the inability to resolve individual stars numerically), we compared the result with a pure N-body simulation and found that the heating rate in our full model is about 30% higher than in the pure N-body calculation. The velocity dispersion of the clouds decreases at  $0.6~\ensuremath{{\rm km}~{\rm s}^{-1}} ~\ensuremath{{\rm Gyr}^{-1}}$ due to the energy dissipation by cloud-cloud collisions. As a result younger stars are born with a smaller velocity dispersion. This means that a part of the observed "heating'' in the age-velocity-relation might indeed be due to cooling as suggested by Rocha-Pinto et al. (2004).

	Figure 7: The evolution of the SFR with time (full line). The SFR decreases with time due to consumption of cloud mass. The average SFR (dashed line) is $\approx$ $1.6~\ensuremath{{M}_{\odot}} ~\ensuremath{~{\rm yr}^{-1}}$. The dotted line is the result of a simple analytic model accounting for the depletion of the cloudy medium by SF.
Open with DEXTER

	Figure 8: The evolution of SF efficiency. SFE is up to 0.4 in the centre of the galaxy. The average over the whole disk is about 0.06.
Open with DEXTER

	Figure 9: The open diamonds show the SFR per $\ensuremath {{\rm pc}^2}$ averaged over $100~\ensuremath{{\rm Myr}}$ as a function of cloud surface density. The full line is a recent determination of the Schmidt Law by Kennicutt (1998) with a slope of -1.4. The vertical dashed line indicates a cut-off density for SF in disk galaxies that was also found by Kennicutt (1998).
Open with DEXTER

The average SFR is $\sim$ $1.6~\ensuremath{{M}_{\odot}} ~\ensuremath{~{\rm yr}^{-1}}$ (Fig. 7). Over the $3~\ensuremath{{\rm Gyr}}$ of evolution the SFR decreases slowly from $\sim$ $4.0~\ensuremath{{M}_{\odot}} ~\ensuremath{~{\rm yr}^{-1}}$ to $\sim$ $0.5~\ensuremath{{M}_{\odot}} ~\ensuremath{~{\rm yr}^{-1}}$. A simple analytic model indicates that this is due to the consumption of the clouds. The mean SFE is roughly constant at a level of $\epsilon = 0.06$. The maximum values of the local SFE are highest in the centre with $\epsilon = 0.2{-}0.4$ (Fig. 8). It can be seen in Fig. 9 that the SFR also compares well to a recent determination of the Schmidt law by Kennicutt (1998), who has also shown that the SF law in disk galaxies steepens abruptly below a threshold gas density of $\sim$ $10~\ensuremath{{M}_{\odot}} ~\ensuremath{{\rm pc}^{-2}}$. We also find in our simulations that the SF law steepens for lower gas densities. This agrees with 1d CD models of the settling of a gas disk, which show a similar threshold for the formation of a thin disk (Burkert et al. 1992). Excluding the data below $10~\ensuremath{{M}_{\odot}} ~\ensuremath{{\rm pc}^{-2}}$, we obtain a least-square fit of our data to Eq. (6) with $n \approx 1.7$ $\pm$ 0.1, which is in fair agreement with $n \approx 1.4$ $\pm$ 0.15 of Kennicutt (1998).

The mass exchange by condensation and evaporation reaches a quasi-equilibrium after $\sim$ $200~\ensuremath{{\rm Myr}}$. Before, evaporation is the dominant process due to the heating of the diffuse gas induced by the initial collapse of this component. In the later evolution, the condensation rate decreases from $\sim$ $2~\ensuremath{{M}_{\odot}} ~\ensuremath{~{\rm yr}^{-1}}$ to $0.1~\ensuremath{{M}_{\odot}} ~\ensuremath{~{\rm yr}^{-1}}$. The evaporation rate is slightly lower after  $1.5~\ensuremath{{\rm Gyr}}$ due to the decreasing SN rate, resulting in a reduced heating of the diffuse ISM.

The total cloud mass slowly decreases, showing the depletion by SF. The total mass of the diffuse gas remains almost constant. A detailed look shows that more than half of the diffuse gas (65%) is located in the disk ( $R<20~{\rm kpc}$ and $z<1~{\rm kpc}$), while the rest is in the halo.

In the initial model all clouds have the same mass. However, due to fragmentation and coalescence, a stationary cloud mass spectrum is built up within $1~\ensuremath{{\rm Gyr}}$ (Fig. 10). This mass spectrum can be fitted by a power-law with index $\alpha \approx -2$ with an uncertainty of 0.2. Once the mass spectrum is established it remains unchanged until the end of the simulation.

	Figure 10: The cloud mass spectrum. The different lines correspond to the mass spectrum after $1~\ensuremath{{\rm Gyr}}$ (dotted), $2~\ensuremath{{\rm Gyr}}$ (dashed), and $3~\ensuremath{{\rm Gyr}}$ (full) of evolution. The fitted power-law mass function is with $\alpha \approx -2.0$ comparable to the observed mass spectrum.
Open with DEXTER

   
5 Discussion

In the previous section we presented a model for a Milky-Way-type galaxy. Since our description depends on several phase transitions and exchange processes, we want to discuss their influence in this section, looking at SF, cloud mass spectrum, and metallicities.

The model has a number of parameters that can be chosen within the constraints of related observations or theory. We therefore performed simulations to investigate the influence of some of those parameters using the model presented in Sect. 4.2 as a reference. We find that most models do not deviate much from the reference model.

5.1 Star formation

The mean SFR is $\sim$ $1.6~\ensuremath{{M}_{\odot}} ~\ensuremath{~{\rm yr}^{-1}}$ in most models with a mean SF efficiency of $\sim$$6\%$. However, the SFR decreases in all models from a few to below one solar mass per year. Observations, on the other hand, reveal a much more constant SF history (Rocha-Pinto et al. 2000). The decrease in the SFR can be reproduced in an analytic model taking the depletion of the cloudy medium by SF into account. An explanation for the discrepancy might be that our simulation does not include any infall to sustain a constant SFR. An infall rate of the order of $\sim$ $1~\ensuremath{{M}_{\odot}} ~\ensuremath{~{\rm yr}^{-1}}$ would be sufficient for this, provided the infalling gas can be used in the SF process. The necessary infall rate could be provided by high velocity clouds (Wakker et al. 1999; Blitz et al. 1999).

The SF law agrees well with observations. Our data can be fitted to Eq. (6) with $n \approx 1.7$, which is well within the observed range of n=1.4-2.0 (Kennicutt 1998; Wong & Blitz 2002; Boissier et al. 2003). It should be noted that this result is widely independent of the choice of our model parameters. Only models with a constant SF efficiency show a different behaviour with n=1.2. Köppen et al. (1995) have shown that stellar feedback leads to a self-regulated SF following a Schmidt-type law, and our simulations confirm this result. We conclude that the observed index n of the SF law is a consequence of the self-regulation process that is enabled by allowing for a variable SFE.

   
5.2 Cloud mass function

The cloud mass function in our reference model has a slope of $\alpha = -2$. In most observations the mass spectrum is shallower with $\alpha \approx -1.5$ (e.g. Sanders et al. 1985), although there are a few observations with higher values, e.g.  $\alpha \approx -1.8$ (e.g. Brand & Wouterloot 1995). Taking the uncertainties in both our model and observations into account we are able to resolve the cloudy phase of the ISM; i.e. each cloud particle in our simulation represents an individual molecular cloud.

Theoretical models of coagulation by Field & Saslaw (1965) also predict $\alpha = -1.5$. These models include a fragmentation by SF, but (different to our model) SF occurs only when clouds reach a maximum mass. The resulting fragments are then added to the lowest mass bin. In a more similar model, which also takes a mass-radius relation (Eq. (1)) into account, Elmegreen (1989) found a mass spectrum with $\alpha = -2$, showing that our simulations can reproduce these results.

   
5.3 Metallicity

We compared two models with metallicities $Z=1/2\ Z_\odot$ and $Z=1.5\ Z_\odot$. The change of metallicity is only moderate and does not affect the dynamical evolution of the galaxy significantly. This result verifies that using a constant metallicity is a valid approximation on short time scales. Using a simple estimate we find that for  $3~\ensuremath{{\rm Gyr}}$ the expected change in metallicity would be of the order $1/10~Z_\odot$, i.e. well within the range given by the two test models. Therefore, a model that includes stellar yields would not change the results, as long as the late evolution (i.e. the last few Gyr) of galaxies is simulated. However, stellar yields would have to be included when the early formation stages are considered.

5.4 The multi-phase ISM

The ISM in our models can be divided into three phases: molecular clouds and diffuse components that can be either warm or hot. The total gas mass fraction after  $3~\ensuremath{{\rm Gyr}}$ of evolution is about $7\%$ in all our models. The ratio of the mass of clouds to the mass of diffuse gas is mainly determined by condensation and evaporation. If the evaporation rate is higher, for example due to higher feedback energy, this ratio is closer to unity. Between $80\%$ and $90\%$ of the diffuse gas is found in the warm gaseous disk, while only $\sim$$15\%$ remain in the halo in a hot gas phase. In the low-metallicity model this fraction is larger by a factor of 2 due to less efficient cooling. Samland & Gerhard (2003) have also presented a model of a disk galaxy, where the ratio of cloud and diffuse gas mass is 3, which is comparable to our results. They also find a hot gaseous halo. In a similar calculation, Semelin & Combes (2002) find equal masses in clouds and gas independent of the feedback energy, but their models do not include phase transitions between the gas phases by condensation and evaporation.

We found a filling factor f of 0.3-0.4 for the hot gas in the disk. McKee & Ostriker (1977) predicted 0.7, a higher value, but our value agrees better with the more recent findings of Norman & Ikeuchi (1989), who predicted a filling factor of about 0.2. Furthermore, three-dimensional simulations of a region of the galactic disk by de Avillez (2000) show a filling factor for gas with $T>10^4~\ensuremath{{\rm K}}$ of $\sim$0.4 in good agreement with our results. We find an increase of 0.2 in f if the energy injection by SNII is higher by a factor 4, which is comparable to the results of de Avillez & Breitschwerdt (2004).

One shortcoming in our treatment of the ISM is the lack of a detailed cooling below  $10^4~\ensuremath{{\rm K}}$. As a consequence, we do not allow clouds to form directly by thermal instabilities. However, due to low mean densities and the turbulent state of the diffuse gas, we expect the formation of clouds in expanding SN shells (which is included) to be much more important. Again, a more sophisticated treatment of thermal instabilities would be required for modelling the early evolution of (gas rich) galaxies.

5.5 Live DM halo

So far, in all presented models the DM halo is treated as a static potential. A static halo is reasonable for the case of the late evolution of an isolated galaxy. However, DM particles have to be used, e.g., for interacting galaxies or early galactic evolution. In one model we therefore used a DM halo made of particles in order to test how a live halo changes our results. The DM halo was modelled with 10⁵ particles, thus DM particles are twice as massive as the initial masses of star particles. As expected we do not find significant changes in the evolution when using a live DM halo. The CPU-time needed in this run is higher by a factor of 2 (10 cpu-days instead of 5).

   
6 Summary and conclusions

We have presented a new particle-based code for modelling the evolution of galaxies. In addition to the collisionless stellar and DM components, it incorporated a multi-phase description of the ISM combining standard SPH with a sticky particle method to describe diffuse warm/hot gas and cold molecular clouds, respectively. We included a new description for SF with a local SF efficiency, depending on the local pressure of the ISM and the mass of star-forming clouds. Stellar life times and an IMF are taken into account, along with SNII and PNe. The gas phases can mix by means of condensation and evaporation, and clouds are subject to a drag force.

We use this description to model $3~\ensuremath{{\rm Gyr}}$ of evolution for a Milky-Way-type galaxy. Overall, the evolution of our model galaxy is quite stable. Furthermore, a parameter study showed that the evolution does not change significantly when varying model parameters like feedback energy, cloud radius, SF time scale or metallicity in a reasonable range.

Our models produced an SF law with an index n=1.7 that agrees well with observations. For models with a constant SF efficiency $\epsilon$, the index is with n=1.2 smaller than observations suggest. We conclude, in agreement with Köppen et al. (1995), that the SF law is the result of a self-regulation process enabled by our new treatment of SF. The average SFR is also comparable to observations. Over time the SFR decreases by a factor of 4-8 because we have not included any gas infall. To get a more constant SFR, a gas infall should be taken into account at a rate of $\sim$ $1~\ensuremath{{M}_{\odot}} ~\ensuremath{~{\rm yr}^{-1}}$. The cloud mass spectrum in our models (with $\alpha = -2$) is a little bit steeper than most observations suggest. Nonetheless, each cloud particle in our code essentially represents an individual molecular cloud.

In all models most (85%) of the diffuse gas forms a warm ($\sim$ $10^4~\ensuremath{{\rm K}}$) rotating disk. A hot gaseous halo with temperatures of  $10^6{-}10^7~\ensuremath{{\rm K}}$ and densities of $10^{-4}~\ensuremath{{\rm cm}^{-3}}$ has been built up. Hot gas is also found in the disk in bubbles accompanied by SNII. The filling factor of the hot disk gas is about 0.4, comparing well to observations as well as models of the ISM. This result could be used to find more realistic conditions for the initial model.

In a next step we will apply our model to interacting galaxies in order to follow the evolution of the ISM and the SF history in perturbed galactic systems. The chemical evolution should be embedded to allow for a more detailed confrontation with observations.

Acknowledgements

The authors would like to thank P. Berczik for fruitful discussions and the anonymous referee for her/his helpful comments that improved this manuscript. We thank W. Dehnen (falcON), K. Kujken & J. Dubinski (GalactICS), and P. Teuben (NEMO) for making their software publicly available. This work was supported by the German Deutsche Forschungsgemeinschaft, DFG project number TH 511/2. S. Harfst is grateful for financial support from the Australian Department of Education, Science and Training (DEST).

References

• Aarseth, S. J. 1985, in Multiple Time Scales, ed. J. U. Brackbill, & B. I. Cohen, 377 (In the text)
• Alongi, M., Bertelli, G., Bressan, A., et al. 1993, A&AS, 97, 851 [NASA ADS]
• Barnes, J. E. 2004, MNRAS, 350, 798 [NASA ADS] [CrossRef] (In the text)
• Barnes, J., & Hut, P. 1986, Nature, 324, 446 [NASA ADS] [CrossRef] (In the text)
• Bertelli, G., Bressan, A., Chiosi, C., Fagotto, F., & Nasi, E. 1994, A&AS, 106, 275 [NASA ADS]
• Blitz, L., Spergel, D. N., Teuben, P. J., Hartmann, D., & Burton, W. B. 1999, ApJ, 514, 818 [NASA ADS] [CrossRef]
• Böhringer, H., & Hensler, G. 1989, A&A, 215, 147 [NASA ADS] (In the text)
• Boissier, S., Prantzos, N., Boselli, A., & Gavazzi, G. 2003, MNRAS, 346, 1215 [NASA ADS] [CrossRef]
• Brand, J., & Wouterloot, J. G. A. 1995, A&A, 303, 851 [NASA ADS] (In the text)
• Bressan, A., Fagotto, F., Bertelli, G., & Chiosi, C. 1993, A&AS, 100, 647 [NASA ADS]
• Brown, J. H., Burkert, A., & Truran, J. W. 1995, ApJ, 440, 666 [NASA ADS] [CrossRef]
• Burkert, A., Truran, J. W., & Hensler, G. 1992, ApJ, 391, 651 [NASA ADS] [CrossRef]
• Carlberg, R. G. 1984, ApJ, 286, 416 [NASA ADS] [CrossRef]
• Combes, F., & Gerin, M. 1985, A&A, 150, 327 [NASA ADS]
• Cowie, L. L., McKee, C. F., & Ostriker, J. P. 1981, ApJ, 247, 908 [NASA ADS] [CrossRef] (In the text)
• de Avillez, M. A. 2000, MNRAS, 315, 479 [NASA ADS] [CrossRef] (In the text)
• de Avillez, M. A., & Breitschwerdt, D. 2004, A&A, 425, 899 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Dehnen, W. 2002, JCP, 179, 27 [NASA ADS] (In the text)
• Elmegreen, B. G. 1989, ApJ, 347, 859 [NASA ADS] [CrossRef] (In the text)
• Elmegreen, B. G. 2000, ApJ, 530, 277 [NASA ADS] [CrossRef]
• Elmegreen, B. G., & Efremov, Y. N. 1997, ApJ, 480, 235 [NASA ADS] [CrossRef]
• Evrard, A. E. 1988, MNRAS, 235, 911 [NASA ADS] (In the text)
• Field, G. B., & Saslaw, W. C. 1965, ApJ, 142, 568 [NASA ADS] [CrossRef] (In the text)
• Gingold, R. A., & Monaghan, J. J. 1983, MNRAS, 204, 715 [NASA ADS] (In the text)
• Kennicutt, R. C. 1998, ApJ, 498, 541 [NASA ADS] [CrossRef] (In the text)
• Köppen, J., Theis, C., & Hensler, G. 1995, A&A, 296, 99 [NASA ADS] (In the text)
• Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545 [NASA ADS] (In the text)
• Kuijken, K., & Dubinski, J. 1995, MNRAS, 277, 1341 [NASA ADS] (In the text)
• McKee, C. F., & Ostriker, J. P. 1977, ApJ, 218, 148 [NASA ADS] [CrossRef] (In the text)
• Mihos, J. C., Richstone, D. O., & Bothun, G. D. 1992, ApJ, 400, 153 [NASA ADS] [CrossRef]
• Mihos, J. C., Bothun, G. D., & Richstone, D. O. 1993, ApJ, 418, 82 [NASA ADS] [CrossRef]
• Monaghan, J. J. 1992, ARA&A, 30, 543 [NASA ADS] (In the text)
• Monaghan, J. J., & Lattanzio, J. C. 1985, A&A, 149, 135 [NASA ADS] (In the text)
• Noguchi, M. 1988, A&A, 203, 259 [NASA ADS]
• Norman, C. A., & Ikeuchi, S. 1989, ApJ, 345, 372 [NASA ADS] [CrossRef] (In the text)
• Palous, J., Jungwiert, B., & Kopecky, J. 1993, A&A, 274, 189 [NASA ADS]
• Pearce, F. R., Jenkins, A., Frenk, C. S., et al. 1999, ApJ, 521, L99 [NASA ADS] [CrossRef]
• Polyachenko, V. L., Polyachenko, E. V., & Strel'Nikov, A. V. 1997, Astro. Lett., 23, 483 [NASA ADS] (In the text)
• Raiteri, C. M., Villata, M., & Navarro, J. F. 1996, A&A, 315, 105 [NASA ADS] (In the text)
• Rivolo, A. R., & Solomon, P. M. 1988, Lecture Notes in Physics (Berlin: Springer Verlag), 315, 42 [NASA ADS] (In the text)
• Rocha-Pinto, H. J., Scalo, J., Maciel, W. J., & Flynn, C. 2000, A&A, 358, 869 [NASA ADS] (In the text)
• Rocha-Pinto, H. J., Flynn, C., Scalo, J., et al. 2004, A&A, 423, 517 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Samland, M., Hensler, G., & Theis, C. 1997, ApJ, 476, 544 [NASA ADS] [CrossRef]
• Samland, M., & Gerhard, O. E. 2003, A&A, 399, 961 [EDP Sciences] [NASA ADS] [CrossRef]
• Sanders, D. B., Scoville, N. Z., & Solomon, P. M. 1985, ApJ, 289, 373 [NASA ADS] [CrossRef] (In the text)
• Schmidt, M. 1959, ApJ, 129, 243 [NASA ADS] [CrossRef] (In the text)
• Semelin, B., & Combes, F. 2002, A&A, 388, 826 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Shu, F. H., Milione, V., Gebel, W., et al. 1972, ApJ, 173, 557 [NASA ADS] [CrossRef] (In the text)
• Teuben, P. 1995, in Astronomical Data Analysis Software and Systems IV, ASP Conf. Ser., 77, 398 (In the text)
• Thacker, R. J., Tittley, E. R., Pearce, F. R., Couchman, H. M. P., & Thomas, P. A. 2000, MNRAS, 319, 619 [NASA ADS] [CrossRef]
• Theis, C., & Hensler, G. 1993, A&A, 280, 85 [NASA ADS]
• Theis, C., Burkert, A., & Hensler, G. 1992, A&A, 265, 465 [NASA ADS]
• Theis, C., Ehlerová, S., Palous , J., & Hensler, G. 1998, Lecture Notes in Physics (Berlin: Springer Verlag), 506, 409 [NASA ADS]
• Wakker, B. P., Howk, J. C., Savage, B. D., et al. 1999, Nature, 402, 388 [NASA ADS] [CrossRef]
• Weidemann, V. 2000, A&A, 363, 647 [NASA ADS] (In the text)
• Wong, T., & Blitz, L. 2002, ApJ, 569, 157 [NASA ADS] [CrossRef]

 

  
Online Material

   
Appendix A: Determination of the star formation efficiency

Elmegreen & Efremov (1997) derived an implicit relation between cloud mass , pressure  and SFE $\epsilon$ ( ${M}_5 = 10^5~\ensuremath{{M}_{\odot}}$)

$$% 1 = \epsilon + \frac{1}{2}A_0\left(\frac{\ensuremath{m_{\rm... ...}\right)^{-5/8}\epsilon^2, {\rm when}\ \frac{\epsilon}{\xi}<1,$$ (A.1)

and

$$% \begin{array}{ll} 1 = \epsilon\ & +\ A_0\left(\displaysty... ... \; {\rm when} \displaystyle\frac{\epsilon}{\xi}>1. \end{array}$$ (A.2)

They use $\ensuremath{P_{\odot}} = \ensuremath{3 \times 10^{4}} ~\ensuremath{{\rm cm}^{-3}} ~\ensuremath{{\rm K}}$ for the pressure in the solar neighbourhood, and $\xi$ is defined as

$$% \xi \equiv \xi_0 \left(\frac{\ensuremath{m_{\rm cld}} }{{\r... ...math{P_{\rm gas}} }{\ensuremath{P_{\odot}} }\right)^{3/8}\cdot$$ (A.3)

The constants A₀ = 180 and $\xi_0 = 0.1$ are determined from observations of star-forming clouds in the solar neighbourhood.

In each SF event Eqs. (A.1) and (A.2) are solved numerically to determine $\epsilon$. The resulting SFE is highest for high-mass clouds in a high-pressure environment. Typical SFE are between 5% and 10%.

   
Appendix B: SPH implementation

The idea of SPH is that a fluid is represented by N particles. Any physical property A of the fluid at the position  ${\vec r}_p$is evaluated by means of a locally weighted average:

$$% A({\vec r}_p) = \sum_q m_q {A_q \over \rho_q} W\left({\vec r}_{pq},h\right).$$ (B.1)

The summation is performed over particles q with mass m_q, density $\rho_q$, while A_q is the value of A at  ${\vec r}_q$ and the distance  ${\vec r}_{pq}$ is ${\vec r}_p - {\vec r}_q$. We apply the spherically symmetric smoothing kernel W(r,h) suggested by Monaghan & Lattanzio (1985)

$$% W(r,h) = {1 \over \pi h^3} \left\{ \begin{array}{lllll} ... ... 0, & \mbox{ for ${r\over h} > 2$ .} \\ \end{array} \right.$$ (B.2)

A variable smoothing length h is applied to keep the number  $N_{\rm sm}$ of nearest neighbours nearly constant at about 32 for all SPH particles.

The equation of motion of an SPH particle p is given by

$$% {{\rm d}^2{\vec r}_p \over {\rm d}t^2} = - {1 \over \rho_p} \nabla P_p + {\vec a}_p^{\rm visc} - \nabla \Phi_p,$$ (B.3)

where $\Phi_p$ is the gravitational potential, P_p the pressure, and  ${\vec a}_p^{\rm visc}$ an artificial viscosity term to handle shocks. We use a symmetric expression to estimate $(\nabla P)/\rho$:

$$% {\nabla P_ p \over \rho_p} = \sum_q m_q \left( {P_p \over \rho_p^2} + {P_q \over \rho_q^2}\right) \nabla W_{pq},$$ (B.4)

where W_pq indicates that the arithmetic mean of the kernel function for the two individual smoothing lengths h_p and h_q is used. For the artificial viscosity term  ${\vec a}_p^{\rm visc}$, we use (Gingold & Monaghan 1983)

$$% {\vec a}_p^{\rm visc} = - \sum_q m_q \Pi_{pq} \nabla W_{pq},$$ (B.5)

with

$$% \Pi_{pq} = { -\alpha_{\rm v}\mu_{pq}c_{pq} + \beta_{\rm v}\mu_{pq}^2\over \rho_{pq}},$$ (B.6)

where $\rho_{pq}$ and c_pq are the arithmetic means of density and sound speed of particles p and q. Then $\mu_{pq}$ is defined by

$$% \mu_{pq} = \left\{ \begin{array}{ll} {\displaystyle{\vec v}_{pq... ... & \mbox{for ${\vec v}_{pq}\cdot{\vec r}_{pq} \geq 0$ } \\ \end{array}\right.,$$ (B.7)

where ${\vec v}_{pq}$ is the relative velocity of particle p and q. The viscosity parameters are set to  $\alpha_{\rm v}=2$, $\beta_{\rm v} = 1$, and $\eta_{\rm v}=0.1$. The evolution of the specific internal energy u is computed by

$$% {{\rm d}u_p \over {\rm d}t} = {1\over2}\sum_q m_q \left({P_... ...}_{pq}\cdot\nabla W_{pq} + {\Gamma - \Lambda \over \rho}\cdot$$ (B.8)

Here, $\Gamma$ and $\Lambda$ account for heating and cooling processes.

The time integration is done by a leap-frog integrator. Individual time steps are used for each particle. The time step for an SPH particle p is chosen to satisfy a modified form of the Courant condition as suggested by Monaghan (1992)

$$% \Delta t_p \le C { h_{p} \over h_{p}\left\vert\nabla\cdot{\... ... \beta_{\rm v} \max_q \left\vert \mu_{pq}\right\vert \right)},$$ (B.9)

where C = 0.4 is the Courant number, h_p the smoothing length, c_p the sound speed, and  ${\vec v}_p$ the velocity of particle p.

Because cooling time scales are usually much shorter than the dynamical times, impractically short time steps might occur. Therefore, Eq. (B.8) is often solved semi-implicitly. However, we chose a different approach: a fourth-order Runge-Kutta scheme is used to integrate the cooling function beforehand with a sufficiently small time step, and a look-up table with times and temperatures is created. The number density  $\ensuremath{n_{\rm gas}}$ is fixed to one value for this, but cooling times scale linearly with  $\ensuremath{n_{\rm gas}}$ and can be easily rescaled. Starting from a given temperature, the new temperature after a (rescaled) time step is then simply looked up on the table. This provides a fast way to compute the change in temperature, i.e. the cooling rate, for any given temperature, density, and time step. The cooling rate $\Lambda$ for each particle is updated at each individual time step.

A more detailed description of the SPH method can be found in Monaghan (1992).