We hope that the Illustris simulation will be of interest to scientists of all disciplines as well as to the broader non-scientific community. Below we describe the 'what', 'why', and 'how' of the project at a level accessible to those without any expertise in physics, astronomy, or computer science.

Motivation & Big Ideas

The standard model of cosmology posits that the mass-energy density of the Universe is dominated by unknown forms of dark matter and dark energy. Testing this extraordinary scenario requires precise predictions for the formation of structure in the visible matter, which is directly observable as stars, diffuse gas, and accreting black holes. These components of the visible matter are organized in a 'Cosmic Web' of sheets, filaments, and voids, inside which the basic units of cosmic structure - galaxies - are embedded. To test our current ideas on the formation and evolution of galaxies, we strive to create simulated galaxies as detailed and realistic as possible, and compare them to galaxies observed in the real universe. By probing our successes and failures, we can further enhance our understanding of the galaxy formation process, and thereby perhaps realize something fundamental about the world in which we live.

The Illustris project is a set of large-scale cosmological simulations, including the most ambitious simulation of galaxy formation yet performed. The calculation tracks the expansion of the universe, the gravitational pull of matter onto itself, the motion or "hydrodynamics" of cosmic gas, as well as the formation of stars and black holes. These physical components and processes are all modeled starting from initial conditions resembling the very young universe 300,000 years after the Big Bang and until the present day, spanning over 13.8 billion years of cosmic evolution. The simulated volume contains tens of thousands of galaxies captured in high-detail, covering a wide range of masses, rates of star formation, shapes, sizes, and with properties that agree well with the galaxy population observed in the real universe. We are currently working to make detailed comparisons of our simulation box to these observed galaxy populations, and some exciting promising results have already been published.

What are cosmological hydrodynamical simulations and why are they useful?

The Lambda Cold Dark Matter (Lambda-CDM) paradigm of cosmology, currently favored by observations of the large-scale distribution of galaxies in space, implies that the cosmos is filled with three distinct components: normal matter (which astronomers term 'baryons'), dark matter, and dark energy. The mathematical models which govern the physical behavior of these components are sufficiently complex that they can only be solved exactly for very particular, simplified "test" problems. Understanding how the nearly uniform, primordial universe evolved into the many diverse phenomena we observe in the night sky today therefore requires the use of computer simulations, which can numerically evolve a representation of some fraction of the universe forward in time.

Simulations of the combined evolution of dark matter and dark energy, which only include the force of gravity, have been run to great success over the past few decades. Recently, such simulations have reached staggering scale, including on the order of 1 trillion particles, each of which exerts a gravitational force on every other. However, such "DM-only" simulations cannot predict the distribution of galaxies made up of normal matter, severely limiting their utility as a means to directly connect with the observations. Our approach for establishing this link is through directly accounting for the baryonic component (gas, stars, supermassive black holes, etc.) in cosmological simulations that calculate fluid motion ("hydrodynamics") as well as gravity, in principle offering a self-consistent and fully predictive methodology. For example, the image below shows two galaxies (out of thousands of similar systems), one on each row, evolving in time from left to right, from when the universe was a quarter its current age, to the present. The top galaxy shows the massive, red, elliptical-shaped galaxy forming after a series of mergers with other systems, whereas the bottom galaxy reveals the formation of a smaller, bluer, disk-shaped galaxy forming after a less violent history of interactions.

However, the vast computational challenges related to simulations of the baryonic component have thus far precluded wide-spread adoption of this approach. Large-scale predictions of the galaxy population have thus, in the past, mainly been obtained with a different method, whereby galaxies are assigned to live in specific dark matter objects after the simulation has already been completed, and given properties and characteristics based on relatively simple descriptions of the relevant physics. In the Illustris simulation we pursue the computational approach, with a large-scale, cosmological volume, modeling directly the dark and baryonic components of matter.

Project Description

What physics does the simulation include?

