Innovatives Supercomputing in Deutschland
inSiDE • Vol. 11 No. 2 • Autumn 2013
current edition
about inSiDE
index  index prev  prev next  next

Simulating the Life Cycle of molecular Clouds

Star formation takes place in the densest and coldest gas in a star-forming galaxy, in so-called molecular clouds (MCs). MCs do not evolve in isolation. They are highly dynamical objects, which are born, fed, heated, and stirred from their turbulent environment into which they eventually dissolve. They form in regions where the hot or warm, ionized and atomic interstellar medium (ISM) condenses into cold (T< 300K), molecular gas. Often concentrated to the midplane of galactic disks, this process involves metallicity-dependent, nonequilibrium chemistry and molecule formation, heating and cooling, turbulence, self-gravity, and magnetic fields. Once formed, MCs further collapse to form stars and star clusters.

Observations in the Milky Way and nearby galaxies show that individual MCs convert only a few percent of their mass into stars [5]. Less than 1% of all new-born stars are more massive than 8 solar masses, but these are particularly important for galaxy formation and evolution. The life and death of massive stars differ intriguingly from those of their low-mass counterparts. Such stars affect their environment dramatically through their strong UV radiation, their energetic stellar wind, and their final explosion as supernovae. These ’feedback’ processes generate turbulence in the parental molecular cloud, dissociate, ionize, and eventually destroy them from within, thereby quenching further star formation. Stellar feedback is thus thought to regulate the star formation efficiency in molecular clouds leading to a self-regulation of star formation on galactic scales (see [9] and references therein).

At the same time feedback from massive stars can drive powerful galactic outflows, removing gas from the host galaxy. In this way star formation may control the galactic gas reservoir available for star formation and regulate galaxy formation on cosmic scales. In support of this idea, high star formation rates linked to powerful galactic outflows were recently observed for gas-rich, star-forming disk galaxies at the peak of cosmic star formation at z=2-3 [12]. Galactic winds may also be of particular importance for driving the Galactic fountain [10] and removing gas from dwarf galaxies which, today, have a lower than average ratio of baryons to dark matter [11].

Understanding how feedback effects propagate from the milli-parsec scales, characterizing the sites of massive star formation, to the kilo-parsec scales of molecular cloud complexes and galactic winds is a physical and numerical challenge. In the framework of the Gauss Tier-0 project “SILCC” we model representative regions of disk galaxies using adaptive, three-dimensional simulations at unprecedented resolution and with the necessary physical complexity to follow the full life-cycle of molecular clouds. These simulations include self-gravity, magnetic fields, heating and cooling at different gas metallicities, molecule formation and dissociation, and stellar feedback in form of stellar winds and supernovae. The simulations use novel, massively parallel techniques, and are computationally very demanding. They are only feasible on a peta-scale supercomputer like SuperMUC at the Leibniz Supercomputing Centre (LRZ) in Garching near Munich.

The ultimate goal of the SILCC project is to provide a self-consistent answer as to how stellar feedback regulates the star formation efficiency of a galaxy, how molecular clouds are formed and destroyed, how galactic outflows are driven and what parameters their mass-loading depends upon. Here we discuss first results showing that the fraction of gas, which is in molecular form, is steadily increasing with the total gas surface density ∑gas of the disk.

Numerical Method

The simulations are being carried out with the massively parallel, Eulerian, adaptive mesh refinement code FLASH 4.0.1 [6]. FLASH is highly modular and includes many physics units, like MHD, an external gravitational potential, or recently also a Barnes & Hut tree method for self-gravity, in the standard release. However, much more physics is required in order to carry out the proposed simulations. In preparation for this project we added new physics modules, which are not present in the standard release.

In particular, we implemented a chemical network that allows us to follow all three phases in the ISM, and which includes non-equilibrium ionization effects as well as H2 and CO formation [7]. Also radiative heating and cooling are computed by this module. The heating and cooling terms, and also the chemical source and sink terms, are operator split from the advection and pdV work terms and are solved implicitly using the DVODE solver [2]. In cases where the chemical or cooling timesteps are much shorter than the hydrodynamic timestep, sub-cycling is used to treat the cooling and chemistry, thereby avoiding the need to constrain the global timestep. Chemical abundances are represented as tracer fields and are advected using the standard FLASH infrastructure for such fields.

