We have performed the largest ever particle simulation of a Milky Way sized dark matter halo, and present the most comprehensive convergence study for an individual dark matter halo carried out thus far. We have also simulated a sample of six ultrahighly resolved Milky Way sized haloes, allowing us to estimate the halo-to-halo scatter in substructure statistics. In our largest simulation, we resolve nearly 300 000 gravitationally bound subhaloes within the virialized region of the halo. Simulations of the same object differing in mass resolution by factors of up to 1800 accurately reproduce the largest subhaloes with the same mass, maximum circular velocity and position, and yield good convergence for the abundance and internal properties of dark matter substructures. We detect up to four generations of subhaloes within subhaloes, but contrary to recent claims, we find less substructure in subhaloes than in the main halo when regions of equal mean overdensity are compared. The overall substructure mass fraction is much lower in subhaloes than in the main halo. Extrapolating the main halo's subhalo mass spectrum down to an Earth mass, we predict the mass fraction in substructure to be well below 3 per cent within 100 kpc, and to be below 0.1 per cent within the solar circle. The inner density profiles of subhaloes show no sign of converging to a fixed asymptotic slope and are well fitted by gently curving profiles of Einasto form. The mean concentrations of isolated haloes are accurately described by the fitting formula of Neto et al. down to maximum circular velocities of 1.5 km s −1 , an extrapolation over some five orders of magnitude in mass. However, at equal maximum circular velocity, subhaloes are more concentrated than field haloes, with a characteristic density that is typically ∼2.6 times larger and increases with decreasing distance from halo centre.

We introduce the Illustris Project, a series of large-scale hydrodynamical simulations of galaxy formation. The highest resolution simulation, Illustris-1, covers a volume of (106.5 Mpc) 3 , has a dark mass resolution of 6.26 × 10 6 M , and an initial baryonic matter mass resolution of 1.26 × 10 6 M . At z = 0 gravitational forces are softened on scales of 710 pc, and the smallest hydrodynamical gas cells have an extent of 48 pc. We follow the dynamical evolution of 2 × 1820 3 resolution elements and in addition passively evolve 1820 3 Monte Carlo tracer particles reaching a total particle count of more than 18 billion. The galaxy formation model includes: primordial and metal-line cooling with self-shielding corrections, stellar evolution, stellar feedback, gas recycling, chemical enrichment, supermassive black hole growth, and feedback from active galactic nuclei. Here we describe the simulation suite, and contrast basic predictions of our model for the present day galaxy population with observations of the local universe. At z = 0 our simulation volume contains about 40, 000 well-resolved galaxies covering a diverse range of morphologies and colours including early-type, late-type and irregular galaxies. The simulation reproduces reasonably well the cosmic star formation rate density, the galaxy luminosity function, and baryon conversion efficiency at z = 0. It also qualitatively captures the impact of galaxy environment on the red fractions of galaxies. The internal velocity structure of selected well-resolved disk galaxies obeys the stellar and baryonic Tully-Fisher relation together with flat circular velocity curves. In the well-resolved regime the simulation reproduces the observed mix of early-type and late-type galaxies. Our model predicts a halo mass dependent impact of baryonic effects on the halo mass function and the masses of haloes caused by feedback from supernova and active galactic nuclei.

We introduce an updated physical model to simulate the formation and evolution of galaxies in cosmological, large-scale gravity+magnetohydrodynamical simulations with the moving mesh code AREPO. The overall framework builds upon the successes of the Illustris galaxy formation model, and includes prescriptions for star formation, stellar evolution, chemical enrichment, primordial and metal-line cooling of the gas, stellar feedback with galactic outflows, and black hole formation, growth and multi-mode feedback. In this paper we give a comprehensive description of the physical and numerical advances which form the core of the IllustrisTNG (The Next Generation) framework. We focus on the revised implementation of the galactic winds, of which we modify the directionality, velocity, thermal content, and energy scalings, and explore its effects on the galaxy population. As described in earlier works, the model also includes a new black hole driven kinetic feedback at low accretion rates, magnetohydrodynamics, and improvements to the numerical scheme. Using a suite of (25 Mpc h −1 ) 3 cosmological boxes we assess the outcome of the new model at our fiducial resolution. The presence of a self-consistently amplified magnetic field is shown to have an important impact on the stellar content of 10 12 M haloes and above. Finally, we demonstrate that the new galactic winds promise to solve key problems identified in Illustris in matching observational constraints and affecting the stellar content and sizes of the low mass end of the galaxy population.

The IllustrisTNG project is a new suite of cosmological magneto-hydrodynamical simulations of galaxy formation performed with the AREPO code and updated models for feedback physics. Here we introduce the first two simulations of the series, TNG100 and TNG300, and quantify the stellar mass content of about 4000 massive galaxy groups and clusters (10 13 M 200c /M 10 15 ) at recent times (z 1). The richest clusters have half of their total stellar mass bound to satellite galaxies, with the other half being associated with the central galaxy and the diffuse intra-cluster light. Haloes more massive than about 5 × 10 14 M have more diffuse stellar mass outside 100 kpc than within 100 kpc, with power-law slopes of the radial mass density distribution as shallow as the dark matter's ( −3.5 α 3D−3). Total halo mass is a very good predictor of stellar mass, and vice versa: at z = 0, the 3D stellar mass measured within 30 kpc scales as ∝ (M 500c ) 0.49 with a ∼ 0.12 dex scatter. This is possibly too steep in comparison to the available observational constraints, even though the abundance of TNG less massive galaxies ( 10 11 M in stars) is in good agreement with the measured galaxy stellar mass functions at recent epochs. The 3D sizes of massive galaxies fall too on a tight (∼0.16 dex scatter) power-law relation with halo mass, with r stars 0.5 ∝ (M 200c ) 0.53 . Even more fundamentally, halo mass alone is a good predictor for the whole stellar mass profiles beyond the inner few kpc, and we show how on average these can be precisely recovered given a single mass measurement of the galaxy or its halo.

Hydrodynamical simulations of galaxy formation have now reached sufficient volume to make precision predictions for clustering on cosmologically relevant scales. Here we use our new IllustrisTNG simulations to study the non-linear correlation functions and power spectra of baryons, dark matter, galaxies and haloes over an exceptionally large range of scales. We find that baryonic effects increase the clustering of dark matter on small scales and damp the total matter power spectrum on scales up to k ∼ 10 h Mpc −1 by 20%. The non-linear two-point correlation function of the stellar mass is close to a power-law over a wide range of scales and approximately invariant in time from very high redshift to the present. The two-point correlation function of the simulated galaxies agrees well with SDSS at its mean redshift z 0.1, both as a function of stellar mass and when split according to galaxy colour, apart from a mild excess in the clustering of red galaxies in the stellar mass range 10 9 − 10 10 h −2 M . Given this agreement, the TNG simulations can make valuable theoretical predictions for the clustering bias of different galaxy samples. We find that the clustering length of the galaxy auto-correlation function depends strongly on stellar mass and redshift. Its power-law slope γ is nearly invariant with stellar mass, but declines from γ ∼ 1.8 at redshift z = 0 to γ ∼ 1.6 at redshift z ∼ 1, beyond which the slope steepens again. We detect significant scaledependencies in the bias of different observational tracers of large-scale structure, extending well into the range of the baryonic acoustic oscillations and causing nominal (yet fortunately correctable) shifts of the acoustic peaks of around ∼ 5%.