In addition to gravity and hydrodynamics, very complex physical processes such as chemical processes in the diffuse gas, radiation, and magnetic fields also affect cosmic structure formation. Moreover, structure formation is a self-regulating process, in the sense that the structures that form, in particular stars and black holes, affect their surroundings and the subsequent evolution of the next generation of structures. In Illustris, a comprehensive (even if not complete) set of physical processes such as star-formation-driven galactic winds, and black hole thermal energy injection, are modeled throughout cosmic history. These sophisticated models are crucial for achieving a realistic population of modeled galaxies.

What computational methods does the simulation leverage?

Over the past few decades, computer simulations of the evolution of the universe, including only the effects of gravity, have reached a point of maturity and accuracy, the Millennium simulation standing as a prime example. Those which attempt to also include a treatment for gas, such as Illustris, have proven to be significantly harder. A number of fundamentally different methods exist for simulating gas on a computer. In astrophysics, most researchers have used one of two approaches: (i) “smoothed particle hydrodynamics”, or SPH, where the mass of the gaseous fluid is parceled out to a discrete number of particles. These particles move in response to the combined forces of gravity and hydrodynamics, and their position at any time indicates where the gas is, and how it is moving. (ii) The second approach of “Eulerian” or “mesh-based” methods, typically utilizing a scheme called “adaptive mesh refinement” or AMR. In this method, space itself is divided up into a grid, and the flow of gas between neighboring cells of this grid is computed over time.

The Illustris simulation uses a different approach, as implemented in the AREPO code, which we typically refer to as employing a “moving, unstructured mesh”. Like in AMR, the volume of space is discretized into many individual cells, but as in SPH, these cells move with time, adapting to the flow of gas in their vicinity. As a result, the mesh itself, called a Voronoi tessellation of space, has no preferred directions or grid-like structure. Over the past few years we have shown that this new type of approach for simulating gas has significant advantages over the other two methods, particularly for large, cosmological simulations like Illustris. In addition to being accurate, the AREPO code is also efficient – it can run on tens of thousands of computer cores simultaneously, leveraging some of the largest computers that currently exist for scientific research within the “high performance computing” (HPC) community. The Illustris simulations were run on supercomputers in France, Germany, and the US. The largest was run on 8,192 compute cores, and took 19 million CPU hours (the equivalent of one computer CPU running for 19 million hours, or about 2,000 years).

What are we learning?