In order to treat the formation of molecules, which is essential for comparing the simulations with observations, we implemented TreeCol [4] to compute shielding and self-shielding of H2 and CO as well as the dust attenuation and dust temperature. TreeCol is a cost-efficient way to treat diffuse radiation. From every cell within the computational domain, a set of rays is casted along which the e.g. column densities are computed during the self-gravity tree-walk. The ray-casting algorithm is based on HealPix [8], an equal surface area pixilation of the sphere. Therefore, the minimum number of rays is 12, and we typically use 48 rays.

Furthermore, we implemented Supernova and stellar wind feedback. Supernovae can be placed randomly with a given Supernova rate as used in the presented simulations. Alternatively, the explosions might be launched from gas density peaks within molecular clouds or from sink particles. Each Supernova injects 1051 erg and a few solar masses of gas into the ISM, which we add in form of internal energy and additional mass to cells close to the explosion center, i.e. within the injection region. Before the explosion, the injection region is always refined to the finest AMR level in order to limit grid effects and over-cooling. We force a minimum number of 4 cells per radius of the injection region, but the radius may be larger if the mass within the injection region is small (≤ 30 Msun). With these parameters the temperature within the injection region is typically raised to ~107 K, and therefore the hydrodynamic timestep is recomputed and limited accordingly. Stellar winds are implemented in a similar way.

Last but not least, we use sink particles with an adequate sub-grid model to describe star clusters for simulations with lower resolution, or small N-clusters of accreting stars in case of high resolution.

We simulate representative regions of disk galaxies (500 pc x 500 pc x ± 2 kpc) with different gas surface densities, i.e. 10, 30, 100, and 500 Msun/pc2, to mimic the conditions of star-forming galaxies at low and high redshift. The initial density distribution of the gaseous disk is Gaussian. The boundary conditions are periodic in x- and y-, and allow for outflow in z-direction. To account for the aforementioned physical processes on spatial and temporal scales ranging from several 100 AU and 10 years to kpc-scales and 107 years, each run requires more than 5 million CPU hours.

The typical Milky Way Disk

Close to the solar radius (at a distance of 8.5 kpc from the Galactic Center) the Milky Way disk has an average gas surface density of 10 Msun/pc2. We present first results from simulations at an intermediate resolution, where we include turbulence driving from randomly placed Supernova explosions. The random positions are weighted with a Gaussian with a scale height of 50 pc. Assuming a certain stellar initial mass function (IMF; e.g. [3]) to estimate the number of massive stars per unit mass that is forming stars, the mean Supernova rate (SNR) for such a disk can be derived from the observed star formation rate – gas surface density relation. Observations yield nearly linear scaling of ∑SFR ~ ∑gas [1]. Within the simulated volume, the SNR for a Milky Way-type disk is ~10 Myr-1. Here we use a SNR of 14 Myr-1.

The simulations show that a dynamical equilibrium is only established if the SNR does not significantly exceed or fall below the observed SNR. In Fig. 1, we show projections and slices through the (y=0)-plane. In the top row we show (a) total column density, (b) density slice, (c) temperature slice, and (d) dust temperature slice; in the bottom row we plot the column densities of (e) atomic hydrogen (HI), (e) ionized hydrogen (HII), (f) molecular hydrogen (H2), and (g) CO. In the disk mid-plane the gas is fully molecular unless a Supernova explodes within or nearby a molecular cloud. The point-like explosions heat up the surrounding gas to 10 million Kelvin. The expanding blast wave then compresses the ISM into a denser shell, within which the dust temperature is temporarily increased. On the other hand, new molecular gas forms in very cold (< 10 K) molecular clouds within the midplane. Within the molecular clouds the gas is fully molecular and the abundances of H2 and CO are well correlated. However, there is molecular hydrogen at higher altitudes, which is not well traced by CO. This gas, which has been driven out of the molecular clouds by Supernovae, is CO-dark and cannot be seen in observations. The reason that CO does not survive at higher altitudes is that it requires a higher column density to shield itself from the interstellar (ultraviolet) radiation field.

Figure 1: Projections and slices through the y=0 plane at t=20 Myr. From top left to bottom right the panels show (a) the total column density, (b) a density slice, (c) a temperature slice, (d) a dust temperature slice, and the column density of (e) atomic hydrogen (HI), (f) ionized hydrogen (HII), (g) molecular hydrogen (H2), and (h) CO.

Disk Galaxies at different Gas Surface Densities: from low to high Redshift

