We introduce the first two simulations of the IllustrisTNG project, a next generation of cosmological magnetohydrodynamical simulations, focusing on the optical colors of galaxies. We explore TNG100, a rerun of the original Illustris box, and TNG300, which includes 2×2500 3 resolution elements in a volume twenty times larger. Here we present first results on the galaxy color bimodality at low redshift. Accounting for the attenuation of stellar light by dust, we compare the simulated (g-r) colors of 10 9 < M /M < 10 12.5 galaxies to the observed distribution from the Sloan Digital Sky Survey (SDSS). We find a striking improvement with respect to the original Illustris simulation, as well as excellent quantitative agreement with the observations, with a sharp transition in median color from blue to red at a characteristic M ∼ 10 10.5 M . Investigating the build-up of the color-mass plane and the formation of the red sequence, we demonstrate that the primary driver of galaxy color transition is supermassive blackhole feedback in its low-accretion state. Across the entire population the median color transition timescale ∆t green is ∼ 1.6 Gyr, a value which drops for increasingly massive galaxies. We find signatures of the physical process of quenching: at fixed stellar mass, redder galaxies have lower SFRs, gas fractions, and gas metallicities; their stellar populations are also older and their large-scale interstellar magnetic fields weaker than in bluer galaxies. Finally, we measure the amount of stellar mass growth on the red sequence. Galaxies with M > 10 11 M which redden at z < 1 accumulate on average ∼ 25% of their final z = 0 mass post-reddening; at the same time, ∼ 18% of such massive galaxies acquire half or more of their final stellar mass while on the red sequence.

We present the full public release of all data from the TNG100 and TNG300 simulations of the Illus-trisTNG project. IllustrisTNG is a suite of large volume, cosmological, gravo-magnetohydrodynamical simulations run with the moving-mesh code Arepo. TNG includes a comprehensive model for galaxy formation physics, and each TNG simulation selfconsistently solves for the coupled evolution of dark matter, cosmic gas, luminous stars, and supermassive blackholes from early time to the present day, z = 0. Each of the flagship runs -TNG50, TNG100, and TNG300 -are accompanied by halo/subhalo catalogs, merger trees, lower-resolution and darkmatter only counterparts, all available with 100 snapshots. We discuss scientific and numerical cautions and caveats relevant when using TNG.The data volume now directly accessible online is ∼750 TB, including 1200 full volume snapshots and ∼80,000 high time-resolution subbox snapshots. This will increase to ∼1.1 PB with the future release of TNG50. Data access and analysis examples are available in IDL, Python, and Matlab. We describe improvements and new functionality in the webbased API, including on-demand visualization and analysis of galaxies and halos, exploratory plotting of scaling relations and other relationships between galactic and halo properties, and a new Jupyter-Lab interface. This provides an online, browserbased, near-native data analysis platform enabling user computation with local access to TNG data, alleviating the need to download large datasets.

Accurate numerical solutions of the equations of hydrodynamics play an ever more important role in many fields of astrophysics. In this work, we reinvestigate the accuracy of the moving-mesh code Arepo and show how its convergence order can be improved for general problems. In particular, we clarify that for certain problems Arepo only reaches first-order convergence for its original formulation. This can be rectified by simple modifications we propose to the time integration scheme and the spatial gradient estimates of the code, both improving the accuracy of the code. We demonstrate that the new implementation is indeed second-order accurate under the L 1 norm, and in particular substantially improves conservation of angular momentum. Interestingly, whereas these improvements can significantly change the results of smooth test problems, we also find that cosmological simulations of galaxy formation are unaffected, demonstrating that the numerical errors eliminated by the new formulation do not impact these simulations. In contrast, simulations of binary stars followed over a large number of orbital times are strongly affected, as here it is particularly crucial to avoid a long-term build up of errors in angular momentum conservation.

We discuss new methods to integrate the cosmic ray (CR) evolution equations coupled to magneto-hydrodynamics (MHD) on an unstructured moving mesh, as realised in the massively parallel arepo code for cosmological simulations. We account for diffusive shock acceleration of CRs at resolved shocks and at supernova remnants in the interstellar medium (ISM), and follow the advective CR transport within the magnetised plasma, as well as anisotropic diffusive transport of CRs along the local magnetic field. CR losses are included in terms of Coulomb and hadronic interactions with the thermal plasma. We demonstrate the accuracy of our formalism for CR acceleration at shocks through simulations of plane-parallel shock tubes that are compared to newly derived exact solutions of the Riemann shock tube problem with CR acceleration. We find that the increased compressibility of the post-shock plasma due to the produced CRs decreases the shock speed. However, CR acceleration at spherically expanding blast waves does not significantly break the self-similarity of the Sedov-Taylor solution; the resulting modifications can be approximated by a suitably adjusted, but constant adiabatic index. In first applications of the new CR formalism to simulations of isolated galaxies and cosmic structure formation, we find that CRs add an important pressure component to the ISM that increases the vertical scale height of disk galaxies, and thus reduces the star formation rate. Strong external structure formation shocks inject CRs into the gas, but the relative pressure of this component decreases towards halo centres as adiabatic compression favours the thermal over the CR pressure.