Before leveraging the output of the simulation to better understand the different types of galaxies that form, and the evolutionary paths they take towards their present day state, we need to ensure that the simulation is a faithful representation of reality. Due to the complexity of the problem, numerical calculations can never fully capture the answer across all scales of space and time. A finite resolution (the size of the smallest details which we can include) means that some processes, such as the birth of individual stars, cannot be directly captured in a cosmological simulation. As a result, we implement many physical approximations, and the burden then arises to ensure that the particular problems of galaxy formation we wish to study are well represented under these approximations. Therefore, we introduce some of the first insights from the Illustris simulation, by first describing this essential comparison against observations:

  • We successfully reproduce a wide range of observable properties of galaxies and the relationships between these properties. Our model is specifically tuned to match two quantities: the present day ratio of the amount of stars to dark matter (DM), for galaxies of all masses, and the total amount of star formation in the universe as a function of time. These two relationships are well constrained by observations, and by construction, we have good agreement with both. Going beyond, we find that the simulation also matches the galaxy stellar mass and luminosity functions (i.e., how many galaxies of each mass and brightness exist), at the present day as well as back in time. The "specific star formation rates" (the rate of new stars being formed in a galaxy, divided by the amount of already-existing stars) also agrees well with observations, not only at the present epoch, but also for all time back to when the universe was thirteen billion years younger than it is today.
  • We precisely measure the gas content of the universe, and where it resides. Galaxies in the simulation have the appropriate amounts of neutral hydrogen gas, molecular gas, and "metals" (all elements other than hydrogen and helium), where observations exist. Where data does not yet exist, we can make specific predictions. For instance, we measure that the molecular gas fractions of galaxies increase going back in time until a specific point, when the age of the universe was 2.7 billion years, at which point the increase halts, and its relationship to the amount of stars within the same galaxy becomes significantly stronger. Outside of individual galaxies, Illustris also predicts that at the present time, the majority of gas (~81%) remains in the "intergalactic medium" (the space between galaxies), but that this gas contains only a minority (~34%) of the metals so-far produced in the universe. We encourage you to investigate the relationship between the gas, metals, stars, and dark matter at one instant in time in The Explorer.
  • We investigate the number of "satellite" galaxies, their properties, and their connection to cosmology. Satellite galaxies are those which orbit in the gravitational influence of a much larger, nearby galaxy. Such associates can be relatively small "groups" like our own Local Group, which contains a few dozen galaxies, or larger "clusters", which can contain thousands or tens of thousands of members. We measure the distributions of satellites as a function of distance to their cluster center, for the 10 most massive clusters in our simulation. This profile has a substantially different shape than dark matter only results - more satellites survive in the central region, because the inclusion of stars and gas makes them more resistant to being torn apart by strong tidal forces. We predict that the shape and amplitude of the profile measuring the number of satellites has little evolution further back in time than it has currently been measured. Illustris also demonstrates the process by which environment quenches star formation, i.e. the fact that galaxy star formation 'shuts off' when in the close vicinity of a massive cluster, when compared to an otherwise identical galaxy far from any massive neighbors.
  • We study changes in internal structure as galaxy populations evolve in time. A primary measure of galactic structure is the "circular velocity" or "rotation curve", which measures the speed required for a star in a galaxy to maintain a circular orbit (against the pull of gravity) at a given distance from the center. Observationally, the circular velocity profiles of galaxies were found to be nearly constant out to very large distances, indicating that there must also be a significant amount of mass at similarly large distances. As the number of stars drops rapidly, this is one of the key arguments in support of dark matter on the scale of galaxies. In Illustris we generally find such flat, or gradually rising, rotation curves within the distance where gas and stars are a larger contribution to the gravitational pull than DM. We find that the maximum velocity reached is larger when gas and stars are included in the simulation. However, we find that when the universe was about half its current age, the stellar component of galaxies was substantially more concentrated, causing an even larger maximum circular velocity, at even smaller distance from the center. Finally, we also study particular populations from when the age of the universe was a quarter of its current age (coinciding with the peak rate of star formation). Investigating how they grow, we find that galaxies with star formation rates higher than average at this earlier time are substantially more likely to undergo a major merger event (i.e. a collision with another galaxy of similar mass).
  • The impact of gas on the structure of dark matter. Some of the largest planned observational campaigns of the near future will attempt to measure the set of "cosmological parameters" at an extremely precise level, by mapping the distribution of galaxies throughout the universe. These parameters constrain the currently accepted "Lambda-CDM" cosmological model, and have implications for e.g. the accelerating expansion of the universe, the age of the universe, the amount of matter in the universe and its division among components such as gas, dark matter, and dark energy. However, expectations for these experiments have so far been generated based on dark matter only simulations. It is already well known that the presence of gas can affect dark matter properties on different scales, for instance causing the centers of dark matter halos to contract to higher densities. We measure the impact of including gas on a function called the "matter power spectrum", and show that the change is up to ~40% when compared to DM-only simulations, and that this change is a complex function of the scale of space you consider. Fully understanding this and related effects of gas on dark matter is an absolute requirement for future cosmological experiments.
  • Producing "mock" observations. An important step in bridging the gap between theory and observation is the ability to perform any analysis or comparison in exactly the same manner. For example, in the real universe we observe the light emitted from stars at various wavelengths, but in the simulation we calculate the physical properties of the stars such as mass, age, and chemical composition. Taking these intrinsic properties of a star and calculating how it would look like if observed is a more straightforward process, with less ambiguity, than going in the opposite direction. We create synthetic or "mock" observations of the starlight from galaxies in this manner. First, we re-create one of the most iconic images in astronomy, the Hubble Space Telescope "Ultra Deep Field" (the image below is split in half, to the left and right - one half is real, and one is simulated, can you tell which?), which is the first time a hydrodynamic simulation could construct a faithful deep UDF-like image. To begin to quantify how our galaxies would look as observed with real telescopes, we have also begun to "image" individual objects, which are available to browse in the Galaxy Observatory. This will allow us to use exactly the same tools as observers, and answer questions such as "which statistic, X or Y, is a more reliable indicator of a galaxy currently undergoing a merger or interaction with another system?"