Simulations of disks with higher gas surface density show that the mass fraction molecular gas is increasing with increasing surface density. This is the case even though the Supernova rate is correspondingly higher in high surface density disks, following a roughly linear scaling relation of star formation rate (hence Supernova rate) and molecular gas surface density [13].

In Fig. 2, we show the molecular gas mass fraction as a function of time for the four different initial disk surface densities: ∑gas =10, 30, 100, and 500 Msun/pc2. The corresponding SNRs are 14, 42, 140, and 700 Myr-1.

For the highest gas surface density of 500 Msun/pc2, the H2 mass fraction rapidly approaches values of up to 90 %. At 100 Msun/pc2 the H2 mass fraction saturates at 70 – 80 %. This result is in good agreement with observational estimates of the molecular gas mass fraction in centers of local galaxies as well as gas-rich high redshift disks [13]. For lower initial gas surface densities the molecular fractions are respectively lower and it takes longer to reach the saturation values. The Milky Way value of roughly 50 % (interior to the Sun’s galactic orbit) agrees well with the 10 Msun/pc2 disk simulation, for which peak values of 40 % are achieved.

Even though the mass fractions are high the fractional volume filled by molecular gas is small, as it represents the densest and coldest gas component. Therefore, simulating the life-cycle of molecular clouds requires simulations with very high spatial resolution, that are currently performed on SuperMUC.

Figure 2: Mass fraction of molecular hydrogen as a function of time for four different disk surface densities.

Conclusions and Outlook

With a European team of experts from Garching, Heidelberg, Prague and Zurich we have performed the first simulations of the turbulent ISM in a galactic environment including the self-consistent formation of molecular Hydrogen and CO. In addition, the simulations include a model for energetic feedback during the late evolutionary phases of massive stars including SN explosions. With a first set of challenging simulations we have been able to demonstrate that gaseous galactic disks with properties similar to the Milky Way reach equilibrium molecular gas fractions of about 40 per cent. At higher surface densities the disks settle at higher molecular gas fractions, up to 90 per cent for 500 Msun/pc2, which is a typical value for massive galactic disks at high redshift. Moreover, we are able to formulate new constraints on the amount of CO-dark molecular gas and will be able to constrain gas outflow properties. We have already started the second phase of the SILCC project with simulations at higher spatial resolution, following the formation of small clusters of massive stars. We will resolve the relevant spatial scales to investigate the self-consistent driving of galactic winds. We will also study the important effects of magnetic fields on the turbulent multiphase structure of galactic disks.


[1] Bigiel, F., Leroy, A., Walter, F., Brinks, E., de Blok, W. J. G., Madore, B., Thornley, M. D. 2008, AJ, 136, 2846

[2] Brown, P. N., Byrne, G. D., Hindmarsh, A. C. 1989, SIAM J. Sci. Stat. Comput., 10, 1038

[3] Chabrier, G. 2002, ApJ, 567, 304

[4] Clark, P. C., Glover, S. C. O., Klessen, R. S. 2012, MNRAS, 420, 745

[5] Evans, II, N. J. and et al. 2009, ApJS, 181, 321

[6] Fryxell, B., et al. 2000, ApJS, 131, 273

[7] Glover, S. C. O., Mac Low, M.-M. 2011, MNRAS, 412, 337

[8] Górski, K. M., Hivon, E. 2011, Astrophysics Source Code Library, 7018

[9] Mac Low, M.-M., Klessen, R. S. 2004, Reviews of Modern Physics, 76, 125

[10] Marinacci, F., Fraternali, F., Nipoti, C., Binney, J., Ciotti, L., Londrillo, P. 2011, MNRAS, 415, 1534

[11] Moster, B. P., Somerville, R. S., Maulbetsch, C., van den Bosch, F. C., Macciò, A. V., Naab, T., Oser, L. 2010, ApJ, 710, 903

[12] Newman, S. F., et al. 2012a, ApJ, 752, 111

[13] Tacconi, L., et al. 2013, ApJ, 768, 74

• Stefanie Walch
• Thorsten Naab
• Philipp Girichidis
• Andrea Gatto
Max-Planck-Institut für Astrophysik Garching

• Simon Glover
• Ralf Klessen
• Paul Clark
• Christian Baczynski
Institut für Theoretische Astrophysik Heidelberg

• Richard Wünsch
Institute, Academy of Sciences of the Czech Republic, Prague

• Thomas Peters
Zurich University, Zurich

top  top