The formation of disk galaxies in computer simulations Lucio Mayer1,2, Fabio Governato3 & Tobias Kaufmann4 1-Institute for Theoretical Physics, University of Zurich, Winterthurestrasse 190, 8057 Zurich, Switzerland 8 2-ETH, Physics Department, Wolfgang Pauli Strasse 16, CH-8093 Zurich, Switzerland, 0 [email protected] 0 2 3-Department of Astronomy, University of Washington, Stevens Way, Seattle, WA 98193, USA n 4-Center for Cosmology, Department of Physics and Astronomy, University of California, Irvine, CA a 92697, USA J 5 Abstract 2 Theformation of diskgalaxies isoneof themostoutstandingproblemsin modernastro- ] physics and cosmology. We review the progress made by numerical simulations carried out h p on large parallel supercomputers. These simulations model the formation of disk galaxies - within the current structure formation paradigm in which the Universe is dominated by a o r cold dark matter component and a cosmological constant. We discuss how computer simu- t s lations have been an essential tool in advancing the field further over the last decade or so. a Recent progress stems from a combination of increased resolution and improved treatment [ of the astrophysical processes modeled in the simulations, such as the phenomenological 1 v description of the interstellar medium and of the process of star formation. We argue that 5 high mass and spatial resolution is a necessary condition in order to obtain large disks com- 4 parable with observed spiral galaxies avoiding spurious dissipation of angular momentum. 8 3 Arealistic modelof thestar formation history. gas-to-stars ratio andthemorphology of the . stellar and gaseous component is instead controlled by the phenomenological description of 1 0 the non-gravitational energy budget in the galaxy. This includes the energy injection by 8 supernovae explosions as well as by accreting supermassive black holes at scales below the 0 : resolution. Wecontinuebyshowingthatsimulations ofgas collapsewithincolddarkmatter v halos including a phenomenological description of supernovae blast-waves allow to obtain i X stellar disks with nearly exponential surface density profiles as those observed in real disk r galaxies, counteracting the tendency of gas collapsing in such halos to form cuspy baryonic a profiles. However, the ab-initio formation of a realistic rotationally supported disk galaxy with a pure exponential disk in a fully cosmological simulation is still an open problem. We argue that the suppression of bulge formation is related to the physics of galaxy for- mation during the merger of the most massive protogalactic lumps at high redshift, where the reionization of the Universe likely plays a key role. A sufficiently high resolution during this early phase of galaxy formation is also crucial to avoid artificial angular momentum loss and spuriousbulge formation. Finally, we discuss the role of mergers in disk formation, adiabatic halo contraction during the assembly of the disk, cold flows, thermal instability and other aspects of galaxy formation, focusing on their relevance to the puzzling origin of bulgeless galaxies. keywords:astrophysics, cosmology, computer science, fluid dynamics 1 Galaxy formation in a hierarchical Universe Galaxies occupy a special place in our quest for understanding the Universe. They are large islands in a nearly empty space and contain most of the ordinary baryonic matter, stars and 1 interstellar gas, that emits radiation and can thus be detected by astronomers . Galaxies come 2 essentially in two broad categories , those in which the luminous mass is arranged in a rotating disk of stars and gas, called disk galaxies or spiral galaxies because of the presence of spiral arms of gas and stars (Figure 1), and those in which the luminous mass is distributed in a smooth, featureless spheroidal structure with little or no rotation, also known as elliptical galaxies (Figure 1). Disk galaxies have a radial light distribution I(r) that is well fit by a decaying exponential law3, I(r) ∼ exp(−r/r ), where r is a characteristic scale length (r ∼ 2−4 kpc for typical spiral d d d galaxies). Indeed many disk galaxies contain also a spheroidal stellar component at their center, the stellar bulge, which has structural properties similar to an elliptical galaxies albeit being 2 much smaller in size . Both types of galaxies are known to contain dark matter, namely matter that is not traced by radiation. In disk galaxies dark matter clearly dominates over luminous matter by mass, as inferred from their high rotation speeds which requires the gravitational 4 pull of a massive and extended halo of dark matter . We live in a galaxy of the first kind, the Milky Way. Indeed disk galaxies are ubiquitous in the local Universe, and also at the largest distances and earliest epochs at which the best ground and space-based telescopes have been able 5 to study the morphology of galaxies reliably . Only the most massive galaxies in the Universe do not posess a disk component, while this becomes progressively more dominant compared to the spheroidal component as the mass of the galaxy decreases (Figure 2). 1.1 The theoretical framework The formation of disk galaxies is one of the major unsolved problems of modern astrophysics. The basic theoretical framework states that disk galaxies arise from the gravitational collapse 6 of a rotating protogalactic cloud of gas within the gravitational potential well of the dark halo . The gas cools via radiative processes during the collapse and eventually settles in centrifugal equilibrium at the center of the halo potential well forming a rotationally supported gas disk 7 provided that some angular momentum is retained during the collapse . These ideas were developed two decades agoandtheystillconstitute thebackbone ofdiskgalaxyformationmodels 8,9,10,11,12,13. What has changed dramatically since then is the cosmological context in which such idea is applied, which reflects the remarkable progress that cosmology has undergone in the meantime. After two decades of active debate there is now one cosmological paradigm according to which the energy density of the Universe is dominated by cold dark matter and a cosmological constant, while ordinary baryonic matter contributes only to a few percent level (an even smaller 14 contribution isyielded by neutrinos) . This modelis supportedby observations of thelargescale 15 mass distribution in the Universe traced by galaxies themselves and by the power spectrum 14 of density fluctuations inferred from the cosmic microwave background radiation . Cold dark matter interacts onlyvia gravity withitself andwithordinarybaryonic matter, andisnotsubject to any dissipative force. The governing evolutionary equation for dark matter is the collisionless 2 Boltzmann equation that describes a zero-pressure fluid, also termed a collisionless fluid . In this model, called ΛCDM (CDM stands for ”cold dark matter”, while Λ is the cosmological constant required to explain the observed acceleration of the Universe) structure forms hierarchically in a bottom-up fashion, starting from the amplification via gravitational instability of primordial small density fluctuations in the dark matter16,17. Because of the scale-free nature of gravity and the dissipationless nature of cold dark matter, in such a model one expects the formation of 18 self-similar, ellipsoidal collapsed objects, dark matter halos, at all scales . The largest halos are 17 the last to form . Also, direct three-dimensional simulations of structure formation in a CDM Universe predict that halos of any mass and size should contain a swarm of smaller halos, the so-called substructure19,20, as shown in Figure 3. The model predicts quantitatively the size and mass of the dark halo of a galaxy with a given measured rotational velocity (the rotational velocity of stars and gas probes the depth of the Figure 1: Top: a typical disk-dominated galaxy, the nearby spiral galaxy M102 in the Ursa Major constellation, 27 million light years from the Sun (the image was obtained by Chris & Dawn Schur from Payson, Arizona at 5150 feet elevation with an amateur telescope). A small spheroidal bulge is visible at the center of the disk. Bottom: a typical spheroidal galaxy with no disk component, the elliptcal galaxy M87, located at 60 million light years from us (credits in the picture). Figure 2: The ratio between the mass of the stellar spheroid (S) and the sum of the mass of the stellar spheroid and the stellar disk (T) as a function of galaxy mass from a galaxy sample of 1 the Sloan Digital Sky Surveys (SDSS) 2 (2007 Blackwell Publishing Ltd). Red points with error bars show the median S/T as a function of stellar mass together with the 10 and 90 percentiles of the distribution. galaxy gravitational potential well). For a galaxy like our own Milky Way, for example, its ∼ 12 observed rotational velocity of 220 km/s implies a halo mass of about 10 solar masses and 21 × 10 a halo radius of about 300 kpc (by comparison, the disk of our galaxy contains about 6 10 10 2 solar masses of stars and about 10 solar masses of gas ). Numerical simulations of the growth of dark matter halos also predict that the radial density profiles of such halos diverge near the center and are well described by a power-law, ρ ∼ r−γ (ρ being the density and r the spherically averagedradiusofthehalo), withγ = 1−1.5 22,23,24,25. Since darkmatterdominatesby massover ordinary baryonic matter, gas collapses within such halos, eventually forming a galaxy, because it is pulled inward by their gravitational attraction rather than collapsing due to its own gravity. In this scenario there is no such a thing as a galaxy forming in isolation, rather structure, both dark and baryonic, builds up via continous accretion andmerging of smaller systems containing a 26 mixture of dark and baryonic matter . This highly dynamical picture emerging from the current cosmology is the main difference compared to earlier attempts to study galaxy formation. Halos gain their angular momentum by tidal torques due to asymmetries in the distribution of matter, and also by acquiring the angular momentum originally stored in their relative orbit as they come together and merge into a larger system27,28,29. Prevailing models have then assumed that baryons and dark matter start with the same specific angular momentum before the collapse 8 begins since they are subject to the same tidal torques . 1.2 Computer simulations:the angular momentum problem The modern tool used to study the hierarchical growth of structure driven by gravity and the concurrent collapse of baryons within dark halos is represented by three-dimensional computer simulations that solve the gravitational and hydrodynamical forces between parcels of gas and dark matter. We will discuss the methodology employed by such simulations in the next section. For now it suffices to say that in the most popular simulation method both the gas (i.e. the baryonic component) and the dark matter are represented by particles so that structures are discretized in mass and space. The evolutionary equations, such as the collisionless Boltzmann equation for cold dark matter, the Euler equation for baryonic matter (baryonic matter is treated as an ideal gas) and the Poisson equation, which holds for both, are solved for such discrete rep- resentation of physical reality. Available methods to discretize physical variables and governing equations are constructed in such a way that they should converge to the exact continuum solu- tion for an infinite number of particles. As we will see in the next section, discretization itself, along with other aspects of the current methods, can introduce spurious effects in the computer models. Simulations take advantage of large parallel supercomputers in which hundreds of pro- cessing units are used simultaneously to compute the forces and advance the system to the next timestep. One important prediction of simulations of a CDM Universe is that halos have a rather universal value of the angular momentum (per unit mass) at any given epoch quite irrespective of their precise mass assembly history. This is parameterized via the dimensionless spin parameter λ = JE1/2/GM5/2, where E, M and J are the total energy, mass and angular momentum of the dark halo (G is the gravitational constant). One can show that λ is proportional to the ratio between the rotational kinetic energy and the kinetic energy in disordered motions associated with the halo. Halos have a universal distribution of spin parameters, which peaks at a value λ ∼ 0.035 29 . Simple spherical one-dimensional models that study disk formation in an isolated CDM halo (namely a halo that does not interact or merge with other halos) predict that the size of disks resulting from the infall and collapse of baryons matches the size of observed disks in galaxies 8 very well . The models use mainly two inputs, both coming from cosmological simulations, the halo density profile, which is related to the gravitational pull that drives the gas collapse, and the ∼ 14 Figure3: Exampleofamassivedarkhalo( 10 solarmasses)assembledbyhierarchicalmerging in a numerical simulation of the ΛCDM structure formation model (simulation performed by 44 L.Mayer and F.Governato with the GASOLINE code ). Many small halos (”substructure”) are orbiting inside the larger halo and will eventually merge with it as a result of dynamical friction eroding their orbital energy and orbital angular momentum. These dark matter lumps are the sites in which gas collapses and forms galaxies. initial specific angular momentum of gas as implied by the typical values of the spin parameter. They further assume that angular momentum is conserved during the collapse. Hence this result is simple and remarkable at the same time; it says that CDM halos have the right amount of angular momentum to form observed disk galaxies. For more than a decade researchers have tried to reproduce the latter result with fully three- dimensional computer simulations but have run into several problems. It was soon realized that, once the hypothesis of isolation is removed and hierarchical merging is accounted for, angular momentum can be lost by the gas to the dark matter due to a process known as dynamical friction30,2,31. During mergers, previously collapsed clumps of gas and dark matter fall into a larger dark halo and suffer a drag force as they move through the latter. The loss of angular momentum caused by the drag force, called ”dynamical friction”, is more effective when the 32 gas is distributed into cold and dense lumps rather than being smooth and extended . But gas is expected to be clumpy in a model with collisionless cold dark matter in which collapse can occur at all scales, and large halos grow by accreting smaller halos which bring their own 32 dense collapsed gas. As a result, early simulations were obtaining improbable small disks with ten times less angular momentum than real ones. Two types of solutions for this ”angular momentum problem” have been considered since then. The first is very drastic and calls for revising the cosmological model itself. Alternative models in which the dark matter has a non- negligible thermal velocity rather then being ”cold” would produce collapsed systems only above a characteristic scale because the thermal jittering will tend to smear out short-wavelength 33 density perturbations . These warm dark matter models (WDM) behave like CDM on large scale, thus maintaining its succesful features. The reduced clumpiness of dark matter halos in the WDM model implies that baryons are smoothly distributed rather than arranged in previously collapsed dense lumps when they fall into large galaxy-sized halos, and therefore lose less angular momentum by dynamical friction33,34. The second, less exotic possibility, is that baryons do not just follow the merging hierarchy imposed by dark matter but somehow decouple from it and remain much smoother. This could happen if the thermal energy content of baryons was enough to resist gravitational collapse, at least up to some critical mass scale. This way a fraction of the gas that would have entered a halo in dense clumps within smaller halos would instead enter with a smooth distribution, perhaps avoiding catastrophic loss of angular momentum. Various plausible astrophysical mechanisms can be responsible for increasing the thermal energy content of the baryons, for example the energy injection by supernovae explosions and the ambient radiation field produced by stars, accreting black holes or also external galaxies. There is, however, a third possibility. This is that the baryons clump excessively in computer simulations because the numerical methods adopted can introduce artificial loss of angular momentum. As we argue in this report, a solution lies probably in a combination of the latter two proposals, with no need of revising the standard cosmological structure formation model. The next two sections will be devoted, respectively, to the role of numerical effects in disk formation simulations, and to the modeling of gas thermodynamics and star formation in the simulations. We will then show how the structure of simulated disks is affected by different models of thermodynamics and star formation. Finally, we will summarize the current status of the field and the major problems that remain to be solved, including the puzzling origin of disk galaxies without a bulge. We will attempt to recall the most important contributions by the various groups actively involved in this field of research while at the same time covering in more detail some recent results of the research group to which the authors of this report belong. Figure 4: Disk size as a function of mass resolution in numerical simulations of disk formation in 48 an isolated CDM halo (2007 Blackwell Publishing Ltd). The three panels show density maps of gas in a slice through the centre of the gaseous galactic disk after 5 Gyr; the gas mass resolution decreases from the left to the right by about a factor 8 each snapshot (the maximum resolution is 1 million gas particles). The box side length is 20 kpc for every panel. At all resolutions disks show asymmetries such as central bar-like structures and spiral arms. The bar, however, is a strong and long-lived feature only for sufficiently high spatial resolution (set by the gravitational 48 softening) . 2 Modeling issues and numerical effects in computer sim- ulations of disk formation 2.1 Numerical methods Before discussing numerical effects in cosmological disk galaxy formation simulations we should keep in mind that until the time of writing this field of research has been dominated by one numerical method, smoothed particle hydrodynamics (SPH)35,36,37,38, in which particles are trac- ers of average fluid quantities such as density, velocity and temperature associated with a finite volume of the fluid. The fluid is thus discretized by means of particles and the calculation of the fluid-dynamical equations is carried out in a Lagrangian fashion. The Euler equation for an inviscid ideal gas is solved rather than the more general Navier-Stokes equation. An artificial viscosity term is introduced in the hydrodynamical force equation and in the internal energy equation to compensate for some artifacts resulting from the discrete representation of the fluid 39 equations, namely to avoid particle interpenetration and damp spurious oscillations in shocks . Although it is introduced for these good reasons, artificial viscosity also causes some unwanted numerical effects, suchasdampingtheangularmomentum ofthefluid. Oneofthereasonsbehind the dominance of the SPH method in this field is that it couples naturally with the most efficient and accurate methods to compute gravitational forces in the fluid as well as in the collisionless 37 dark matter component, the so called treecodes . Treecodes are an approximate but fast and 40 accurate way of solving the N-Body problem . In treecodes gravitational forces between all particles, dark matter and baryons, are solved via a type of multipole expansion of the gravi- tational field which reduces to individual particle interactions only at short distances. Eulerian techniques that solve gravity and the fluid equations on a fixed or adaptive grid have been less used in this field of research because they have been generally slower and less accurate when it comes to compute gravity. However, the scenario is changing rapidly with the appearance of several fast, parallel adaptive mesh refinement codes (AMR) that are beginning to have an impact in studies of galaxy formation within a cosmological context41,42,43. Restricting ourselves to the current particle-based methods, we will now briefly review the most important numerical artifacts that can severely affect the angular momentum content and structure of disks forming in the simulations. Many of the results that we will discuss in greater detail in this and in the 44 following sections were obtainedwith thetree+SPH codeGASOLINE insimulations performed on large parallel supercomputers. 2.2 Numerical effects: two-body heating One major problem of all particle-based simulations, both hydrodynamical and collisionless, is numerical two-body heating. Two-body heating is the spurious increase of kinetic or thermal energy of a particle representing gas, stars or dark matter due to a collision with another particle. Particles behave as gravitating point masses and therefore can undergo strong gravitational accelerations in close encounters with consequent large transfers of momentum. Such an effect is clearly an artifact of the discrete representation of a continuum by means of particles. For gas particles, such spurious large accelerations can be partially compensated by pressure or artificial viscosity that tend to deflect particles as they approach one another. In order to partially overcome this problem for all types of particles, gravity is decreased at small distances thanks 45 to the introduction of gravitational softening . Softening the gravity field at scales of order the interparticle separation is consistent with the fact that individual particles should represent a fairly large sub-volume of the hydrodynamical fluid or collisionless continuum (e.g. the dark mattercomponent) ratherthanrealpoint-masses(intypicalsimulationdarkmatterandbaryonic 6 − 7 particles can indeed weight 10 10 solar masses) However, in cosmological simulations of the CDM model dark matter particles are typically much heavier than gas and star particles because dark matter accounts for most of the mass. As a result massive dark matter particles can transfer significantkineticenergyduringgravitationalencounterswithotherparticlesdespitethepresence of gravitational softening. This spurious transfer of energy in two-body encounters increases the kinetic energy in random motions because any component of the velocity can be boosted as a result of a given encounter. Numerical experiments have shown that rotationally supported, thin stellar disks can be grad- ually degraded into a thick spheroidal distribution because the random velocity of the baryonic 46 particles becomes increasingly more important compared to ordered rotation . If two-body heating is moderate and a recognizable disk component survives the angular momentum of the disk along the original axis of rotation still decreases as a result of the randomization of velocity 46 vectors . Hence two-body heating induces artificialangular momentum loss. Gas particles suffer 47 an increase of temperature as a result of two-body heating . This affects the radiative cooling 47 efficiency of the gas because the cooling rate is a function of temperature . The increase of temperature occurs because the kinetic energy gained in two-body collisions is thermalized by artificial viscous dissipation. The only way to reduce these effects is to reduce the graininess of the mass distribution, which is only achieved by increasing the number of particles used in the simulation. This of course calls for more computing power. The reduction of spurious effects induced by two-body heating with increasing resolution has been tested systematically by means of toy-models that represent an isolated, already assembled galaxy34,46. This type of models is idealized but allows to gain insight into the phenomenology of the much more complex cosmological simulations where many interacting objects are simul- 5 taneously modeled. Such studies have indicated that > 10 particles are required in the dark matter of a single galaxy to keep the spurious kinetic energy increase to levels < 10%. In typical cosmological simulations many galaxies are followed at the same time and the latter becomes a tough resolution requirement. Figure 5: Top: growth of the disk mass (expressed relative to the total gas mass within the 48 dark halo) as a function of resolution (2007 Blackwell Publishing Ltd). Bottom: Evolution of the specific angular momentum in the disk as function of resolution48. The resolution of the simulations is increasing from LRLD (3 × 104 gas matter particles and 2.5 × 103 dark matter × 4 × 4 particles) to LR (3 10 dark matter and gas particles), IR (9 10 dark matter and gas × 5 6 particles) and HR (5 10 gas matter particles and 10 dark matter particles) ; IRLSNL differs 48 from the other runs in the prescription of artificial viscosity. See Table 1 in for details on the simulations.