Many areas of the simulation remain unexplored as researchers associated with the Illustris project continue to investigate interesting and unexpected results. The ultimate goal in each case is a deepening of our understanding of the processes by which the observed population of galaxies formed and evolved over cosmic time. Understanding galaxy formation, however, has broader implications for scales both smaller and larger: galaxies are the sites where individual stars, as well as the planets and the possible tracers of life they harbor, form. At the same time, galaxies are the visible tracers of the large-scale structure of the universe itself, and therefore inexorably linked to the Big Bang, inflation, and the mysterious nature of both dark matter and dark energy.

Moving forward, particular attention to specific discrepancies between simulation and observation will allow us to better understand our galaxy formation model. By addressing these issues, and improving in both physical realism and numerical methods, we can start to think about the next generation of cosmological hydrodynamic simulations after Illustris.

We describe the key science goals of the project, relate them to open questions in our theoretical understanding of galaxy formation and evolution, and describe more technical details of the simulation below.


Since the advent of modern observational surveys, the sheer volume of available data on the properties of galaxies, their distribution within the large-scale structure of the universe, and their evolution in time has exploded. Surveys including SDSS, 2DF, DEEP2, and CANDELS, and upcoming projects such as LSST, have provided an increasingly precise observational constraint against which any theoretical idea of galaxy formation must be benchmarked. The LCDM cosmology favored by early-universe CMB experiments – particularly WMAP and now PLANCK – provides an extraordinarily precise measurement of the initial conditions for the problem of cosmic structure formation. Numerical calculations are required to probe past the linear regime of early time, and efforts modeling only the gravitational interaction of dark matter and dark energy, neglecting the role of gas and baryonic physics, have led to significant physical insight. Nevertheless, such DM-only simulations do not directly predict anything about the galaxies themselves, requiring an extra step in order to bridge the gap with observations.

Over the past two decades two dominant approaches have been used to establish this link: (1) the technique of ‘semi-analytical modeling’, whereby baryonic physics are modeled at the scale of an entire galaxy, and applied in post-processing on top of DM simulations, and (2) hydrodynamic simulations, whereby the evolution of the gaseous component of the universe is treated using the methods of computational fluid dynamics. The latter approach enables the complex interaction of the different baryonic components (gas, stars, black holes) to be treated at a much smaller scale, ideally yielding a self-consistent and powerfully predictive calculation.

In the context of previous large simulations

