3.1. The physics of the ISM
Despite the large amount of dark matter that they contain, disk galaxies are identified via their baryonic component, namely gas and stars. There is general consensus that the thermodynamics of the interstellar medium (ISM), the gaseous component in the disks of galaxies, is a crucial aspect of galaxy formation and evolution. Stars indeed form out of the ISM, being the end result of the gravitational collapse of the densest regions of clouds made of molecular gas. The interstellar medium in our Galaxy is multi-phase 62. A minimal ingredient of a realistic model of the ISM should account for the four main phases found in the disk of a galaxy, which, in order of increasing density, are the the warm ionized medium (WIM, ~ 0.5 atoms/cm3, T ~ 104 K), the warm neutral medium (WNM, ~ 0.5 atoms/cm3, T ~ 104 K) , the cold neutral medium (CNM, ~ 10-50 atoms/cm3, T ~ 103 K which comprises clouds of atomic hydrogen) and the cold molecular phase (mostly clouds of molecular hydrogen, > 300 atoms/cm3, T ~ 10 K ), with the first three phases being in approximate pressure equilibrium 62. A hotter, diffuse phase also surrounds the disk (hot intercloud medium, ~ 3 × 10-3 atoms/cm3, T ~ 105 K) and is likely constantly fed by supernovae explosions that inject large quantities of thermal energy and momentum in the ISM. The different phases are the result of thermal balance between radiative heating and cooling at different densities and, at the same time, thermal instability (perhaps coupled with gravitational instability) that determines the emergence of the two densest phases, the CNM and the molecular phase 62. Molecular gas is formed and destroyed via a number of microscopic interactions involving ions, atoms and catalysis on dust grains. These processes become biased towards the formation rather than towards the destruction of molecular hydrogen only at densities > 10 atoms/cm3. A great deal of energy in the interstellar medium is non-thermal; this turbulent energy, which is essentially observed as random gas motions of clouds and inside clouds is supersonic, being several times larger than the thermal energy at the scale of giant molecular complexes. Turbulent kinetic energy is thought to be the main agent that supports the largest molecular clouds (> 105 solar masses) against global collapse 63. The partial suppression of gravitational collapse owing to turbulent support also explains the low efficiency of star formation in our Galaxy (only a few percent of the molecular gas mass present in the Milky Way appears to be involved in forming stars). Magnetic fields also play an important role in resisting gravitational collapse at scales larger than 0.1 pc, while below this scale ambipolar diffusion and Ohmic dissipation give way yield to the action of gravity 62.
Supernovae explosions are a likely driver of ISM turbulence since the blast-waves generated by them can transport energy and momentum to scales as large as several hundred parsecs, perhaps up to kiloparsecs. This gives rise to dramatic outflows of gas above the disk plane in galaxies that are actively forming stars (Figure 10). Other drivers of turbulence in the ISM are probably operating, both at small scales, for example protostellar outflows, and at large scales, such as spiral waves generated by the large scale gravitational instability of the galactic disk at scales of 1 kpc and above. Magnetic fields probably also play a role since they can generate turbulence via magneto-hydrodynamical (MHD) waves 62. This brief summary highlights the complexity of the ISM physics. Things are complicated even further by the fact that the main mechanism of molecular cloud formation is still unclear, and so are the details of how molecular cores collapse into stars.
Figure 10. The disk galaxy NGC3079 which shows clear evidence of outflows triggered by supernovae explosions in the central region. The image was obtained with the Hubble Space Telescope (HST) (courtesy of NASA). The gas in the outflow is heated to large temperatures and shows up in X-rays. The central bubble is more than 3,000 light-years wide and rises 3,500 light-years above the galaxy's disk. The smaller photo at right is a close-up view of the bubble.
3.2. Modeling issues
How can we model the complex phenomenology of the ISM in large scale cosmological galaxy formation simulations? This has always been a major problem because the multi-phase structure of the ISM is also multi-scale, going from kpc scales, where most of the gas is ionized, to scales of tens of parsecs and below, where most of the gas is molecular, with the density changing by about 8 orders of magnitude across the different phases. Cosmological simulations follow the formation of galaxies in a volume of at least several tens of Mpc because, as we explained in section 2.4, the generation of angular momentum and other aspects of structure formation require that the large scale gravitational field is properly modeled. This calls for resolving seven orders of magnitude in length scales, and about ten orders of magnitude in density, going from tens of Mpc scales and densities of 10-7 atoms/cm3 to parsec scales and densities > 100 atoms/cm3. This can only get worse if we aim at following directly the process of star formation; indeed stars form from the collapse of the densest region of molecular clouds, the so called cores, at scales of < 0.1 pc 62. Currently, limitations in computer power as well as in the efficiency of available algorithms, make it possible to resolve at best scales of 100 pc in cosmological simulations (very recent simulations have been able to achieve a resolution of less than 50 pc, but they can only cover the first few billion years of evolution 43). Simulations with resolution adequate to follow directly molecular clouds, interstellar turbulence and star formation do exist, but are restricted to studying an isolated galaxy 64, 65 or a small region of a single galaxy 66. Detailed numerical models of the effects of supernovae explosions also exist, but again they are restricted to a small volume of the ISM 67, 68. For this reasons, the past decade or so has seen a lot of research activity being focused on designing the so called "sub-grid" models for simulations. These models essentially contain a phenomenological description of the processes occurring below the smallest scale resolved in the simulation. The phenomenological model is incorporated into the same parallel codes that compute gravity and hydrodynamics as well as radiative heating and cooling. Sub-grid models, being phenomenological, inevitably contain some free parameters that are tuned to reproduce important observables such as the typical star formation rate for a galaxy of a given mass, namely how much gas is turned into stars over a given amount of time.
The description of star formation is fully sub-grid, while the thermodynamics of the ISM is partially sub-grid. What do we mean by "partially" subgrid? One example is the following; radiative cooling is directly modeled for the range of densities accessible to the simulations while its effect below the grid is only implicitly accounted for in a phenomenological way. Since simulations typically lack resolution below a few hundred parsecs this sets a maximum density that the simulation can resolve of order 0.1 atoms/cm3, which is close to the density of the WNM (T ~ 104 K). For this reason cooling processes that are important below T ~ 104 K are usually neglected. Likewise, radiative heating is partially determined by the thermal and turbulent energy injection from supernovae explosions and/or from accreting massive black holes at the center of galaxies, both processes being not directly resolved, but also by the ultraviolet interstellar or intergalactic radiation background produced by stars, or by cosmic-ray and x-ray heating originating from various astrophysical sources, which can be directly included as constant heating terms in the internal energy equation without the need of a sub-grid model. Finally, interstellar turbulence cannot be resolved in galaxy formation simulations, nor it is accounted for in the sub-grid models. In what follows we will recall the main features of the sub-grid models widely used in galaxy formation simulations, pointing out their differences. First we will cover the star formation recipes and then sub-grid thermodynamics.
3.3. Star formation recipes
Star formation models used by different groups are very similar in essence. They describe the conversion of gas into stars 69, 70 A given star particle, once created, will represent a star cluster rather than a single star because of the limited mass resolution in the simulations. Sub-grid models in this context are all based on the idea that stars can form only when the gas climbs above some threshold density 69, 70. This density is usually taken to be 0.1 atoms/cm3. Indeed this is the typical density of the WNM, and is 6-7 orders of magnitude lower than the density of dense, star forming molecular cloud cores despite being a relatively high density in the context of simulations that model a large cosmological volume. Once the density of a gas parcel is above the latter threshold stars are formed over a timescale comparable to the local dynamical time, as expected from the fact that stars emerge from the gravitational collapse of dense gas 69, 70 The amount of stars formed over such a time span is regulated by an efficiency parameter that can be tuned to match observed star formation rates. This efficiency parameter crudely mimics the effect of unresolved physical processes that regulate the conversion of gas into stars. Among the latter, the formation and destruction rates of molecular hydrogen that determines how much gas is available in the star-forming dense molecular phase, and the effect of turbulence, magnetic fields, protostellar outflows and all those processes that determine how efficient is star formation within single clouds. Different algorithms then are distinguished by (1) how they implement additional conditions, based on e.g. the velocity field or temperature of the gas, as pre-requisites for star formation, and by (2) the details of the method used to convert a single gas particle into a collisionless star particle.
Despite various differences, all the sub-grid star formation recipes designed so far fare quite well at reproducing the slope of the observed correlation between the average star formation rate and gas density, the so called Schmidt-Kennicutt law 71. The latter relation is global in nature since it looks at the total amount of stars formed in an entire galaxy, and is one of the most fundamental observables that simulations use as a benchmark for their validity. In more detail, the Schmidt-Kennicutt law states that the density of cold neutral gas (CNM), mostly atomic hydrogen, correlates with the global star formation rate. This is a non trivial correlation since in reality it is only the molecular phase that is directly related to the production of stars. On the other end, somehow the molecular phase stems from the neutral atomic phase once this is able to achieve a high enough density. The success of simple sub-grid recipes in reproducing the Schmidt-Kennicutt law is thus probably related to the fact that they are all based on a threshold gas density and a dynamical timescale 72. Both numbers are determined by how the gas density evolves with time; this in turn is likely controlled by the large (kpc) scale gravitational instability in galactic disks, which is resolved in the simulations.
Nevertheless, for low mass galaxies large deviations from the Kennicutt relation occur, and such deviations are also evident in single star forming regions of well studied nearby galaxies 73. In order to reproduce the latter, more complex observational scenario high resolution simulations of galaxies have begun to include a sub-grid description of the molecular phase starting from the CNM instead of having to bridge all over from the WIM and WNM to stars 74, 75. These first attempts show that a model that incorporates the formation and destruction of molecular hydrogen as a function of density and ambient temperature allows to reproduce detailed correlations between the local star formation rate and the local density of molecular gas as traced by various molecular tracers in observations of individual galaxies 73. These recent models produce realistic results for the star formation properties of galaxies of varying masses 76, 75, despite the fact that heating by supernovae explosions is neglected. They provide a way to follow the molecular gas fraction when the resolution is not high enough to model molecular gas formation directly by using the gas at CNM densities as a proxy for the amount of molecular gas 75, 77. For this procedure to be sensible the number of gas particles must of course be high enough to allow resolving the density of the CNM (> 50 atoms/cm3). The latest cosmological simulations of galaxy formation are just now starting to approach this regime.
3.4. Modeling ISM thermodynamics with supernovae feedback
Thermal and kinetic feedback Sub-grid models of supernovae feedback come in a wide variety. Early attempts in the 90s were essentially based on two types of implementations. In the first implementation 69, 34 a fraction of the energy of individual supernovae explosions (1051 erg) is damped to the surrounding gas as thermal energy. In the second method a fraction of the energy of the explosion is converted into kinetic energy of the surrounding gas particles rather than into thermal energy. The latter model is motivated by the fact that the blast-wave produced as a result of the stellar explosion will not thermalize immediately but will expand in the interstellar medium for some time as a result of its bulk motion 78, 79, 80. With thermal feedback the gas loses quickly the added energy because the radiative cooling time is very short for the typical temperatures reached by the gas (T ~ 105 K). Eventually the kinetic energy added to the gas in kinetic feedback is also thermalized and radiated away, but on a longer timescale. As a result, kinetic feedback is more effective at increasing the internal energy of the gas compared to thermal feedback; gas cooling is counteracted more efficiently, and in cosmological simulations less gas ends up in dense lumps. Overall the gas component loses less angular momentum while merging, resulting ultimately in larger disks 79, 80, 81. However, this method is strongly resolution-dependent and is not directly informed by the actual dynamics and thermodynamics of the supernova remnant since (the magnitude of the velocity kick given to the particles is arbitrary). Finally, both the thermal and the kinetic method do not account for the multi-phase structure of the ISM.
The rate of supernovae heating depends on how many supernovae type I and type II go off in a given time, which in turn depends on the adopted initial stellar mass function (IMF). The IMF, N(m*), quantifies how many stars of a given mass m* are present in a representative ensemble of stars. Each star particle in the simulation is indeed not a single star but rather represents a star cluster whose unresolved member stars obey the chosen IMF. Most sub-grid models assume a standard stellar mass function based on the observed mass function in galactic star clusters. Standard mass functions are well described by power-laws, N(m*) ~ m*-, < 2, and are dominated by stars comparable or less massive than the Sun (for a review see 82). Recently, some researchers have begun to distinguish between a late epoch, close to present day, in which the IMF is one of the standard forms, and an early epoch, ten billions of years ago or more, in which IMF was likely dominated by stars with masses well exceeding 10 solar masses, called a top-heavy IMF 57, 60. Of course in the latter case the amount of energy damped by supernovae is much larger because a higher number of massive stars, those that rapidly explode into supernovae type II, is produced for a given star formation episode. The stronger supernovae heating rate resulting from a top-heavy IMF has a positive effect on suppressing the formation of cold and dense gas clumps in dark halos 57, 60.
Adiabatic feedback and blast-wave feedback In a third, alternative approach to feedback the energy of the explosions is damped to the gas as thermal but the radiative cooling of such gas is stopped for a timescale of about 20-30 Myrs, based on the assumption that it will take this much time for the (unresolved) interstellar turbulence generated by supernovae explosions to be dissipated 83, 84. Of course this is a very simplistic way of accounting for the non-thermal energy budget of the ISM. In fact it is known that turbulence does not just behave as an effective pressure support for the ISM but can also drive cloud collapse, thus aiding star formation, by creating local density enhancements 85. When applied to cosmological simulations, this "adiabatic" feedback proved even more effective than kinetic feedback in producing an extended disk component 84. Recently, we have developed further this idea 86; in the remainder we will refer to our model as the "blast-wave" model. In the latter model the timescale during which cooling is shut off is self-consistently calculated based on a sub-grid model of the blast-wave produced by a supernova explosion. By temporarily preventing the cooling of the hot phase created by supernovae feedback this type of methods naturally produces a two-phase medium with hot bubbles triggered by supernovae explosions (T > 105 K) surrounded by a colder, filamentary phase (T ~ 104 K) (see Figure 11). Based on the density and temperature range, these two phases roughly represent the hot intercloud medium and a combination of the WNM and WIM.
Figure 11. Gaseous distribution of an isolated galaxy simulated with blast-wave feedback 86 (2006 Blackwell Publishing Ltd). Note the holes punched by the supernovae explosions (hot bubbles) surrounded by much colder and dense gas in filamentary structures.
The blast-wave model stems from a classic model developed in the late seventies 87 which reproduces the main features of the multi-phase ISM in our Galaxy. In the latter the blast-wave sweeping the ISM undergoes an adiabatic expansion phase during which radiative cooling is negligible. The radius of the blast-wave as a function of time can be directly computed from the local physical parameters of the ISM. In the numerical implementation such radius defines the size of the volume of gas particles that are unable to cool during the adiabatic phase. The adiabatic phase lasts of order 30 Myrs, after which radiative losses would become efficient and the gas is again allowed to cool radiatively. The blast-wave can affect volumes with length scales of up to a kiloparsec, as observed in star supernova-driven outflows (Figure 10); the latter are well resolved in the highest resolution cosmological simulations available today 55. The delay introduced in the cooling time has also a direct effect on star formation since the gas becomes hotter than the typical temperature of gas eligible to form stars. The latter is a free parameter in star formation recipes; it is typically set close to 104 K, the lower limit of the radiative cooling function for atomic gas comprising only hydrogen and helium (the standard composition adopted in cosmological simulations). This somewhat mimics the fact that the interstellar turbulence fed by supernovae explosions might temporarily suppress the collapse of dense molecular cores into stars. The blast-wave feedback model has only one free parameter, the efficiency of supernovae feedback, namely what fraction of the energy generated by the supernovae explosions is damped to the gas 86. An individual particle belongs to one phase only and there is no enforcement of pressure equilibrium between different phases. Indeed in the interstellar medium pressure equilibrium applies only to the warm and cold diffuse atomic phase but not to the star forming cold molecular phase or to the hot intercloud medium produced by supernovae explosions 62.
ISM Models with effective equation of state Another method is based on treating the interstellar medium via an effective equation of state that accounts for the relative contributions of a hot and a cold phase in a statistical fashion 88. This method is the modern version of earlier attempts made to model interstellar gas as a two-phase medium with a hot and a cold phase 89. Such earlier methods were splitting gas particles in two sets, cold and hot particles, with the transition from one to the other phase being regulated via radiative cooling and heating by supernovae feedback. Some of these methods were decoupling the hot and cold phase by solving the hydrodynamical equations separately for two sets of particles 51, 90. In the new method gas particles come only in one species; each gas particle in the simulation is effectively representing a finite volume of the ISM with its share of cold and hot phase being determined by the local density 88. The share of hot and cold phase determines how compressible is the gas, namely the form of the function P() (P is the pressure and the density of the ISM). The larger the fraction of hot phase associated with a given particle the stiffer its equation of state, and therefore the higher will be the gas pressure assigned to it for a given density. The hot phase, and thus the pressurization of the gas, is the result of supernova feedback and its effects on the interstellar medium. These effects are treated in a phenomenological way and involve the creation and destruction rates of unresolved star forming molecular clouds. The phenomenological equations contain free parameters. In denser regions the equation of state becomes increasingly stiffer as star formation generates heating via supernovae explosions and augments the fraction of hot phase. When the hot phase increases, however, cooling also becomes more efficient and tends to replenish the cold phase. For reasonable values of the free parameters the simulated galaxies tend to reach quickly self-regulation, in the sense that the equation of state approaches a steady-state regime. A more sophisticated version of this model also includes the additional heating resulting from accretion onto massive black holes at the center of galaxies 91, which typically renders the equation of state stiffer. On the contrary, switching off heating by supernovae feedback produces a softer equation of state that approaches the isothermal regime (as expected based on the efficient cooling rate at the density and temperatures typical of galactic disks). Since the description of the interstellar medium given by the effective equation of state is statistical in nature the two phases are implicitly assumed to be in pressure equilibrium locally, which is not the case in the blast-wave feedback model. This is a major difference between the models, and its consequences should be systematically explored in the future. This method produces very smooth disks in the simulations when feedback is included 56 (without feedback the isothermal disks undergo rapid fragmentation owing to gravitational instability) as opposed to disks rich of filamentary structures and flocculent spiral arms arising in the blast-wave model (Figure 11).
Recently, the effective equation of state model was modified to allow for strong supernovae feedback due to a top-heavy stellar mass function arising in specific conditions 57, such as during galaxy mergers (there are indications that the high star formation rates produced by mergers might lead to a top-heavy stellar mass function by changing the balance between heating and cooling that controls the collapse of molecular cores 92). This produces a "burst mode" in which heating is so efficient that it can lead to partial evaporation of the ISM in low mass galaxies. This appears to reduce the mass of the bulge and make a galaxy more disk dominated, although galaxies always end up with some bulge component.