Hydrodynamical cosmological simulations, due to their high computational cost, have usually targeted specific problems, including the use of ‘zoom-in’ simulations which resolve only one or a few galaxies. They have been harnessed to study the effects of different models or model variations on a particular problem, without explicitly aiming to reproduce the large number of observational constraints available. Only in the past few years have several groups started projects with similar approaches and objectives as the Illustris simulation. These simulations have been steadily improving in three distinct ways: (1) by increasing in size, both in terms of volume and the number of resolution elements employed, (2) by improving the scope, complexity, and physical fidelity of the sub-grid models required to provide a complete and accurate description of the many processes that govern galaxy formation, and (3) by developing more accurate and efficient numerical methods. This trend can be seen in the figure, which shows simulations similar to Illustris (periodic volumes evolved down to z=0) in terms of the number of resolution elements used, as a function of time (Illustris is #19).

Nearly tracking Moore’s Law, the growth in simulation size is exponential with a doubling time of approximately 20 months. It is only now, with the available computing power, the sophistication of the numerical approach, and the fidelity of the physical models, that we can simulate a statistically significant volume of the universe down to z=0 with sufficient detail to resolve the internal structure of individual galaxies.

Key Science Areas

Below we detail three areas of investigation which have strongly motivated Illustris. However, the project also has the potential to provide important contributions, both in theoretical understanding and in the interpretation of observational data, in a large number of subfields within astrophysics: (1) the low-density IGM, and Lyman-alpha forest, (2) high column density absorbers, LLSs/DLAs, (3) disk galaxies, (4) spheroids, (5) galaxy populations, (6) the CGM, (7) galaxy interactions, (8) galaxy mergers, (9) super massive black holes, (10) starbursts and star-formation modes, (11) AGN and quasars, (12) satellite galaxies, (13) hot halo gas, X-ray, (14) groups and clusters, (15) observed high redshift phenomena, (16) background radiation fields, (17) large scale structure, (18) impact of baryons on dark matter, (19) reionization, and (20) gravitational lensing.

Structure, kinematics and morphology of galaxies

A long-standing problem in galaxy formation simulations has been that simulated galaxies are overly compact compared to the observed population, indicating that the angular momentum distribution between the galaxies and their surroundings was incorrect. Our preliminary findings using the AREPO code indicate that the moving-mesh approach may help solve this problem. Early simulations with AREPO yielded galaxies that were more disk-like, more spatially extended, and more rotationally supported than those in otherwise identical simulations performed with the SPH code GADGET-3. However, for the sake of comparison, these simulations did not include feedback processes that are believed to be responsible for removing a large fraction of baryons from galaxies. These were also relatively small periodic volumes, 25 Mpc/h in extent per dimension.

The Illustris simulation volume is 27 times larger, allowing us to obtain a statistical representation of the cosmological galaxy population. By comparing the morphological distribution of this population to locally observed galaxies we can calibrate uncertainties in our modeling of star formation and feedback. We can study the evolution of the gas and metal content of galaxies over time, and compare against the observed galaxy mass-metallicity relation. We can investigate the properties and nature of damped Lyman-alpha absorbers. Unlike in DM-only simulations, we can make a prediction for the galaxy-galaxy merger rate as a function of mass ratio, redshift, and mass. These results are important in understanding the growth of galaxies through mergers, and for calibrating the prescriptions used to quantify merger rates in semi-analytical models. Our previous simulations have already demonstrated that exponential disks are a natural consequence of the hierarchical galaxy assembly process. We can now use Illustris to understand in detail how these profiles arise, by tracing the evolutionary history of the simulated galaxies back to higher redshift, and following the acquisition of angular momentum in the baryon component.

Gas flows in galactic halos and the physical state of the CGM

The accretion of gas is one of the main drivers of the internal growth and evolution of galaxies, yet the process by which gas from the IGM interacts with the halo/CGM regime and finally makes it into the galactic ISM remains uncertain. Our numerical method has demonstrated significant advantages in the hydrodynamic treatment, particularly important in this regime. Additionally, the Monte Carlo tracer particle scheme and the inclusion of many such tracers in the main Illustris simulation allow us to reconstruct the time evolution of gas which is accreting into halos and galaxies.

Observations have revealed that galaxies drive high velocity outflows into their surroundings. This material can shock heat, adding hot gas to the halo, while denser parts of the outflows can drive small shocks and possibly survive to large distances. Outflows can compress halo gas, increasing its density and radiative cooling rate, which can cause it to 'rain back' onto the galaxy. Outflows will also interact with cosmological inflow, possibly modulating the nature or rates of this supply of new gas for the central galaxy. We can address the consequences and detectability of both the infalling and outflowing material in galactic halos, and compare to observational signatures of gas in halos and metals in the CGM as probed by methods such as quasar absorption line spectroscopy. By accounting for chemical enrichment processes - and the energetic feedback which redistributes both gas and metals - we can make accurate synthetic observations of the circumgalactic material in the simulations.

Co-evolution of galaxies and their central black holes

In the past few decades, a compelling picture has emerged indicating that black holes are common inhabitants of the vast majority of galaxies. The observed scaling relationships between BH masses and their host galaxy properties imply that we need to understand how BHs grow and affect their surroundings. We can explore the evolution of the black hole accretion rate as a function of time, and the link to the host galaxy. Current observational efforts seek to constrain evolution in the scaling relations at z>0, which we can extract directly from the simulations, as well as how the scaling relations depend on e.g. the morphological type of the host galaxy. We can also explore the impact of AGN feedback on the dichotomy between blue and red galaxies, the link to the quenching of star-formation, and how this picture changes with redshift.

It is commonly assumed that active SMBHs at high redshift reside in rare, massive halos embedded in gas-rich, over-dense environments, while recent observational efforts are still working to unambiguously demonstrate that this is the case. Although our simulated volume is not sufficiently large to capture the formation of the brightest quasars at z~6, the volume can be used for determining, and with statistical robustness, the typical environments where less extreme quasars form. Their evolution down to z=0 will allow us to address several open questions. For instance, in analogy to the downsizing of galaxies, it is observed that at high redshifts more massive BHs are powering bright AGN, while at lower redshifts the peak of quasar activity shifts towards smaller mass BHs. With a representative sample of simulated BHs, whose accretion rates are measured directly from the surrounding gas, it will be possible to probe the downsizing of the AGN population as a whole, and why and in which hosts SMBHs enter radiatively inefficient regimes.

Simulation Details

We follow the coupled dynamics of DM and gas with the robust, accurate, and efficient quasi-Lagrangian code AREPO. In this approach, an unstructured Voronoi tessellation of the simulation volume allows for dynamic and adaptive spatial discretization, where a set of mesh generating points are moved along with the gas flow. This mesh is used to solve the equations of ideal hydrodynamics using a second order, finite volume, directionally un-split Godunov-type scheme, with an exact Riemann solver. The gravitational force is calculated with a split Tree-PM approach, where long-range forces are calculated from a particle-mesh method, and short-range forces are calculated with a hierarchical octree algorithm. Our galaxy formation model is based on the inclusion of several additional astrophysical processes:

  • Gas cooling and photo-ionization: the cooling function is calculated as a function of gas density, temperature, metallicity, UV radiation field, and AGN radiation field. The UV background is a spatially uniform, time dependent field for which reionization completes at z~6.
  • Star formation and ISM model: our subgrid model for the ISM puts dense gas above 0.13 particles per cubic cm on an effective equation of state, assuming a two-phase medium of cold clouds embedded in a tenuous, hot phase. Star formation occurs stochastically following the KS-law and with a Chabrier IMF.
  • Stellar evolution: stellar populations return mass to the gas phase through stellar winds and supernova. We track and follow the rates associated with SN Ia, SN II, and AGB stars. The evolving abundances of nine elements (H, He, C, N, O, Ne, Mg, Si, Fe) are independently considered and advected.
  • Stellar feedback: we employ a kinetic stellar feedback scheme, which generates a wind with velocity scaled to the local DM dispersion, and mass loading inferred from the available SN energy for energy-driven winds (1051 erg per SNII). The metal enrichment of winds is modulated such that only 40% of local ISM metals are driven away by these galactic winds.
  • Black holes and SMBH feedback: we seed new BHs with a mass of 105 solar masses in all FoF groups which become more massive than 7x1010 solar masses. The feedback of SMBHs includes both quasar-mode and radio-mode regimes, and a consideration of the radiation field of AGNs, which heats surrounding halo gas, modifying its ionization state and net cooling rate.

The initial conditions assume a LCDM cosmology consistent with WMAP-9 measurements, from which a linear power spectrum is used to create a random realization in a periodic box with side length 75 Mpc/h = 106.5 Mpc, at a starting redshift of 127. A series of simulations are run at different resolutions, and a second set is run with only dark matter. The main simulation initially has 18203 = 6,028,568,000 hydrodynamic cells, and the same number of DM particles and MC tracers (see table for more details, including mass resolutions and gravitational softening lengths). Evolving the main simulation to z=0 used 8,192 compute cores, a peak memory of 25 TB, and 19 million CPU hours.

Preliminary Results

For details, please see the three companion papers which present the first key results: Vogelsberger et al. (2014a+b) and Genel et al. (2014), in addition to the many other topics which have been investigated using the Illustris simulation, listed on the Results page. We highlight some of the results from those papers below:

  • Comparison to stellar observables: by construction, our model agrees well with the build-up of stellar mass measured observationally. In particular, the cosmic rate star formation density (SFRD) shows good agreement up to z~10, although we slightly over predict the amount of present day star formation. The z=0 galaxy stellar mass and stellar luminosity functions agree well with SDSS-based measurements over the stellar mass range of 109 to 1012.5 solar masses, and the r-band luminosity range of -15.0 to -24.5 magnitudes. Furthermore, the stellar mass functions and the relation between stellar and halo mass agree well with observations from z=0 to z=7. The baryonic conversion efficiency is most efficient around halo masses of 1012 solar masses, and drops rapidly for both lower and higher mass systems, due to SN and AGN feedback, respectively. Illustris predicts a maximum efficiency of ~20%, occurring at the Milky Way halo mass scale, in reasonable agreement with abundance matching results. Finally, the evolution of the galaxy specific star formation rates (sSFR) is found to be in good agreement with observations up to z=8.
  • Neutral hydrogen, molecular gas and metals in galaxies: Compared to HI observations with ALFALFA, the simulated galaxy population shows the same trend of decreasing HI richness with increasing stellar mass. Contrasting the satellites of our most massive cluster at z=0 to HI measurements of the galaxies in the Virgo cluster, we similarly observe lower HI content than the global population. When considering the molecular gas fraction of galaxies, as a function of stellar mass and redshift, we find good agreement with observations at 0 < z < 3. We measure increasing gas fractions up to z ~ 2.5, before the increase halts and the dependence on redshift flattens, while the mass dependence steepens, a prediction testable by upcoming observations. Finally, we agree well with the observed mass-metallicity relationship over nearly five orders of magnitude in stellar mass, whereby more massive galaxies contain a larger fraction of metals, including the flattening of this relation above 1011 solar masses.
  • Neutral hydrogen characteristics on large scales: we compare the HI column density distribution predicted from Illustris (z=3), as well as the metallicity distribution of DLAs (2 < z < 4), with observations. We find excellent agreement, reproducing in detail the shape and location of the metallicity peak, as well as the shape and normalization of the CDDF. The simulated amount of metals in DLAs, for instance, results from the interplay between gas accretion, outflow, and ionization from the UV background. Probing the large-scale gas distribution at z=0 remains an observational challenge. The Illustris model predicts that the majority of gas (~81%) remains in the IGM at z=0, i.e. is gravitationally unbound from any halo, but that this gas contains only a minority (~34%) of the heavy metals in the universe.
  • Satellite galaxies: we measure the abundance and spatial distribution of satellite galaxies – in particular, the radial distribution of satellites within our 10 most massive galaxy clusters at z=0.3. This profile is in good agreement with the profile shape and normalization as observed in SDSS groups and clusters at 0.15 < z < 0.4. The profile lies above the DM-only result in the central region, as the stellar component is more resistant to tidal disruption at small radii due to the inclusion of dissipational baryonic physics. In physical units, the evolution of the number density profiles is weak, in agreement with existing observations out to z~1.6, and we predict that this non-evolution further extends to z>3. Studying galaxy colors and red fractions as a function of stellar mass and over-density, we see trends of both mass and environment quenching, in analogy to observations. We also note the environmental independence of the SFR of star forming galaxies, implying that high density environments increase the quenching probability, but do not directly modify the star formation rates of non-quenched systems.
  • Evolution of galaxy populations: we follow the evolutionary tracks of galaxies across cosmic time and, as a first application, study two distinct groups. (i) Selecting galaxies at z=2 and investigating their mass growth and merger histories down to z=0, we observe that z=2 galaxies lying above the star-formation main sequence are more likely to undergo a major merger by z=0, than z=2 galaxies on the SFMS. However, the most quenched galaxies at z=2 also have an increasing probability of a future major merger, indicating a relation between extreme star-formation rates and quenching. (ii) In the plane of star formation rate vs. compactness, we find a strong mass segregation in the evolutionary tracks of size evolution and quenching time scale, pointing towards the simultaneous importance of distinct physical processes.
  • Galaxy morphologies: we present a quantitative picture of the evolution of morphology from z=0 to z=5 by creating synthetic observations matching the characteristics of the Hubble Extreme Deep Field (XDF). In agreement with observations, simulated galaxies become more irregular and clumpy towards higher redshift, with increasing gas surface densities. The galaxy population in the mock UDF appears strikingly similar to the observed field, in terms of number density, colors, sizes, and morphologies. Visual classification indicates that this single simulation produces both star-forming, blue disk-like galaxies and quenched, red elliptical galaxies. A kinematic bulge-disk decomposition of all systems with stellar mass greater than 1011 stellar masses indicates that the relative fractions of early-type and late-type galaxies as a function of mass are in good agreement with recent observational samples. Studying intrinsic (g-r) and (u-i) colors, neglecting dust, we find a bimodal color distribution, albeit with an under abundance of red objects when compared to observations.
  • Galaxy circular velocity profiles: we find that z=0 galaxies in Illustris have rising or flat rotation curves inside the radius where the baryonic component dominates, and that the peak circular velocity is larger than V200 of the halo, as a result of baryon-induced contraction. This is in good agreement with observations. At redshifts greater than one, however, we predict that the stellar component is significantly more concentrated, which drives an even larger peak in the circular velocity at smaller radii, around one percent of R200 of the halo. The gas contribution to the circular velocity is always sub-dominant. Finally, based on a small subset of spiral galaxies, our model reproduces the slope and spread of the observed stellar and baryonic Tully-Fisher relations reasonably well.
  • How baryons impact DM-only results: we find that the dark matter halo responds to baryonic effects in a non-monotonic and mass dependent manner. The peak induced contraction of halos occurs for galaxy stellar masses of approximately 1011 solar masses at z=0, and for higher stellar masses at higher redshifts. On larger scales, baryonic processes modify the matter power spectrum P(k), particularly relevant on the scales targeted by upcoming weak lensing surveys such as EUCLID, which aim to make a 1% level measurement. We find that on scales smaller than k ~ 0.7/Mpc, AGN-driven outflows reduce the total power by up to 40% compared to DM-only simulations, while on even smaller scales (k ~ 70/Mpc) the power is enhanced by a factor of a few. Neither DM-only simulations nor empirical fitting models for the evolution of P(k) accurately describe the result obtained with large volume, hydrodynamic simulations.

Future Directions

Particular attention to specific discrepancies between the simulation and well-measured observational quantities will allow us to better understand our galaxy formation model. In the course of understanding the preliminary results of the Illustris simulation, we have identified three issues requiring further work: (1) even given the very energetic feedback processes that we employ, additional suppression of galaxy masses is required both in low-mass and high-mass halos. We believe this motivates the development of improved models for sub-grid feedback prescriptions, both stellar and AGN driven. (2) We find that at low redshift, massive halos around 1013 solar masses are essentially devoid of gas as a result of radio-mode AGN feedback, in disagreement with observations. Given the first point, there then remains no way to match both the stellar and gaseous content of such massive systems by adjusting the energetics of our model alone. We believe a resolution to this issue will require modifications to e.g. the duty cycle or burst-character of the radio-mode, the development of a fundamentally more sophisticated model, and/or the inclusion of presently neglected physics. (3) The stellar ages of low-mass galaxies are incorrect. Although the general trend is captured, galaxies with stellar masses below approximately 1010.5 solar masses are a factor of 2-3 too old when compared to observations. These systems form their stars too early - the challenge is how to create relatively young stellar populations in these galaxies by delaying the bulk of their star formation until significantly later times.

Although it is possible that additional fine-tuning of model parameters could alleviate some of the tensions listed above, these issues likely point to key assumptions or methods in our models which require revision. They therefore generate a path forward, identifying the areas in which we can push the physical realism and fidelity of our models even further, and hinting at the direction of the next generation of cosmological hydrodynamic simulations.