We introduce the first two simulations of the IllustrisTNG project a next generation of cosmological magnetohydrodynamical simulations focusing on the optical colours of galaxies. We explore TNG100 a rerun of the original Illustris box and TNG300 which includes 2  2500SUP3SUP resolution elements in a volume 20 times larger. Here we present first results on the galaxy colour bimodality at low redshift. Accounting for the attenuation of stellar light by dust we compare the simulated g  r colours of 10SUP9SUP lt MSUBSUBMSUBSUB lt 10SUP12.5SUP galaxies to the observed distribution from the Sloan Digital Sky Survey. 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 colour from blue to red at a characteristic MSUBSUB  10SUP10.5SUP MSUBSUB. Investigating the buildup of the colourmass plane and the formation of the red sequence we demonstrate that the primary driver of galaxy colour transition is supermassive black hole feedback in its low accretion state. Across the entire population the median colour transition timescale tSUBgreenSUB 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 star formation rates gas fractions and gas metallicities their stellar populations are also older and their largescale interstellar magnetic fields weaker than in bluer galaxies. Finally we measure the amount of stellar mass growth on the red sequence. Galaxies with MSUBSUB gt 10SUP11SUP MSUBSUB which redden at z lt 1 accumulate on average 25 per cent of their final z  0 mass postreddening at the same time 18 per cent of such massive galaxies acquire half or more of their final stellar mass while on the red sequence.
['galaxies: evolution', 'galaxies: formation', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Cosmology and Nongalactic Astrophysics']
First results from the IllustrisTNG simulations the galaxy colour bimodality
{'First results from the IllustrisTNG simulations: the galaxy color bimodality': 'Dylan Nelson 1 ? , Annalisa Pillepich 2 , Volker Springel 3 ; 4 , Rainer Weinberger 3 , Lars Hernquist 5 , Rudiger Pakmor 3 , Shy Genel 6 , Paul Torrey 7 , Mark Vogelsberger 7 , Guinevere Kauffmann 1 , Federico Marinacci 7 , Jill Naiman 5  \n- 1 Max-Planck-Institut fur Astrophysik, Karl-Schwarzschild-Str. 1, 85741 Garching, Germany\n- 2 Max-Planck-Institut fur Astronomie, Konigstuhl 17, 69117 Heidelberg, Germany\n- 3 Heidelberg Institute for Theoretical Studies, Schloss-Wolfsbrunnenweg 35, 69118 Heidelberg, Germany\n- 4 Zentrum fur Astronomie der Universitat Heidelberg, ARI, Monchhofstr. 12-14, 69120 Heidelberg, Germany\n- 5 Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA, 02138, USA\n- 6 Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA  \n7 Kavli Institute for Astrophysics and Space Research, Department of Physics, MIT, Cambridge, MA, 02139, USA  \n1 February 2018', 'ABSTRACT': '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 GLYPH<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 GLYPH<12> < 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 ? GLYPH<24> 10 10 : 5 M GLYPH<12> . 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 GLYPH<1> t green is GLYPH<24> 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 GLYPH<12> which redden at z < 1 accumulate on average GLYPH<24> 25% of their final z = 0 mass post-reddening; at the same time, GLYPH<24> 18% of such massive galaxies acquire half or more of their final stellar mass while on the red sequence.  \nKey words: galaxies: evolution - galaxies: formation', '1 INTRODUCTION': "Numerical simulations of cosmological galaxy formation have entered a new era. From the pioneering beginnings of dark-matter only simulations (Frenk et al. 1983; Davis et al. 1985; Navarro et al. 1997), through milestone projects of the past decade (Evrard et al. 2002; Springel et al. 2005; Klypin et al. 2011), it has now become possible to also model the hydrodynamics of cosmic gas together with the principal baryonic physics driving the formation and evolution of galaxies. Increasingly sophisticated physical models are used to realize simulations of ever increasing scope. Larger vol-  \numes enable statistically robust comparisons of entire galaxy populations, while improved numerical methods combined with higher resolution realizations provide increasingly realistic results. The emergence of cosmological hydrodynamical simulations as powerfully predictive theoretical models is embodied in the recent Illustris (Vogelsberger et al. 2014a,b; Genel et al. 2014; Sijacki et al. 2015) and EAGLE (Schaye et al. 2015; Crain et al. 2015) simulations. In concert with other large-volume efforts (Dubois et al. 2016; Dolag et al. 2016; Dav'e et al. 2016) these projects have convincingly demonstrated that hydrodynamical simulations of structure formation at kilo-parsec spatial resolution can reasonably reproduce the fundamental properties and scaling relations of observed galaxies. This zeroth order agreement has enabled many  \nfruitful investigations; it has, at the same time, revealed many informative shortcomings.  \nHere we present IllustrisTNG, a simulation campaign which (i) retains the fundamental approach and physical model flavor of Illustris, (ii) alleviates many Illustris model deficiencies with respect to benchmark observations, and (iii) significantly expands the scope with simulations of larger volumes, at higher resolution, and with new physics. In doing so we continue to heavily rely on subresolution models for physics below scales of a few hundred parsecs: by definition, this includes the process of star formation, the details of supernovae feedback, and the growth and evolution of supermassive blackholes.  \nIn order to improve upon our existing physical model, one of our first tasks is therefore to address regimes in which it fails to provide a reasonable match to observational data. To this end we identified a set of clear model deficiencies in the original Illustris simulation and specifically addressed them in the new TNG model (Pillepich et al. 2017; Weinberger et al. 2017). The principal model changes are revised treatments of both stellar feedback driven galactic winds and supermassive blackhole feedback. In this new context, we now investigate a key issue which existed in the previous simulations: the lack of a clear color bimodality between two distinct red and blue galaxy populations at low redshift. In Vogelsberger et al. (2014a) we presented the distributions of intrinsic (u-i) and (g-r) colors of simulated Illustris galaxies. Broadly, we concluded that the red peak for quiescent galaxies was either entirely absent, or not sufficiently pronounced, depending on stellar mass, while the transition from predominantly blue to red systems occurred at too large M ? , and evidence for a bimodal color distribution at any mass was weak. We principally attributed this failure to the inability of Illustris blackhole feedback to entirely quench residual star formation in massive halos.  \nThere are surprisingly few other investigations of color across the full galaxy mass range in cosmological hydrodynamical simulations. Analyzing the EAGLE simulation, Trayford et al. (2015) is by far the most comprehensive, demonstrating favorable comparison to low-redshift color-magnitude and luminosity function GAMAresults, although this involves a combination of the fiducial volume with smaller box, higher resolution, re-calibrated models; despite good qualitative agreement, the coherent level of quantitative agreement is by no means perfect. In Trayford et al. (2016) the authors explore the evolution of the color-mass space, making the connection to both environmental-driven and blackhole feedbackdriven color transition. Previously, Cen (2014) mainly studied the reddening of satellite galaxies due to hydrodynamic processes, without including blackholes or comparing to observed color distributions. Gabor & Dav'e (2012) used a phenomenological halo quenching model to study the formation of the red sequence, likewise omitting blackholes and attributing both central and satellite quenching to environments dominated by virialized gas. An update in Dav'e et al. (2017) where the enforced halo quenching mass varies with redshift also produces red central galaxies at the expense of excessively red satellites.  \nThe emergence of a color bimodality was quickly identified as a critical challenge for semi-analytical models of galaxy formation as well (Cole et al. 2000). Specifically, models which introduced a central heating process from a supermassive blackhole were then able to recover the red colors of the most massive galaxies (Kang et al. 2005; Menci et al. 2005; Croton et al. 2006; Cattaneo et al. 2007). The co-evolution of a central blackhole and its host galaxy is an important ingredient (Bower et al. 2006; Somerville et al. 2008; Benson 2012), and more recent efforts continue to explore  \ngalaxy color as a fundamental observational constraint (e.g. Henriques et al. 2015; Zu & Mandelbaum 2017).  \nObservationally, the color of a galaxy is a fundamental distillation of the emergent spectral distribution of light emitted by its constituent stellar populations. Measurements of galaxy colors are one of the most direct probes of the physical characteristics of extragalactic stars, with a correspondingly rich observational history. Modern large-scale surveys - namely, the Sloan Digital Sky Survey (SDSS) - enabled the discovery (Strateva et al. 2001) and early characterization (Kauffmann et al. 2003a; Blanton et al. 2003) of the bimodal distribution of galaxy colors at low redshift. Based on visual classification, Strateva et al. (2001) demonstrated the rough relationship between color and morphology: that members of the red population are generally more spheroidal, while members of the blue population are more disk-like. Baldry et al. (2004) modeled the color distribution of SDSS galaxies using a double-Gaussian formalism, measuring the separate red and blue luminosity and mass functions. Beyond the local universe, Bell et al. (2004) confirmed the dual population picture at 0 : 2 < z < 1 : 0 using the COMBO17 survey, constraining the evolution of the red population. Subsequently, Faber et al. (2007) used COMBO-17 and DEEP2 to contrast the z = 0 and z = 1 red and blue galaxy populations, constraining the evolution of galaxies between, and within, each. More recently, Taylor et al. (2015) presented a sophisticated modeling of the overlapping red and blue galaxy populations from z < 0 : 1 results of the GAMA survey, quantifying their respective properties and highlighting the care needed in their separation.  \nBeyond the fundamental relation of galaxy color and stellar mass, observations have also explored correlations between, and implications arising from, other galactic properties. Studying environment, Kauffmann et al. (2004) and Baldry et al. (2006) used SDSS to demonstrate that while the location of the two populations in the D4000/color-mass plane does not correlate with (projected neighbor density) environment, the relative fraction in each population does. Thereafter Peng et al. (2010) advocated for the separability of mass and environmental (i.e. satellite) driven color transformation to z = 1, quantifying an increasing red fraction with overdensity at fixed stellar mass. Considering structural properties, Kauffmann et al. (2003b) found that stellar surface mass density and stellar concentration vary strongly with mass, separating at a characteristic threshold of M ? GLYPH<24> 3 GLYPH<2> 10 10 M GLYPH<12> likely established by baryonic feedback processes. Focusing on morphology, Schawinski et al. (2014) used Galaxy Zoo classifications and found two distinct transition pathways, with different physical timescales, correlated with early versus late-type morphology. Leveraging the CANDELS survey up to z = 3, Pandya et al. (2016) measured the prevalence of the intermediary population, as well as the occupation timescales and physical properties of its members. Recently, Powell et al. (2017) measured the color-mass distributions at z = 1 separated by morphological type and argued for two distinct timescales and the possible role of blackhole feedback induced transformation.  \nIn this paper we have two goals: to present the new TNG simulations, and to quantitatively assess their agreement with the z = 0 observed color distribution, as a benchmark for improvement over the original Illustris simulation. Given this newly accessible theoretical regime, we then proceed to study the characteristics and origin of the galaxy color bimodality. In Section 2 we describe the simulations and their methods, while Section 3 details our modeling of galaxy colors and the effects of dust. Section 4 compares the results of TNG to the observational data. In Section 5 we provide a theoretical picture of the color transition process and the origin and buildup of the galaxy color bimodality. We conclude with a dis-  \ncussion in Section 6 and summarize with Section 7. Finally, Appendix A evaluates numerical convergence and details of our color modeling.", '2 THE TNG SIMULATIONS': "The IllustrisTNG project 1 is the successor of the Illustris simulation (Vogelsberger et al. 2014a,b; Genel et al. 2014; Sijacki et al. 2015). It uses an updated 'next generation' galaxy formation model which includes both new physics as well as refinements to the original Illustris model. The complete description of the TNG galaxy formation model, including details of the physical processes we include and explorations of the model and model variations on smaller test volumes are presented in the two TNG methods papers (Weinberger et al. 2017; Pillepich et al. 2017). For clarity, we avoid an exhaustive repetition of the details and here enumerate only the principal features.  \nTNG is a series of large, cosmological, gravo-magnetohydrodynamical simulations incorporating a comprehensive model for galaxy formation physics. We use the AREPO code (Springel 2010) to solve the coupled equations of self-gravity and ideal, continuum magneto-hydrodynamics (MHD; Pakmor et al. 2011; Pakmor & Springel 2013), the former computed with the TreePM approach and the latter employing a Godunov/finite-volume method with a spatial discretization based on an unstructured, moving, Voronoi tessellation of the domain. The scheme is quasiLagrangian, second order in both space and time, uses individual particle time-stepping, and has been designed to efficiently execute large, parallel astrophysical simulations on modern supercomputer architectures - our largest run used 24000 compute cores.  \nOn top of this numerical framework the TNG simulations implement models for the key physical processes relevant for galaxy formation and evolution. These are: (i) microphysical gas radiative mechanisms, including primordial and metal-line cooling and heating with an evolving ultraviolet background field, (ii) star formation in the dense interstellar medium, (iii) stellar population evolution and chemical enrichment following supernovae Ia, II, and AGB stars, individually tracking H, He, C, N, O, Ne, Mg, Si, and Fe, (iv) stellar feedback driven galactic-scale outflows, (v) the formation, merging, and accretion of supermassive blackholes, (vi) multimode blackhole feedback operating in a thermal 'quasar' mode at high accretion states and a kinetic 'wind' mode at low accretion states.  \nWe describe the complete model used to run the TNG cosmological simulations in Pillepich et al. (2017) - see Table 1 of that work for the salient differences with respect to the original Illustris model - together with Weinberger et al. (2017), focusing on the blackholes and the high-mass end. Every aspect of the physical model and detail of its numerical implementation, including all parameter values and the simulation code version, are documented therein and entirely unchanged for the simulations we present here.  \nThe TNG project is made up of three simulation volumes: TNG50, TNG100, and TNG300. Here we present the two larger volumes: the first simulation, TNG100 , includes 2 GLYPH<2> 1820 3 resolution elements in a GLYPH<24> 100 Mpc (comoving) box. The baryon mass resolution is 1 : 4 GLYPH<2> 10 6 M GLYPH<12> , the gravitational softening length of the dark matter and stars is 0.7 kpc at z =0, and the gas component  \nTable 1. Most important details of the two primary simulations presented here: TNG100 and TNG300. These are: the simulated volume and box sidelength (both comoving), the number of initial gas cells, dark matter particles, and Monte Carlo tracers. The mean baryon and dark matter particle mass resolutions, in solar masses. The minimum allowed adaptive gravitational softening length for gas cells (comoving Plummer equivalent), and the redshift zero softening of the dark matter and stellar components in physical kilo-parsecs. In Tables A1 and A2 we provide additional details, including the resolution series and dark-matter only analogs of each box.  \nhas an adaptive softening with a minimum of 185 comoving parsecs. The phases are the same as original Illustris, enabling objectby-object comparison. The second simulation, TNG300 , includes 2 GLYPH<2> 2500 3 resolution elements in a GLYPH<24> 300 Mpc box, a volume GLYPH<24> 20 times larger. Its baryon mass resolution, collisionless softening, and gas minimum softening are 1 : 1 GLYPH<2> 10 7 M GLYPH<12> , 1.5 kpc, and 370 parsecs, respectively. The key details of each volume are given in Table 1. We also run two lower resolution realizations for each volume and dark-matter only analogs of every configuration, a total of twelve simulations which are detailed, and studied, in Appendix A. The high-resolution TNG50 run will be introduced in future work (Nelson et al. in prep, Pillepich et al. in prep). It includes 2 GLYPH<2> 2160 3 resolution elements in a GLYPH<24> 50 Mpc box, with a baryon mass resolution of 8 : 5 GLYPH<2> 10 4 M GLYPH<12> , z = 0 collisionless softening of 280 parsecs, and a minimum gas softening of 74 parsecs.  \nThis is one in a series of five papers that introduce IllustrisTNG, each of the other four with a distinct scientific focus. In Naiman et al. (2017) we study neutron star-neutron star mergers and r-process chemical enrichment, and in Pillepich et al. (2017) we decompose the different stellar components at group and cluster scales and investigate the universality of stellar mass profiles. In Marinacci et al. (2017) we measure radio synchrotron emission from massive galaxy clusters, while in Springel et al. (2017) we examine the spatial clustering of matter and of galaxies.  \nFor TNG we adopt a cosmology consistent with recent constraints (Planck Collaboration 2016), namely GLYPH<10>GLYPH<3> ; 0 = 0 : 6911, GLYPH<10> m ; 0 = 0 : 3089, GLYPH<10> b ; 0 = 0 : 0486, GLYPH<27> 8 = 0 : 8159, ns = 0 : 9667 and h = 0 : 6774. Herein we identify gravitationally bound substructures using the SUBFIND algorithm (Springel et al. 2001) and link them through time with the SUBLINK merger tree algorithm (RodriguezGomez et al. 2015).", '3 MODELING GALAXY COLORS': 'Weconsider three models for computing the colors of galaxies: (A), (B), and (C). Each builds off the previous with increasing complex-', '4 D. Nelson et al.': "ity. Model (A) sums up the intrinsic emission of each star particle in a galaxy, with no attenuation. Model (B) adds 'unresolved' dust attenuation due to dense gas birth clouds surrounding young stellar populations. Model (C) adds 'resolved' dust attenuation due to the simulated distribution of neutral gas and metals. Therefore, the first two models are invariant to orientation, while colors from the third depend on viewing angle. In this paper all results use model (C) exclusively , while the three models are compared in the appendix.", '3.1 Model (A) - no dust': 'Each star particle in the simulation is modeled as a single-burst simple stellar population (SSP), whose birth time, metallicity, and mass are recorded. We use the FSPS stellar population synthesis code (Conroy et al. 2009; Conroy & Gunn 2010; Foreman-Mackey et al. 2014) with the Padova isochrones, MILES stellar library, and assuming a Chabrier IMF. The isochrones contain 22 metallicity steps from GLYPH<0> 3 : 7 to GLYPH<0> 1 : 5 (log total metal mass ratio), and 94 age steps from GLYPH<0> 3 : 5 to 1 : 15 (log Gyr). For each point on this 2D grid, we convolve the resulting population spectrum with each of the observed pass-bands (SDSS u,g,r,i,z airmass 1.3), saving the total magnitude in each. Default models for dust emission and nebular line and continuum emission are included.  \nFor each star particle we then use bicubic interpolation in this ( Zi ; t ssp) plane to obtain m 1, the band magnitude of a 1 M GLYPH<12> equivalent population. Multiplying by the initial mass Mi gives the total band magnitude of the actual population m = m 1 GLYPH<0> 2 : 5 GLYPH<3> log 10 ( Mi = M GLYPH<12> ). For the total magnitude of each galaxy, we sum all the subhalo member star particles, with no additional restrictions. Galaxy colors are the difference between two of these total band magnitudes.', '3.2 Model (B) - unresolved dust': 'Following Charlot & Fall (2000) we add a simple powerlaw extinction model for the attenuation by finite lifetime birth clouds surrounding young stellar populations as well as the ambient diffuse ISM. The intrinsic luminosity L i( GLYPH<21> ) of each stellar population is attenuated as L obs( GLYPH<21> ) = L i( GLYPH<21> ) = e GLYPH<0> GLYPH<28> B GLYPH<21> where the wavelength dependent optical depth GLYPH<28> B is then given as  \nGLYPH<21>  \nGLYPH<28> B GLYPH<21> = 8 > > < > > : GLYPH<28> 1( GLYPH<21>=GLYPH<21> 0) GLYPH<0> GLYPH<11> 1 ; t ssp 6 t age GLYPH<28> 2( GLYPH<21>=GLYPH<21> 0) GLYPH<0> GLYPH<11> 2 ; t ssp > t age 9 > > = > > ; : (1)  \nAll parameters are unchanged from their commonly adopted values. We take GLYPH<28> 1 = 1 : 0, GLYPH<28> 2 = 0 : 3 (although see Wild et al. 2011), GLYPH<11> 1 = GLYPH<11> 2 = 0 : 7, t age = 10 Myr, and GLYPH<21> 0 = 550 nm. The resulting colors remain independent of viewing angle.', '3.3 Model (C) - resolved dust': "We follow the distribution of neutral gas in and around each simulated galaxy. Using its metal content we can then model the subsequent dust attenuation of the optical light emitted by stars. We reiterate that this obscuration is in addition to the previously described Model (B) component. Following Guo & White (2009) and Xu et al. (2017) we likewise adopt the results of Calzetti et al. (1994) and take their internal dust model, where the ionized gas and dust are co-spatial and uniformly mixed (see the Appendix for  \nan exploration of the alternative geometrical assumption of a foreground screen). In the present case  \nL obs( GLYPH<21> ) = L i( GLYPH<21> ) GLYPH<1> 1 GLYPH<28>GLYPH<21> GLYPH<0> 1 GLYPH<0> e GLYPH<0> GLYPH<28>GLYPH<21> GLYPH<1> (2)  \nwhere GLYPH<28>GLYPH<21> is the total dust attenuation arising from both absorption and scattering. We use the resolved gas distribution of the simulation to calculate GLYPH<28> a GLYPH<21> , the contribution from absorption. The scattering contribution is then incorporated by deriving a total GLYPH<28>GLYPH<21> as  \nGLYPH<28>GLYPH<21> = f e GLYPH<11> ( GLYPH<28> a GLYPH<21> ; h GLYPH<21> ; !GLYPH<21> ) = GLYPH<28> a GLYPH<21> h h GLYPH<21> (1 GLYPH<0> !GLYPH<21> ) 1 = 2 + (1 GLYPH<0> h GLYPH<21> )(1 GLYPH<0> !GLYPH<21> ) i (3)  \nwhere h GLYPH<21> and !GLYPH<21> are taken from fits to the scattering anisotropy weight factor and albedo over GLYPH<21> 2 [100 nm, 1000 nm], respectively (Calzetti et al. 1994). For the absorption optical depth, we take  \nGLYPH<28> a GLYPH<21> = A GLYPH<21> AV ! GLYPH<12> (1 + z ) GLYPH<12> ( Zg = Z GLYPH<12> ) GLYPH<13> ( N H = N H ; 0) : (4)  \nHere, the first term is the solar neighborhood extinction curve which we take from Cardelli et al. (1989). Note that in the optical the differences between the empirical extinction curves of the SMC, LMC, and Milky Way are negligible. The second and third terms parameterize a redshift and metallicity dependent dust-to-gas ratio (see also McKinnon et al. 2016). We adopt GLYPH<12> = GLYPH<0> 0 : 5 and a broken power-law GLYPH<13> = 1 : 6 ( GLYPH<21> > 200 nm) GLYPH<13> = 1 : 35 ( GLYPH<21> < 200 nm, outside the SDSS bands). We take Z GLYPH<12> = 0 : 0127 and for the normalization of the hydrogen column N H ; 0 = 2 : 1 GLYPH<2> 10 21 cm GLYPH<0> 2 . All parameters are taken unchanged from their literature values.  \nThe quantities Zg and N H are both computed independently for each star with the following procedure. For each galaxy, we trace through the resolved gas distribution from twelve viewing angles, corresponding to the verticies of the Ns =1 Healpix sphere (G'orski et al. 2005) oriented in simulation coordinates (i.e. randomly with respect to each galaxy). For each direction, an orthographic projection grid with a pixel scale of 1 : 0 kpc is created covering all the member star particles. All gas cells and star particles in the subhalo are z-sorted by distance from the view point, then gas cells are progressively deposited on to the grid front-to-back using the standard cubic spline kernel with h = 2 : 5 r cell = 2 : 5 (3 V cell = 4 =GLYPH<25> ) 1 = 3 . When the z-position of a star is reached it samples with bilinear interpolation the accumulated gas values in the four grid cells nearest to its projected location. These values are the neutral hydrogen column density and the neutral hydrogen mass weighted gas metallicity, which are used for N H and Zg above. Note that the neutral hydrogen fraction of each gas cell is computed during the simulation using the atomic network of Katz et al. (1996) with the time variable though spatially uniform ionizing UV background of Faucher-Gigu'ere et al. (2009) given the density-based self shielding prescription of Rahmati et al. (2013).  \nEach full stellar population spectrum is convolved with the GLYPH<21> -dependent attenuation and then by the observed passband. Finally, we impose a 3D radial restriction of at most 30 physical kpc from the subhalo center. This approximates a Petrosian aperture, and includes all subhalo particles for low-mass galaxies.  \n8  \n6  \n6  \nFigure 1. The distribution of simulated versus observed (g-r) colors. In black the SDSS z < 0.1 sample is shown, while red shows the result from the old Illustris simulation, and blue the result from our new TNG100 simulation of equivalent box size and resolution. In all cases, galaxies with stellar masses from 10 9 to 10 12 M GLYPH<12> are included, divided into six mass bins as indicated in the six panels (log M ? values shown). One-dimensional kernel density weighted PDFs are shown (with a half-Scott bandwidth of 0 : 5 N GLYPH<0> 1 = 5 ). In each case galaxies are shown regardless of if they are centrals or satellites. In addition, in the lowest mass panel we also decompose the simulated results and show centrals only (dotted lines), demonstrating that the somewhat excessive red peak at these masses arises from the red satellite population. Illustris and TNG100 galaxies are treated with the same dust model C (see Section 3). Evident in both SDSS and TNG100: the dominant blue population at low masses, the emergence of the red peak and the strong bimodality around 10 10 : 5 M GLYPH<12> , and the dropout of blue systems with the red population dominant at high masses.  \n<!-- image -->", '4 COMPARISON TO OBSERVATIONS': "To assess the outcome of the simulations we compare to the observed colors of galaxy populations in the local universe. Specifically, we draw an observed sample from the Sloan Digital Sky Survey (SDSS DR12; Alam et al. 2015) by selecting all galaxies with spectroscopic redshifts z < 0.1, a sample of 378,287 objects. For their absolute magnitudes in the five ugriz SDSS broadband filters we adopt the modelMag values. These are computed from the best fit of either an exponential or de Vaucoulers profile in the rband, whose amplitude (but not shape) is then fit in the remaining bands. They are in excellent agreement with Petrosian magnitudes for the majority of galaxies, with offsets due to color gradients and aperture corrections at the massive end (Strauss et al. 2002). We add the correction for Galactic extinction following Schlegel et al. (1998) as well as a k-correction using the method of Chilingarian et al. (2010). For stellar mass estimates for each observed galaxy we adopt the Granada FSPS wide SED fit values (Ahn et al. 2014).  \nIn Figure 1 we show normalized PDFs of (g-r) color in six equal bins of stellar mass from 10 9 to 10 12 M GLYPH<12> . The black curves in each panel show the observed distribution of galaxy colors in that mass range. In red we show the results of the original Illustris simulation, while in blue we show the outcome of the new TNG100 simulation. We adopt the same dust model ('resolved dust' model C) for both, and include all z = 0 galaxies with 10 9 < M ?= M GLYPH<12> < 10 12 . In both the observations and in TNG100 there is a consistent over-  \nall picture: at low stellar masses, the blue population dominates the sample with (g-r) GLYPH<24> 0.4, while the red peak emerges towards higher M ? at a characteristic (g-r) GLYPH<24> 0.8. The two reach roughly equal magnitude, displaying the strong observed color bimodality, at M ? ' 10 10 : 5 M GLYPH<12> . Above this mass blue systems quickly disappear, leaving the red galaxy population dominant at high mass.  \nThe agreement between TNG100 simulated galaxies and the observational dataset is in general excellent - see Section 4.2. For stellar masses below 10 11 M GLYPH<12> where there are clearly two distinct peaks, we recover their respective locations in (g-r) and the position of the minimum between them to ' 0 : 05 mag or better. The evolution of the relative amplitude of the blue and red distributions as a function of M ? is likewise well captured, as is the stellar mass of maximal bimodality, where the red fraction of the galaxy population is half. On the other hand, the original Illustris simulation does not exhibit a color bimodality sufficiently strong, or at the correct location in M ? , having generally too few red systems at intermediate masses between 10 9 : 5 < M ?= M GLYPH<12> < 10 11 . Above 10 11 M GLYPH<12> Illustris agrees well with SDSS in the peak location of the red sequence, although the distribution may be too narrow in color at the highest masses.  \nThe clearest potential tension between TNG100 and the z < 0.1 SDSS sample is that the simulation appears to contain too many red galaxies at the lowest masses, below 10 9 : 5 M GLYPH<12> . In this panel (upper left) we also show the separate color distribution for central galaxies only, which demonstrates that this overly popu-  \n<!-- image -->  \nFigure 2. The simulated stellar age versus mass (left panel) and stellar metallicity versus mass (right panel) relations, compared to several observational data sets (gray). The ages are compared to observational results from Gallazzi et al. (2005) and Bernardi et al. (2010), while the stellar metallicity relation is compared to data from Gallazzi et al. (2005) with dwarf points from Woo et al. (2008) and Kirby et al. (2013). In both cases, the simulated lines are taken at z = 0 : 1 by considering the mean age and/or metallicity of all star particles within a 4 kpc mock SDSS fiber aperture, weighted by their r-band luminosities. For Z stars we include approximate corrections from Guidi et al. (2016) for the difference between our luminosity weighted fiber-aperture technique and the spectral index fitting technique actually used in Gallazzi et al. (2005) (dashed blue and orange lines). In green we show the TNG100-1 trend and its scatter derived by fitting mock SDSS spectra with matched noise and resolution characteristics (see footnote for details).  \n<!-- image -->  \nlated red component is entirely made up of satellites. TNG100 may genuinely produce too many low-mass red satellite galaxies. These are however some of the most marginally resolved systems, containing only GLYPH<24> 1000 star particles, and in Appendix A we show that with increasing resolution the low-mass galaxy population tends to become bluer. This may be, therefore, an unresolved regime for TNG100. On the other hand, the observations are not complete at these low masses. A bias in the effective volume sampled by the SDSS spectroscopic survey as a function of color, such that intrinsically fainter red galaxies are underrepresented, could explain the discrepancy. Similarly, a bias from the spectroscopic targeting due to fiber collisions in high-density environments, such that intrinsically redder satellites are underrepresented, could also contribute. Amore sophisticated analysis - namely, a full mock SDSS spectral survey - could definitely demonstrate the impact of both issues.", '4.1 Stellar Ages and Metallicities': "The principal ingredients of the colors are the ages and metallicities of the stars. In Figure 2 we therefore independently compare the predictions for TNG100 and TNG300 at z = 0 : 1 against observed stellar age and metallicity relations. The primary constraint for both is from Gallazzi et al. (2005), which is also based on lowredshift SDSS galaxies. From the simulations we compute luminosity weighted r-band 2D fiber-aperture 4 kpc restricted values as an approximation of the SDSS fiber derived observations.  \nThe simulated stellar ages are in excellent agreement with the observed relation, both in the median values (solid lines) and in the scatter (16th-84th percentiles). We caution that the simulated scatter is intrinsic only - we have not convolved by an estimate of the observational uncertainty. The simulated scatter is therefore consistent with the observed scatter in the sense of a lower limit. The average stellar ages transition from a relatively constant ' 8 GLYPH<0> 10 Gyr at the highest masses to a relatively constant ' 2 Gyr at the lowest masses, with intermediate values over 10 10 < M ?= M GLYPH<12> < 10 11 . This same trend, and the transition point at Milky Way mass,  \nis seen in the observed ages. Tension exists in the value of the mean age of stars in the lowest mass galaxies, below 10 9 : 5 M GLYPH<12> . While Gallazzi et al. (2005) favor an asymptotic value of 1 Gyr, we find a somewhat higher mean age of 2 Gyrs. This lack of present day low-mass dwarf systems formed within the past billion years in the simulation may be a genuine failure of the model. As we discuss next, it may also arise from an overly simplistic comparison.  \nThe stellar metallicities do not reproduce the observed trend to the same degree. At 10 10 : 5 M GLYPH<12> and above the simulations agree with the observed population quite well, exhibiting a weak scaling with M ? normalized at a super-solar metallicity value. However, the trend with stellar mass is too flat, and this leads to too high metallicities at low mass. Specifically, the maximal discrepancy is ' 0 : 5 dex in stellar metallicity at log( Z ?= Z GLYPH<12> ) ' GLYPH<0> 0 : 5 (observed), at the level of the 1 GLYPH<27> (upper 84th percentile) range of the observed sample. We hypothesize that this discrepancy is driven by the rather different method for deriving the simulated versus observed stellar metallicities. Specifically, we show (orange and blue dashed lines) a proposed correction between the two from Guidi et al. (2016), which involves a dust radiative transfer computation of the emergent spectrum, followed by measurement of individual spectral features (e.g. D4000 n , H GLYPH<12> , [Mg2Fe]) and the application of the specific procedure of Gallazzi et al. (2005). This is rather distant from the simple luminosity-weighted fiber-aperture values we present, which do not involve any spectral measurements. Although we would like to use this correction directly, it has been derived from a small simulated sample over a limited mass range, and with a different physical model. A robust comparison of TNG stellar metallicities (and ages) to SDSS fiber spectrally-derived values therefore requires more sophisticated forward modeling. 2  \n2 A preliminary analysis - whereby we have created mock SDSS fiber spectra for every simulated galaxy with noise and instrumental resolution characteristics sampled randomly from real SDSS fibers, identically performed a MCMC fit on both the mock spectral catalog and the actual SDSS spectral sample, using FSPS as the generative model (Leja et al. 2017, John-  \n4.0  \n1.0  \nFigure 3. The distribution of observed (black) versus simulated (blue) galaxies in the color-mass plane at z =0, shown for three band combinations beyond the (g-r) explored in this work: (u-i), (u-r), and (r-i). Conditional 2D kernel density estimates (KDEs) are shown with contours at { 0.2,0.5,0.75,0.98 } . On the right side of each panel the marginalized one-dimensional PDF is shown, drawing a simulated sample with the same stellar mass distribution as the observations using a discrete inverse transform sampling method. We see that the (u-i) and (u-r) colors show the same level of excellent quantitative agreement with respect to the observations, while the (r-i) colors of low mass galaxies are slightly too red.  \n<!-- image -->  \nIndeed, the colors of stellar populations are sensitively dependent on both age and metallicity. Yet, we have demonstrated excellent quantitative agreement between simulated and observed galaxy colors down to, at least, stellar masses of 10 9 : 5 M GLYPH<12> . Therefore, the apparent tension in the stellar age and metallicity relations at these same masses suggests that the direct comparison of age and metallicity is presently less informative than the more indirect comparison of color, because of the difficulty of making a rigorous comparison of the simulated results with SDSS fiber-derived quantities.  \nson et al. in prep) with the MILES libraries over 3750 < GLYPH<21> < 7000 in the seven-dimensional parameter space of f z residual, M ? , Z ? , GLYPH<28> SFH, t age, GLYPH<28> 1 ; dust, GLYPH<27> disp g with broad priors on each - suggests that the correspondingly derived median stellar metallicities are offset lower by GLYPH<24> 0 : 2 dex at M ? = 10 10 M GLYPH<12> and by GLYPH<24> 0 : 5 dex at M ? = 10 9 M GLYPH<12> with respect to the presented simulation values. This implies that the apparent tension in the Z ? GLYPH<0> M ? relation arises solely from an overly simplistic comparison (Nelson et al. in prep). Note that we truncate this mock Zstars measurement at M ? ' 10 10 : 5 M GLYPH<12> , after which our current modeling implies a similar and roughly constant offset extends up to the more massive galaxies. A more robust conclusion in this regime, however, would require an improved fitting of systems with significant stellar velocity dispersions.", '4.2 Characterizing the Color-Mass Plane': 'Color combinations with the u-band will span the 4000Å break at low redshift, making them more sensitive indicators of quiescence and star formation history (Bruzual 1983). We therefore augment the (g-r) colors presented thus far with several additional comparisons against the SDSS sample. In Figure 3 we show contours of the full two dimensional distributions in the color-mass plane for four band combinations: (u-i), (g-r), (r-i), and (u-r), comparing each to the same z < 0 : 1 SDSS galaxies. We see that each of (u-r), (u-i), and (r-i) share the same level of quantitative agreement with the observations as (g-r). The blue population in (r-i) color is excessively red by roughly 0.1 dex, and since (r-i) is sensitive to dust reddening this discrepancy may indicate tension in our dust modeling.  \nMore interestingly, we see the emergence of several second order discrepancies. These would have been impossible to observe, much less characterize, without this excellent base level of agreement. Two are immediately obvious.  \nThe first is the fact that the sharp transition from blue to red occurs at slightly too high M ? in the simulations, by about 0.1 dex. The immediate conclusion would be that galaxies transition to a red-dominated population slightly too late. We caution, though, that  \nwe have not modeled mock SDSS stellar masses of galaxies at the same level of rigor as their colors (see Torrey et al. 2015; Guidi et al. 2015; Bottrell et al. 2017, for work in this direction); 0.1 dex is likely also less than the observational uncertainty in M ? .  \nThe second discrepancy is that the slope of the red sequence is flatter in the simulations than in the observations. This implies that there may be some physical effect tilting the observed color-mass relation at the high-mass end, which TNG either misses or captures less strongly. The physical origin is likely the residual trends of age and Z ? within the passive population (Loubser et al. 2009). It is possible that in TNG the mean stellar age and Z ? do not increase as strongly with stellar mass as they should. On the other hand, we find that the integrated colors of the massive red population are quite sensitive to aperture. We propose that galactic color gradients act to wash out this intrapopulation effect: because passive galaxies have measurable radial gradients in stellar age, metallicity, and color (Gonz\'alez Delgado et al. 2015; Cook et al. 2016) integrated colors are sensitive to the radius within which they are measured, and our fixed 30 kpc aperture restriction cannot precisely correspond to the SDSS model magnitudes. Smaller apertures exclude a progressively larger fraction of (generally bluer) ICL, reddening the integral color of the simulated galaxy. Similar considerations may explain why simulated galaxies in the most massive bin of Figure 1 (at M ? > 10 11 : 5 M GLYPH<12> ) do not reach the maximal (g-r) values seen in SDSS, although the lack of this tail at (g-r) > 0.9 may also reflect the inability of our dust modeling to represent the most extremely reddened observed galaxies.  \nTo avoid ambiguous statements about having \'reasonable\', \'good\', or \'excellent\' agreement between observed and simulated galaxy properties we quantify the respective color distributions. Inspired by Baldry et al. (2004) and Taylor et al. (2015) we proceed without making a priori cuts on where red versus blue galaxies populate the color-mass space. Instead we assume that, at a given stellar mass, two distinct (possibly overlapping) populations exist, each with a gaussian distribution G ( C ) in color C . In a discrete set of mass bins, we therefore fit a series of 2-gaussian mixture models  \n1 X i = 0 G i = 1 X i = 0 i + ( GLYPH<0> 1) i Ai GLYPH<27> i (2 GLYPH<25> ) GLYPH<0> 1 = 2 exp " GLYPH<0> ( C GLYPH<0> GLYPH<22> i ) 2 2 GLYPH<27> 2 i # : (5)  \nWe consider the stellar mass range spanning 9.0 < log( M ? /M GLYPH<12> ) < 12.0 with a binsize GLYPH<1> log M ? = 0.15 giving 20 total mass bins. Similarly, a color binsize of GLYPH<1> ( g GLYPH<0> r ) = 0.04 over the color range 0.0 < (g-r) < 1.0 gives 25 color bins per mass bin. Rather than let the five parameters of each P G i vary independently, we constrain their center positions and widths to be continuous functions of stellar mass according to a global linear+tanh relation, namely  \nT ( M ? ) = p 0 + p 1 log M ? + q 0 tanh " log M ? GLYPH<0> q 1 q 2 # : (6)  \nWe have 20 relative amplitudes for the red versus blue populations across the mass bins. In addition, ten T ( x ) function parameters for the red and blue GLYPH<22> i , and likewise the ten T ( x ) function parameters for the red and blue GLYPH<27> i , resulting in a forty parameter vector GLYPH<2> to fit. Rather than a simultaneous global fit, we first implement and carry out the parameter-split iterative procedure as described in Baldry et al. (2004). We then take this solution GLYPH<2> 0 as the initial guess of a global MCMC fit to explore the 40-dimensional parameter space and estimate the uncertainties of each (using EMCEE; Foreman-Mackey et al. 2013).  \n0.16  \nFigure 4. Comparison of the red and blue population fit parameters: widths GLYPH<27> i (top panel) and center locations GLYPH<22> i (bottom panel). The same procedure has been applied to the SDSS z < 0 : 1 sample and the TNG100-1 simulation at z = 0. We visually focus on the low-mass blue and high-mass red parameters by showing them as dark lines and shaded regions, which give the median and 10-90 percentile model ranges of the final chains, respectively.  \n<!-- image -->  \nApplying the same procedure to the simulated and observed data sets, Figure 4 shows the result. In it we compare the locations and widths of the red and blue populations as a function of stellar mass, in SDSS (dotted) versus TNG100-1 at z = 0 (solid). The mean color of blue galaxies is in excellent agreement, within 0.05 mag. Its location in (g-r) may increase more slowly with M ? in the observations, although we caution that this could be driven largely by the balance between the pi and qi components of the T function (i.e., the amount of curvature is not well constrained). Most interestingly, we can quantify the red sequence discrepancies noted above: in the simulation the mean of this is slightly too blue (by . 0 : 05 mag). This tension arises because the dependence of GLYPH<27> red on M ? is steeper than observed, while the dependence of GLYPH<22> red on M ? is shallower than observed. The simulated color widths are slightly less than observed, indicating that the intrinsic scatter is consistent with the observed scatter. We reiterate that we have not convolved the intrinsic colors with estimated observational uncertainties, which would be required to compare the scatter in detail. In addition, we here use (as always) the fiducial model (C) \'resolved dust\' colors for the simulations, and note that the exact level of quantitative (dis)agreement with respect to SDSS must depend at some level on the details of our dust modeling assumptions, which we explore in the Appendix.  \nFigure 5. The cosmological context offered by the TNG100 and TNG300 simulations. The main figure shows three views of TNG300-1 at z = 0, centered on the most massive cluster in the box: gas column density (left), gas temperature (center), and magnetic field amplitude (right). The first two are projected through a depth equal to a third of the volume, while the magnetic field is shown in a thin 100 kpc slice. The 50 most massive halos in the thick slice are indicated by circles with radii of r 200 ; crit. In the inset in the upper right, we show the TNG100-1 volume at z = 0 in projected dark matter density (left) and gas density (right).  \n<!-- image -->', '5 THEORETICAL INTERPRETATION': 'Together, the TNG100 and TNG300 simulations are a powerful theoretical tool. In Figure 5 we show a visual overview of the two simulation volumes at the present day, z = 0, including several properties of the gas as well as the dark matter distributions. The light of stars is entirely absent from this composite view, while all  \nthe galaxies (together with their colors) which we have thus far studied exist within this full cosmological context. By construction these large volumes include galaxies of a great variety of types, environments, and formation histories. Simultaneously, each system has a large information content, whereby a wide range of properties of its stars, gas, dark matter, or blackholes can be measured.  \nFigure 6. Six views of the color-mass plane of all central galaxies in TNG300-1 at z =0. In each case colors show a physical, median value of all the galaxies in that bin, and measurements are always within 2 r ?; 1 = 2 unless otherwise noted. We show (i) the specific star formation rate, (ii) the gas metallicity normalized by solar, (iii) the gas fraction - gas mass divided by total baryon mass, (iv) the mean r-band luminosity weighted stellar age within a 4 kpc fiber-aperture, (v) the mean magnetic field strength in the ISM, and (vi) the GLYPH<12> GLYPH<0> 1 parameter of the mean ratio of magnetic to thermal pressure in the halo (0 : 15 < r = r vir < 1 : 0). Gray indicates a zero (or n/a) value for a given property of all the galaxies in that bin. Contours show conditional 75% and 95% levels.  \n<!-- image -->  \n0.0  \nAll these components have arisen self-consistently as the result of a coupled physical evolution across the history of the universe. In this section we use these relationships to elucidate the origin of the galaxy color distribution. We frequently include the TNG300 large volume in order to take advantage of its superior statistics, and refer  \nto the Appendix for a discussion of resolution convergence and the differences between TNG300 and the higher resolution TNG100 simulation.', '5.1 Relating Color to Other Galaxy Properties': 'In Figure 6 we show six views of the same TNG300-1 color-mass plane at z = 0. In each case, we map color to one median property of all the galaxies falling in that specific bin. These are: the specific star formation rate (upper left), gas metallicity (upper right), galaxy gas fraction (center left), galaxy stellar age (center right), ISM magnetic field amplitude (lower left), and halo magnetic to gas pressure ratio (lower right). Gray colors indicate a uniformly zero value for that bin, while contours show the conditional 75% and 95% levels of the galaxy population as a function of M ? .  \nIn broad comparison of the lower left to the upper right of each panel - the locations of the blue and red populations, respectively - we recover several expected results. Red galaxies are passive, gas poor, and old. With respect to otherwise increasing trends with stellar mass, they also have lower gas metallicities and smaller magnetic field strengths, the latter in both the galaxy and the halo. Stellar age is the only property shown which depends almost entirely on color and not at all on M ? , indicated by the purely vertical gradient. The fact that more massive galaxies tend to be older therefore imprints itself as an increasing red fraction with stellar mass. On the other hand, galactic stellar metallicity has a largely horizontal gradient in this plane (not shown), reflecting the stronger dependence of color on age.  \nThe other five properties we show here, all related to the gas, have more complex behavior. Each appears driven by an underlying scaling with galaxy mass, which is then strongly modulated by color. For instance, the mean galactic gas metallicity increases with M ? following the mass-metallicity relation only while the majority of galaxies are star-forming. Above a threshold stellar mass of GLYPH<24> 10 10 : 5 M GLYPH<12> we see that Z gas is suppressed. Similar trends are clearly visible in f gas as well as j B j ISM. The emergent sSFR lacks any significant mass dependence, but otherwise has a similar trend. The halo GLYPH<12> -ratio is the most complex, with non-monotonic behavior along both axes which we shortly explore further.  \nIn the transition regime of M ? GLYPH<24> 10 10 : 5 M GLYPH<12> and at fixed stellar mass, the statistics of TNG300 allow us to see by eye the correlations between galaxy color and all five of the gas properties. For instance, in the median specific star formation rate, a vertical slice at this mass scale transitions from highly star forming (in yellow) at a (g-r) color of 0.3 to complete quiescence (in purple/gray) at a (gr) color of 0.7. Similarly for the gas fraction, which also drops from near unity to essentially zero at a fixed stellar mass of GLYPH<24> 10 10 : 5 M GLYPH<12> . We point out an explicit prediction of the TNG model: at the Milky Way mass scale blue galaxies should have a large-scale interstellar magnetic field strength of order 10-30 GLYPH<22> Gwhile red galaxies, at the same stellar mass, should host a much weaker field of . 1 GLYPH<22> G.  \nAn intriguing feature evident in Figure 6 is the diagonal boundary dividing the red and blue regimes. It implies the existence of a mass-dependent color threshold, or equivalently stated, a colordependent mass threshold, beyond which galaxies enter a distinct physical population. We show later that this diagonal separation in the CM plane does not correspond to a boundary orthogonal to the bulk flux of galaxies as they evolve, which is more vertical. Still, its existence, location in M ? , and slope are all signatures of the color transition mechanism, driven by the specifics of the quenching process of the TNG blackhole feedback model.  \nTo identify the underlying physical mechanism of color transition, Figure 7 shows a diagnostic of the integrated energy input from feedback by supermassive blackholes. Specifically, the cumulative energy released in the low accretion state (kinetic mode) is compared to the cumulative energy released in the high accretion  \n0  \nFigure 7. The role of supermassive blackhole feedback in setting the characteristic stellar mass of color transition: the logarithm of the ratio of the integrated energy injection in the low versus high accretion states. These two quantities are summed across all blackhole progenitors. The median across all central galaxies is given for TNG100-1 (blue) and TNG300-1 (orange). This ratio is always less than one, increasing rapidly between 10.0 < log( M ?= M GLYPH<12> ) < 10.5, similar to the transition from a blue to red dominated galaxy population. In addition we show the value of GLYPH<31> ( M BH) defining the fraction of the Eddington accretion rate below which a blackhole is in the low state (dashed gray), and the actual z = 0 blackhole accretion rates (dashed green), which intersect at this same characteristic mass.  \n<!-- image -->  \nFigure 8. The correlation between the GLYPH<12> GLYPH<0> 1 = P B = P gas pressure ratio of the halo gas and the (g-r) color of the galaxy at a stellar mass of 10 10 : 5 M GLYPH<12> . We include all central galaxies for TNG100-1 (blue median line and individual galaxies indicated by circles) as well as TNG300-1 (orange median line and 10 to 90 percentile shaded region) in a mass bin of width GLYPH<6> 0 : 1 dex. Blue galaxies have GLYPH<12> GLYPH<0> 1 GLYPH<24> 1, which then increases moderately with increasing color for transitional systems, before reversing and plummeting to GLYPH<12> GLYPH<0> 1 GLYPH<24> 10 GLYPH<0> 2 in the halos of the reddest galaxies. This trend is a signature of the physical mechanism of color transformation.  \n<!-- image -->  \nstate (quasar mode), as a function of M ? of the central galaxy. For low-mass galaxies this ratio is effectively zero as these SMBHs have never been in the low-state. For the most massive galaxies this ratio asymptotes to GLYPH<24> 0.5, but is always less than one. A rapid transition between these two regimes occurs between 10 10 M GLYPH<12> and 10 10 : 5 M GLYPH<12> , coinciding with the stellar mass where the median color of galaxies evolves from predominantly blue to predominantly red.  \nFigure 9. The flux of galaxies across the color-mass plane with time in TNG300-1. Both panels include all redshift zero central galaxies and show evolution between z = 0 and z = 0 : 3, corresponding to 3.4 Gyr - the results are not sensitive to this time interval. Contours show the locations of the red and blue populations at redshift zero for reference. Left panel: individual arrows show the median origin M ? and (g-r) color of all the galaxies which arrive to the indicated bin by z = 0. The background grid shows the binning, with green and gray arrows on alternative rows. Right panel: the background color measures the median rate of color evolution, computed from all the galaxies in that bin at redshift zero, in magnitudes per billion years. Overlaid are streamlines (in black) showing the direction of the local vector field ( @ M ?=@ t ; @ ( g GLYPH<0> r ) =@ t ) of galaxy evolution at that point. Several evolutionary phases are evident including (horizontal) mass growth along the blue population and (vertical) reddening towards quiescence. These figures reveal the dynamic and evolving nature of the color-mass plane by showing the bulk evolution of galaxy ensembles, and do not necessarily represent the pathway of any individual system.  \n<!-- image -->  \n<!-- image -->  \n0.10  \nIn the TNG model, blackholes enter the low-state when their accretion rate relative to Eddington drops below a threshold ˙ M = ˙ M edd < GLYPH<31> ( M BH) which depends on the blackhole mass (Weinberger et al. 2017, Eqn. 5). Using the median M BH GLYPH<0> M ? relation we plot the effective threshold in the same figure (dashed gray line), as well as the actual ˙ M BH accretion rates (dashed green line). As the stellar mass of a galaxy grows above GLYPH<24> 10 10 M GLYPH<12> BHaccretion drops while the value of GLYPH<31> also rises, increasingly promoting transition into the low-state. Note that the model input GLYPH<31> was set to address the stellar mass content of massive halos and the location of the knee in the stellar mass function (Pillepich et al. 2017). We therefore conclude that the mass scale where galaxies undergo rapid color transition is set by SMBH feedback in the low-accretion state. Furthermore, we find that an evolution from blue to red color at this particular mass is unavoidably linked to the steep decline of the high-mass end of the low redshift stellar mass function.  \nFinally, in Figure 8 we return to the ratio P B = P gas of magnetic to thermal gas pressure support in the halos of galaxies. Selecting all galaxies in the narrow mass bin 10 10 : 4 < M ? /M GLYPH<12> < 10 10 : 6 we measure the median relation as a function of (g-r) color, effectively at fixed stellar mass. While blue galaxies have a typical magnetic to gas pressure ratio of order unity in the halo gas, this value increases by a factor of a few at intermediate colors, before dropping sharply by a factor of 100 for the reddest objects. That is, the relative importance of the halo magnetic field with respect to the halo thermal energy is enhanced during the color transition phase and becomes negligible after. We generally expect amplification towards equipartition and GLYPH<12> = 1, as seen in the blue population. On the other hand, with P B GLYPH<28> P gas the magnetic field in the halo is effectively controlled by the gas thermodynamics. This sharp decline arises because P B drops much faster than P gas, coinciding with the changing density and thermal energy content of the hot halo. The properties of the circumgalactic gas around a galaxy therefore encode information about the energy release by the central blackhole, the cessation of star formation, and the color transition.', '5.2 The Buildup of the Color-Mass Plane': "To untangle the origin of the z = 0 properties considered thus far, we begin to examine the time evolution of galaxies, focusing on their motion through the two-dimensional space of (g-r) color and stellar mass. In Figure 9 we show two views of the flux of galaxies across the color-mass plane from z = 0 : 3 to redshift zero, a period of roughly 3.4 Gyr. In the left panel we show, for all the galaxies in each bin at z = 0, where this ensemble of galaxies originated from. Arrows start at the median mass and color at z = 0 : 3 and end at the redshift zero location. In the right panel we show the rate of color evolution in the background color, superimposing streamlines derived from the vector field of GLYPH<1> M ? and GLYPH<1> (g-r) taken across this same time period.  \nAlong the lower edge of the z = 0 blue population, (g-r) ' 0.35, we see mass growth largely in the horizontal direction. From the right panel it is clear that this line acts as a separatrix in the population between either a positive or negative mean change of galaxy (g-r) color with time. Below this line, the population becomes more blue as stellar mass increases, while above the opposite is true, leading to the butterfly pattern in color-mass space. This line of demarcation extends all the way to M ? = 10 12 M GLYPH<12> , although the origin at high mass differs. Above a stellar mass of & 10 11 M GLYPH<12> it moves up to a (g-r) color of ' 0.7, at the bottom edge of the red population. Here, galaxies above this line redden slowly due to passive stellar evolution, while those below undergo episodes of rejuvenated star formation which temporarily pollute their integrated colors with the light from young stars. From the left panel we see these galaxies can traverse the full color range, although the sparsity of high-mass blue systems at z = 0 indicates these events are rare.  \nOne feature not immediately evident is strong passive mass growth along the passive red population. The most massive galaxies must either grow in mass through a succession of dry mergers after becoming red at lower mass, or directly join the red sequence at essentially their final mass.", '5.3 Timescale of the Color Transition and Formation of the Red Population': "To measure the evolution of individual galaxies in to, and out of, the red and blue populations, we first repeat the characterization of the color-mass plane as described in Section 4.2 over the range 0 6 z 6 1. Colors are always derived from the rest-frame spectra of galaxies; changing colors indicate underlying physical evolution. At each discrete redshift and in every stellar mass bin this yields independent fits of GLYPH<22> red ; GLYPH<22> blue ; GLYPH<27> red ; GLYPH<27> blue. At each M ? we then fit for the redshift evolution of each of these four parameters assuming a powerlaw in (1 + z ). For example, at M ? = 10 11 M GLYPH<12> (in the red population) we find GLYPH<22> red / (1 + z ) GLYPH<0> 0 : 3 and GLYPH<27> red / (1 + z ) 0 : 8 such that from z = 0 to z = 1 the location of the red peak drops by ' 0.15 mag in (g-r) and its width increases by almost a factor of two. At M ? = 10 10 M GLYPH<12> (in the blue population) we find GLYPH<22> blue / (1 + z ) GLYPH<0> 0 : 6 and GLYPH<27> blue / (1 + z ) GLYPH<0> 0 : 4 such that, over the same redshift range, the blue peak moves down by ' 0.2 mag in (g-r) while its width decreases by GLYPH<24> 30%.  \nThis quantitative characterization of the evolving blue/red galaxy distributions in color-mass space leads to a straightforward way to measure the color evolution of individual galaxies. We define two boundaries in color, C blue and C red, such that a galaxy with C < C blue is classified as blue, C > C red as red, and C blue < C < C red as an intermediate ('green valley') system. Each boundary is taken as some multiple of the width of that population away from its center location, both varying with stellar mass as well as redshift. Specifically,  \nC blue( M ?; z ) = GLYPH<22> blue( M ?; z ) + f blue GLYPH<1> GLYPH<27> blue( M ?; z ) ; C red( M ?; z ) = GLYPH<22> red( M ?; z ) GLYPH<0> f red GLYPH<1> GLYPH<27> red( M ?; z ) : (7)  \nThe parameters f blue and f red set the precise locations of the color boundaries. We take f blue = f red = 1 : 0 as a compromise between (i) the separation between the red and blue populations, and (ii) contamination of the intermediate group by non-transitional galaxies. The timescales we measure below can be modulated by these values, which is to say, the time it takes galaxies to evolve through a region of color space depends on the definition of that region.  \nFor every galaxy at z = 0 we follow the evolution of its M ? and (g-r) color back to z = 1. We locate all crossings of the evolving red and blue boundaries, if applicable. For each galaxy we then know the first time (highest redshift) it left the blue population t blue, the time it subsequently joined the red population t red ; ini, and the timescale of its color transition GLYPH<1> t green - its green valley occupation time. We also know its stellar mass upon exiting the blue population M ? ( t blue) and when entering the red population M ? ( t red ; ini). Any mass growth after joining the red sequence is given by GLYPH<1> M ?; red = M ? ( z = 0) GLYPH<0> M ? ( t red ; ini). After becoming red, a galaxy may at some later point drop out of the red population by satisfying C < C red, in which case we record the number N rejuv and time(s) t i rejuv ( i = 0 ; 1 ; 2 ; ::: ) of these 'rejuvenation' events, as well as the mass growth GLYPH<1> M ? ( t i rejuv ) and timescale(s) GLYPH<28> i rejuv of each, if the galaxy rejoins the red population thereafter.  \n8  \nFigure 10 (top panel) presents the distribution of color transition timescales GLYPH<1> t green for all galaxies with M ? > 10 9 M GLYPH<12> at z = 0. Approximately 60% of all galaxies transition between the blue and red populations in less than two billion years, with the distributions peaked at GLYPH<24> 1.5 Gyr. An extended tail reaches to much longer timescales, with GLYPH<24> 10% of galaxies taking 4 Gyr or longer to transition in color. Although broad, the distributions of GLYPH<1> t green are unimodal (roughly log-normal) and do not necessarily require the ex-  \n<!-- image -->  \nFigure 10. Top. The timescale of color transition GLYPH<1> t green measured as ( t red ; ini GLYPH<0> t blue), corresponding to the time difference between leaving the blue population and entering the red population. We show for TNG100/TNG300 the timescale PDF for all galaxies with M ? > 10 9 M GLYPH<12> (solid) and the separate distributions for z = 0 central versus satellite galaxies only (dotted and dashed, respectively). In all cases the timescale distribution is broad, approximately log-normal, with a peak at GLYPH<24> 1.5 Gyr. Bottom. Relation between GLYPH<1> t green and stellar mass when the galaxy became red M ? ( t red ; ini). At M ? . 10 10 M GLYPH<12> galaxies require GLYPH<24> 2.5 Gyr to transit the green valley. Above this mass, we find a decreasing trend with mass: galaxies which arrive into the red population with more stellar mass do so more rapidly. Shaded areas show the 10-90th percentiles for the 'All Galaxies' samples.  \n<!-- image -->  \nistence of multiple color transition pathways with distinct physical timescales. Quantitatively, in TNG100-1 the median GLYPH<1> t green across all galaxies is 1.6 Gyr with the 10-90th percentiles spanning 0.7 Gyr to 3.8 Gyr. There is only a minor difference between z = 0 central versus satellite galaxies: satellites have a slightly narrower distribution of color transition timescales, P 90 10 =[0.8 Gyr, 3.7 Gyr]. This reflects a non-trivial similarity between central and satellite quenching timescales, the latter arising only after some threshold of environmental influence.  \nIn Figure 10 (bottom panel) we show the timescale of color transition GLYPH<1> t green as a function of the 'initial red' galaxy stellar mass M ? ( t red ; ini). Low mass systems which enter the red population with M ? . 10 10 M GLYPH<12> take on average 2 GLYPH<0> 3 billion years to become red (albeit with significant variation). Galaxies arriving with progressively larger stellar masses do so quicker, taking only GLYPH<24> 1 Gyr at M ? & 10 11 M GLYPH<12> . We suggest that the increasingly common onset of  \n<!-- image -->  \nFigure 11. Left. The initial stellar mass of galaxies joining the red population, M ? ( t red ; ini). Only central galaxies at z = 0 are shown to exclude the contribution from low-mass satellite reddening. In dotted lines we show the distributions of z = 0 stellar masses of these same samples of galaxies, which are qualitatively quite similar. Right. The fraction of final z = 0 stellar mass of galaxies which originates from growth after joining the red population GLYPH<1> M ?; red = M ?; z = 0. Only central galaxies with M ?; z = 0 above the indicated minima are included to focus on the massive population. Right Inset: Correlation between this red-phase fractional mass growth and the z = 0 stellar mass of the galaxy. Shaded areas show 10-90th percentiles for the 'All' samples.  \n<!-- image -->  \nlow-state BH feedback at GLYPH<24> 10 10 : 5 M GLYPH<12> provides for the steep drop in GLYPH<1> t green for centrals evident around this mass.  \nUpon completion of a transition from blue to red over this timescale, galaxies will arrive to the red population with a range of stellar masses. Figure 11 shows the distribution of M ? ( t red ; ini), the initial M ? upon first crossing C red and entering the red population, for central galaxies (left panel). 'Initial red' galaxy masses are peaked strongly at ' 10 10 : 5 M GLYPH<12> , and it would be tempting to conclude that the vast majority of massive galaxies join the red population at its low-mass end. However, the drop-off towards large M ? also encapsulates the correspondingly steep stellar mass function. Comparison with the z = 0 stellar mass of these same populations (dashed lines) clarifies that most galaxies join the red population roughly in proportion to where they reside at z = 0, and that we should examine the differential stellar mass growth.  \nEven in the complete absence of any further in-situ star formation, a newly red galaxy may still increase its M ? down to z = 0. In Figure 11 (right panel) we show the distribution of this mass increase GLYPH<1> M ?; red normalized by final z = 0 stellar mass (main panel). We include only central galaxies with M ?; z = 0 > 10 10 : 5 M GLYPH<12> to focus on the massive population and exclude the negative tail arising from satellite mass loss. In this regime 0 < GLYPH<1> M ?; red = M ?; z = 0 < 1 and for TNG100-1 the average fraction of the final stellar mass which accumulates post-quiescence is 0 : 11, while the upper 10th percentile of galaxies exceed 0 : 41. Above 10 11 M GLYPH<12> these values increase to 0 : 25 and 0 : 53, respectively. Here, ' 18% of all galaxies acquire half or more of their present day mass on the red sequence. Asmall number of the most massive galaxies can acquire nearly all of their z = 0 stellar mass while on the red sequence - we find examples as high as 80%. The inset panel (upper right) quantifies this correlation between red-phase mass growth and total stellar mass. The observed differences between TNG100 and TNG300 primarily reflect the lower mass resolution of the latter.  \nWenote that for satellite galaxies with M ? below 10 10 : 5 M GLYPH<12> the mean GLYPH<1> M ?; red = M ?; z = 0 is GLYPH<0> 0 : 22 with the lowest 10th percentile obtaining GLYPH<0> 0 : 44. The lowest mass satellite galaxies can lose > 50% of their stellar mass after becoming red. Post-reddening stellar mass growth (or loss) is clearly a strong function of both stellar mass and environment.", '6 DISCUSSION': "So far we have presented our analysis of the color distribution of TNG galaxies without ever showing either a blue or red galaxy. We therefore turn to Figures 12 and 13 which visualize two z = 0 samples of 35 blue and 35 red galaxies, respectively, making the connection to stellar morphology.  \nStellar light composites in the optical emphasize younger stars in the blue channel and older stars in the red. Galaxies are shown face-on, projections are 60 physical kiloparsecs on a side and each panel includes the galaxy stellar mass, subhalo ID, (g-r) color, and GLYPH<20>?; rot j Jz > 0, a proxy for the rotational support (strong spheroids having GLYPH<20> . 0 : 3 and strong disks having GLYPH<20> & 0 : 6; RodriguezGomez et al. 2017). Both samples are drawn from a halo mass bin 10 12 < M halo/M GLYPH<12> < 10 12 : 2 where overlapping red and blue populations exist at similar stellar masses. Here we select blue galaxies as having (g-r) < 0 : 6, which may also include transitional color objects, while we require (g-r) > 0 : 6 for the red galaxy sample.  \nIn this roughly Milky Way mass bin the bluer galaxies are typically large, star-forming disks. Spiral features are often present, commonly traced by the blue light of younger stellar populations. Spiral arms range from tightly wound configurations, to open grand-design patterns, to flocculent and even more irregular arrangements, although we caution that the we do not expect the physical mechanisms of spiral arm formation to be realistically captured in cosmological simulations like TNG. Some systems have strong central bulges, bar-like features, clumpy disk substructures, and/or merging companions; 10 : 4 . M ? / M GLYPH<12> . 10 : 6 and 0 : 5 . GLYPH<20>? . 0 : 7 are typical values.  \nFor passive centrals this is a relatively low mass scale, and the red galaxy sample is a mix of spheroidal and disk morphologies. We see that some 'red disks' exist, as well as the more common featureless ellipticals. Many of these systems exhibit finegrained phase space structures including stellar shells and elongated stellar streams from tidally destroyed satellites, indicative of recent merger activity - both features become more common at higher M ? . Typically 10 : 5 . M ? / M GLYPH<12> . 10 : 7 and 0 : 2 . GLYPH<20>? . 0 : 4 for the spheroids, while the passive disks continue to have GLYPH<20>? & 0 : 6.  \nFigure 12. Stellar light composites of a sample of z = 0 blue galaxies, selected as having (g-r) < 0.6, and taken from the halo mass bin 10 12 < M halo/M GLYPH<12> < 10 12 : 2 . Imaged in NIRCam f200W, f115W, and F070W filters (face-on), no dust. Every panel is 60 kpc on a side.  \n<!-- image -->  \nFigure 13. Stellar light composites at z = 0 as in Figure 12 except here for red galaxies, selected as having (g-r) > 0.6 and in the same halo mass bin 10 12 < M halo/M GLYPH<12> < 10 12 : 2 , corresponding to the low-mass end of the red sequence.  \n<!-- image -->  \nThese two visual samples demonstrate the existence of multiple, distinct galaxy types spanning a wide range of stellar morphologies and structural characteristics. This diversity of the galaxy population even at fixed mass , as observed to occur in the real universe, represents a fundamental achievement (and stringent requirement) of any theoretical model for galaxy formation.", '6.1 Galaxy Color Evolution and the Formation of the Color-Mass Plane': "We seek to draw together our conclusions and thereby arrive at a coherent picture for the color distribution of low redshift galaxies and the physical mechanisms most responsible. In Figure 14 we therefore present a schematic view of the color-mass plane and the movement of galaxies therein. Gray contours show the distribution of color versus stellar mass for all galaxies at redshift zero, which divides broadly into the blue population at low stellar mass and the red population at high stellar mass. The dotted contours indicate the bulk evolution of these two groups to z GLYPH<24> 0 : 5, each becoming bluer at earlier times. With eleven colored lines we show characteristic evolutionary tracks which individual central galaxies can take through this space. They encompass all prototypical behaviors without being exhaustive.  \nThe red population continually incorporates new members, as galaxies transitioning over the timescale GLYPH<1> t green arrive with a stellar mass M ? ( t red ; ini). Figure 11 quantified the distribution of this 'initial red' mass, which is strongly peaked at GLYPH<24> 10 10 : 5 M GLYPH<12> - the purple and pink tracks which ascend almost vertically from the high-mass end of the blue population are therefore typical in this regard. Although most redshift zero red sequence galaxies, by number, join the population at its low-mass end, this is by and large a reflection of the steep stellar mass function - most z = 0 red sequence galaxies also simply reside at its low-mass end.  \nEven without any further in-situ star formation, passive red galaxies may grow further in M ? down to z = 0 by acquiring new stars through mergers (van Dokkum et al. 1999; Naab et al. 2006), with accreted stellar mass fractions of 50% or higher for galaxies exceeding 10 11 M GLYPH<12> (Rodriguez-Gomez et al. 2015; Pillepich et al. 2017). The multi-merger track (pink) and GLYPH<1> M ?; red GLYPH<29> 0 track (gold) respectively provide typical and extreme examples of this mechanism. The former undergoes two distinct dry merger events which together increase the galaxy stellar mass by a factor of two, with marginal impact on its color. The latter undergoes several dry mergers which increase its mass by a factor of four, becoming more red in the process. Figure 11 quantifies this red sequence mass growth in terms of GLYPH<1> M ?; red = M ?; z = 0, the fraction of the redshift zero stellar mass of a galaxy accumulated after its entrance into the red population - above M ? > 10 11 M GLYPH<12> galaxies on average accumulate 25% of their final mass after reddening.  \nWe naturally contrast our findings for GLYPH<1> M ?; red with the picture proposed in Faber et al. (2007). In particular, their Figure 10 presents three schematic possibilities of galaxy migration to, and growth along, the red sequence. Adopting that nomenclature, we here specifically refer to galaxies which enter the red population at its low mass end ( GLYPH<24> 10 10 : 5 M GLYPH<12> ) as 'early quenching' and those which enter the red population at higher masses as 'late quenching'. We rule out the 'early quenching + dry merging' scenario. Namely, while early quenchers are plentiful (epitomized in our purple and brown tracks), they grow the least through subsequent dry stellar accretion. The 'late quenching + no dry merging' scenario is common and important. We call these GLYPH<1> M ?; red = 0 tracks, typified by our red trajectory. We note that evolutionary paths of this shape  \n(with GLYPH<1> M ?; red ' 0) occur across the full spectrum of M ? ( t red ; ini), and at any redshift. The 'mixed, quenching + dry merging' scenario is also common and important. We label these as GLYPH<1> M ?; red > 0 tracks, typified by the pink and gold trajectories. Both mechanisms are present; therefore, neither alone can capture the manner in which the red sequence builds up from the entire galaxy population.  \nOur contribution to this theoretical picture is given in Figure 11 (inset), which quantifies that high mass galaxies experience more red-phase growth. That is, in the 'mixed' scenario of Faber et al. (2007), arrows denoting GLYPH<1> M ? along the red sequence become longer towards higher redshift zero M ? . In general, we also find that galaxies which enter the red population with more mass do so later: at 10 10 : 5 M GLYPH<12> systems transition at a mean ¯ z red ; ini GLYPH<24> 0.6 while those which enter at 10 11 : 5 M GLYPH<12> do so at ¯ z red ; ini GLYPH<24> 0.3.  \nAs we have already seen, there is significant variation in GLYPH<1> t green even at similar stellar mass. Figure 10 quantifies the distribution of this green valley occupation timescale, and we find a median value across all galaxies of 1.6 Gyr with the 10-90th percentiles spanning 0.7 Gyr to 3.8 Gyr. A large fraction of the population is therefore consistent with passive reddening after a short quenching timescale, consistent with observational demographics of poststarburst galaxies (e.g. Wild et al. 2016) and passive evolution times using single SSPs (e.g. Trayford et al. 2016). At the same time, the presence of existing older stellar population could induce a more rapid blue to red transition, while dust may also be an important factor in color transition timescales (e.g. Gutcke et al. 2017). As the symbols for each track in our schematic figure are placed at z = f 0 ; 0 : 1 ; 0 : 2 ; 0 : 3 ; 0 : 4 ; 0 : 5 ; 0 : 7 ; 1 : 0 ; 1 : 5 ; 2 : 0 g , they are spaced roughly 1 Gyr apart. As an example, while the brown and pink paths begin at the same mass and color at z = 2, they traverse the green valley over different timescales, the multi-merger case requiring GLYPH<24> 2 Gyr less. Trayford et al. (2016) also computed color transition timescales from cosmological simulations, finding a median of GLYPH<1> t green ' 2 : 0 GLYPH<0> 2 : 5 Gyr for 10 10 < M ? / M GLYPH<12> < 10 11 with a tail towards longer values and no strong mass dependence. Despite being qualitatively similar, in detail this result differs from our findings in Figure 10 by up to a factor of two, as well as in the decreasing dependence of GLYPH<1> t green on stellar mass. This is undoubtedly due to a combination of different physical models and so quenching mechanisms between the two simulations, compounded with different definitions of the color transition region and the definition of the transition between blue and red galaxies. We find that the most massive central galaxies of the present epoch have spent the least amount of time in the green valley, partially reflecting the merging loci of the blue and red populations towards the highest masses.  \nAlthough we have not presented comprehensive results on the process of rejuvenated star formation in the red population, we make a few brief observations from TNG100-1. Once red, the majority of galaxies do not ever leave the red population. Across the entire resolved mass range, 94% of redshift zero galaxies have N rejuv = 0. If we restrict our attention to massive, passive galaxies by considering M ?; z = 0 / M GLYPH<12> > 10 11 M GLYPH<12> we find that 90% of galaxies have never undergone a rejuvenation event, while 10% have rejuvenated once and GLYPH<24> 1% have gone through two or more such events. These 'firefly' galaxies (light orange line) continually relight their star formation and so jump along the red sequence by increasing their stellar mass in intermittent flashes. We find a maximum of N rejuv = 5, although such extreme cases represent a negligible fraction of the galaxy population. The distribution of rejuvenation timescales t i rejuv has a similar shape as the GLYPH<1> t green distribution, with a slightly lower normalization. The median time between exiting and re-entering the red population is 1.1 Gyr with a percentile  \nFigure 14. Schematic diagram of galaxy evolution across the color-mass plane. The distribution of all z = 0 galaxies is shown by the gray contours, while the dotted lines indicate the bulk evolution of the red and blue populations to z ' 0 : 5, both moving downward in (g-r). Eleven colored lines show the characteristic evolutionary pathways of individual, central galaxies , from at most z = 2 to the present day. Although drawn from actual tracks in TNG100-1, we intend them as prototypical demonstrations of the various trajectories of central galaxies towards z = 0 and in particular the ways they can ascend to the red population (see text).  \n<!-- image -->  \nspread P 90 10 =[0.5 Gyr, 2.7 Gyr]. The duty cycle of rejuvenation in massive red galaxies should correspond to observable residual star formation in such objects (Fogarty et al. 2015; McDonald et al. 2016). Some massive galaxies can attribute almost all of their latetime stellar growth to a single large rejuvenation event. The light red track provides one example - despite the fact that this galaxy enters the red population already by z = 1 it resumes star formation and does not finally quench until GLYPH<24> 5 Gyr later. Note that the total P GLYPH<1> M ? ( t i rejuv ) is included by definition in the GLYPH<1> M ?; red values we present throughout, and that we could miss rejuvenation events due to the finite time spacing of our simulation snapshots or due to the precise definition of the red color boundary.", '6.2 The Connection to Galaxy Morphology': "We have also not addressed the question of morphological transformation, an intriguing aspect of galaxy formation which is intimately connected to color. Specifically, we have not attempted to discern if color transition and the quenching of star formation is always accompanied by morphological change, and if so, if these occur co-temporal or if one precedes the other. From our blue/red samples in Figures 12 and 13 we already see by eye the relation between morphology and color in this halo mass bin (see also Snyder et al. 2015). We add a few brief comments here. At fixed M ? = 10 10 : 5 M GLYPH<12> we find a trend of morphology as measured through GLYPH<20>? and (g-r) color, such that redder galaxies are on average less rotationally supported. We find no correlation between  \nGLYPH<20>? ( z = 0) and GLYPH<1> t green. On the other hand, GLYPH<20>? exhibits a suggestive scaling with red-phase mass growth, such that GLYPH<20>? ' 0 : 4 for GLYPH<1> M ?; red = M ?; z = 0 GLYPH<24> 0 while GLYPH<20>? . 0 : 2 for GLYPH<1> M ?; red = M ?; z = 0 & 0 : 25. That is, galaxies which acquire more mass on the red sequence are driven towards preferentially spheroidal morphologies; the underlying physical cause is undoubtedly merger activity, which lowers the former while increasing the latter. A similar trend exists with GLYPH<1> M green, hinting that mass growth in the green valley is associated with disturbed, non-steady-state star formation. This may indirectly agree with the observationally based hypothesis of Powell et al. (2017) that a lack of green spheroids results from the rapidity of their quenching via mergers. However, the lack of correlation between z = 0 morphology and GLYPH<1> t green does not appear to support the idea that early-type (or late-type) galaxies should arrive to the red population over characteristically shorter (or longer) transition timescales due to either near instantaneous quenching or slow secular gas exhaustion, respectively (Schawinski et al. 2014). A more detailed study along these lines is clearly needed. The interwoven pathways of stellar reddening and structural change, together with the origin of red disks, blue spheroids, and the connection with central stellar bulges, all promise to be a fruitful future direction of investigation.", '6.3 Final Considerations': "Throughout this paper we have intentionally focused on the massive red population and its emergence from previously star-forming  \ngalaxies. That is, we have neglected the properties and evolution of red satellites, satellite quenching, and environment. The upper left corner of Figure 14 is populated by low-mass red galaxies, which are almost entirely satellites in dense environments. From observations we expect strong environmental dependence (including one-halo conformity) in color related properties such as colorsplit luminosity functions or stellar mass functions (Baldry et al. 2006), red fraction (Peng et al. 2010), color-split clustering (Wang &White 2012; Farrow et al. 2015; Zu & Mandelbaum 2017; LawSmith & Eisenstein 2017), and quiescent fraction (Cucciati et al. 2017). We also expect different physical processes such as tidal and ram-pressure stripping to play important roles in satellite galaxy transformation, not only in color but also in star formation rate, gas content, and morphology (e.g. Ann et al. 2008; Sales et al. 2015). We leave for the future an investigation of the environmental dependency of galaxy color in TNG.  \nFinally, we have not specifically addressed the balance between stellar feedback - that is, our supernova driven galactic winds - and blackhole feedback (for an interesting view see Bower et al. 2017). We have shown that the low-state BH feedback mode is responsible for inducing color transformation in the TNG model and setting its characteristics. Interestingly, we find only a moderate correlation of blackhole mass with galaxy color at our canonical stellar transition mass of 10 10 : 5 M GLYPH<12> . Here the reddest galaxies have blackholes which are ' 0.4 dex more massive than their bluest counterparts. This is because blackhole mass growth in the TNG model is dominated by accretion in the high state (quasar mode). Therefore, even though transition to the low-state quenches galaxies and transforms galaxies from blue to red, it cannot leave a strong impression in the relic BH total mass. We do not see the same magnitude signature of excess SMBH mass as found in Trayford et al. (2016) (their Figure 3). However, in the integrated ratio of mass growth between the low and high states, the integrated ratio of energy injection (as shown in our Figure 7), or the total energy injection in the low-state, we do see a clear marker of blackhole induced quenching. For example, at M ? = 10 10 : 5 M GLYPH<12> blue galaxies have a total integrated E injected ; low ' 10 58 erg, while red galaxies have injected a hundred times more energy in the low accretion-state, with median E injected ; low ' 10 60 erg. There is a pronounced increase in these BH related indicators at a color threshold of (g-r) = 0.55 mag.  \nIn Figure 14 the remaining three lowest mass evolutionary tracks - cyan, orange, and blue - are all below this color threshold, and all still members of the blue population at z = 0. Each has grown in stellar mass through quasi-equilibrium oscillations of star formation activity, which drives their optical color periodically up and down. Indeed the self-regulated nature of star formation should establish the color width of the blue population. Here we identify the origin of the features seen at low M ? in Figure 9. Because of these oscillations, galaxies below a center line will have generally arrived from redder colors, while those above will have arrived from bluer colors. This establishes the separatrix in the streamline field and the butterfly pattern of the flux of galaxies through this region of the color-mass plane.", '7 SUMMARY AND CONCLUSIONS': "In this paper we have introduced the TNG100 and TNG300 simulations, the first two efforts of the IllustrisTNG project, a series of 'next generation' cosmological magnetohydrodynamical simulations using the updated, comprehensive TNG galaxy formation model. Leveraging their unique combination of volume, resolution,  \nand physical fidelity, we model the optical colors of the simulated galaxy populations, incorporating a resolved treatment of attenuation from interstellar dust. 3 We first compare to observational data from SDSS and then, bolstered by the excellent agreement, study the relationships between color and other galaxy properties, the process of color transformation, and the physical mechanisms driving the buildup of the color-mass plane. Our main results are:  \n- GLYPH<15> The simulated (g-r) colors of TNG galaxies at low redshift are in excellent agreement with a quantitative comparison to observational data from SDSS at z < 0 : 1. The color bimodality for intermediate mass galaxies is demonstrably improved with respect to the original Illustris simulation. We recover the locations in color of both the red and blue populations as a function of M ? , as well as the relative strength between the red and blue distributions, the location of the color minimum between the two, and the location of the maximal bimodality in stellar mass. We show the same level of quantitative agreement exists in other colors including (u-r), (u-i), while slightly less favorable comparison to the dust-sensitive (r-i) may indicate shortcomings in our relatively simple dust modeling.\n- GLYPH<15> We use a double gaussian mixture model to characterize the red and blue populations and apply the same MCMC fitting procedure to both the simulations and observations. The loci of the red and blue distributions agree across 10 9 < M ? / M GLYPH<12> < 10 12 to within 0.05 mag, while the scatter is consistent within 0.02 mag for blue galaxies and 0.01 mag for red galaxies. This foundational agreement allows us to identify a second order discrepancy of possible interest: the slope of the red sequence may be too flat, an effect likely driven by the finite aperture choice together with non-zero color gradients in extended stellar light.\n- GLYPH<15> As the two primary ingredients of colors beyond reddening, we directly compare the stellar ages and metallicities in TNG versus observed relations. The luminosity-weighted, fiber aperture restricted stellar ages are in remarkable agreement. The trend of Z ? versus M ? appears too shallow, implying that the simulated stellar metallicities are too high for low-mass galaxies. However, we show that this discrepancy plausibly arises only due to an overly simplistic comparison, and that a more robust analysis based on consistent fits to mock SDSS fiber spectra is required (Nelson et al. in prep).\n- GLYPH<15> Relating redshift zero (g-r) color to other galaxy properties, we find as expected that red galaxies are quiescent, old, and gas depleted. Compared to otherwise increasing trends with M ? , they also have lower gas-phase metallicities and lower interstellar magnetic field strengths. We find a non-monotonic correlation, at fixed stellar mass, between galaxy color and the GLYPH<12> ratio of magnetic to thermal pressure in the gaseous halo, such that the P B = P gas ratio peaks for green valley systems and drops sharply thereafter for red galaxies. We point out an explicit prediction of the TNG model: at M ? ' 10 10 : 5 M GLYPH<12> blue galaxies should have interstellar magnetic field strengths of 10 GLYPH<0> 30 GLYPH<22> G, while red galaxies at the same mass should have much lower j B ISM j . 1 GLYPH<22> G.\n- GLYPH<15> We attribute the primary driver of galaxy color transition in the TNG model to supermassive blackhole feedback in its lowaccretion state. We demonstrate that the ratio of the cumulative energy injection by blackholes in the high (quasar) versus low (kinetic wind) state drops sharply around a characteristic M ? GLYPH<24> 10 10 : 5 M GLYPH<12> , a value set by our M BH-dependent accretion rate threshold. Having  \nbeen fixed during model development in order to match the stellar mass content of halos and in particular the knee of the stellar mass function, we conclude that color transition at this mass scale is inseparably linked to a realistic high-mass galaxy population.  \n- GLYPH<15> Across all z = 0 galaxies, we measure a median color transition timescale GLYPH<1> t green of GLYPH<24> 1.6 Gyr, with a 10-90th percentile spread of 0 : 7 GLYPH<0> 3 : 8 Gyr. This distribution of green valley occupation times is unimodal, approximately log-normal, and does not show direct evidence for more than one characteristic timescale. We find a relation between GLYPH<1> t green and stellar mass, such that galaxies which enter the red population at larger M ? do so with a correspondingly more rapid color transition and shorter green valley occupation.\n- GLYPH<15> We measure the amount of 'red-phase' stellar mass growth, between the time a galaxy enters the red population and z = 0. Massive galaxies with M ? > 10 11 M GLYPH<12> accumulate on average 25% of their mass after reddening; at the same time, the top ' 18% of galaxies accumulate half or more of their z = 0 stellar mass while on the red sequence. The amount of GLYPH<1> M ?; red is found to increase with increasing redshift zero M ? . We therefore advocate for an updated 'mixed' scenario for the formation of the red population whereby galaxies become red across a wide range of stellar mass, some fraction subsequently having GLYPH<1> M ?; red = 0 while the remainder grow by an amount GLYPH<1> M ?; red > 0 which increases for more massive systems.  \nTogether, the TNG100 and TNG300 simulations have a dynamic range of 10 7 per spatial dimension. This enables us to connect the internal characteristics of galaxies - notably, their colors - to details of their host halos, as well as to the statistical properties of the global galaxy population within a full cosmological context. Our present investigation of the low redshift galaxy color distribution demonstrates the predictive and explanatory power of these simulations and their impressive utility as a theoretical tool for understanding the process of galaxy formation.", 'ACKNOWLEDGEMENTS': 'DN would like to thank Charlie Conroy for development of the FSPS code and for helpful comments on the dust methods, Ben Johnson for development of the PROSPECTOR code and assistance in its use, Vicente Rodriguez-Gomez for developing and allowing us to use the SUBLINK merger tree code, Arjen van der Wel for ideas on the spectral modeling, and Simon White for useful comments and suggestions.  \nVS, RW, and RP acknowledge support through the European Research Council under ERC-StG grant EXAGAL-308037 and would like to thank the Klaus Tschira Foundation. SG through the Flatiron Institute is supported by the Simons Foundation. PT and SG are supported by NASA through Hubble Fellowship grants HST-HF-51384.001-A and HST-HF2-51341.001-A , respectively, awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. MV acknowledges support from a MIT RSC award, the Alfred P. Sloan Foundation, and by NASA ATP grant NNX17AG29G. JN acknowledges support from NSF AARF award AST-1402480.  \nIn the execution of the primary simulations presented herein, the authors gratefully acknowledge the Gauss Centre for Supercomputing (GCS) for providing computing time for the GCS Large-Scale Projects GCS-ILLU (2014) and GCS-DWAR (2016) on the GCS share of the supercomputer Hazel Hen at the High Performance Computing Center Stuttgart (HLRS). GCS is the alliance  \nof the three national supercomputing centres HLRS (Universitat Stuttgart), JSC (Forschungszentrum Julich), and LRZ (Bayerische Akademie der Wissenschaften), funded by the German Federal Ministry of Education and Research (BMBF) and the German State Ministries for Research of Baden-Wurttemberg (MWK), Bayern (StMWFK) and Nordrhein-Westfalen (MIWF). Additional simulations were carried out on the Hydra and Draco supercomputers at the Max Planck Computing and Data Facility (MPCDF, formerly known as RZG), as well as on the Stampede supercomputer at the Texas Advanced Computing Center through the XSEDE project AST140063. Some additional computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University.', 'APPENDIX A: THE DUST UNDER THE RUG': "In this appendix we show a few additional figures which explore the differences between the dust models and the resolution convergence of our results.  \nFirst, Figure A1 shows the output of our color modeling on five individual galaxy examples spanning a range of stellar masses and star formation rates, as shown (colored lines). For each, the (g-r) colors derived with Model A and Model B are shown as single circles and squares, respectively. As these models have no viewing angle dependence, each galaxy has one unique color. The y-axis locations of these points are therefore arbitrary. PDFs of the Model C colors are also shown, for N side = 4 (dotted; 192 projections per galaxy) and N side = 8 (solid; 768 projections per galaxy), the view-angle dependent color distributions converging with sufficient sampling of the sphere. Note that in the main text we have used N side = 1 throughout, and taken either all 12 resultant colors per galaxy or a single random color of these twelve, depending on the specific application.  \nFigure A2 compares the different dust models on the full galaxy populations, in three stellar mass bins from TNG100-1. The unresolved dust correction of Model B is significant for low-mass (blue) galaxies, while the impact of Model C is relatively small in comparison. The modeled attenuation is also important in setting the precise location of the red peak for more massive galaxies, which is shifted blueward in the case of no dust. We note that Model C aims to account for attenuation by 'resolved' gas structures in the simulation, and high density clumpy substructure in the ISM is far from resolved in the TNG simulations, implying that in higher resolution simulations we may be better able to model small-scale dust effects than in the current work. Including all twelve Model C colors for each galaxy or including only one at random has little impact on these global distributions, although this may not hold if we preferentially selected colors based on edge-on or face-on projections. We also include an as of yet unexplored Model D, a variation on Model C which considers an alternative geometrical assumption for the dust distribution along each line of sight. Specifically, it implements the 'uniform scattering model' of Calzetti et al. (1994) instead of the internal dust model. In this case, the calculation of GLYPH<28>GLYPH<21> is identical, and the attenuated luminosity is given by L obs( GLYPH<21> ) = L i( GLYPH<21> ) = e GLYPH<0> GLYPH<28>GLYPH<21> instead of by Equation (2). The differences between Model C and Model D are small, and reflect only a minor change in the positive contribution of dust scattering to the emergent luminosity as a function of how the dust column is distributed along each line of sight through the gas distribution.  \nFinally, we comment on resolution convergence. As a reminder, we do not adjust in any way our model parameters as a function of resolution (for details see Pillepich et al. 2017). This is in contrast to the approach adopted in the EAGLE project where feedback mechanisms are specifically adjusted to counterbalance any (un)physical change with resolution and thereby to achieve a roughly, or at least more-so, constant result as a function of numerical resolution (Schaye et al. 2015). Different properties of our simulated galaxies have different convergence properties, the most fundamental of which are discussed in Appendix A of Pillepich et al. (2017). In Figure 2 we add to that discussion by showing the TNG100-1 and TNG300-1 predicted stellar ages and metallicities. In doing so, we are effectively looking at the resolution convergence of these quantities, since TNG300-1 is 'one step lower' in resolution (a factor of two in spatial resolution, eight in mass) and therefore equivalent to the intermediate TNG100-2 run. We find that the stellar age and metallicity trends have somewhat different  \nFigure A1. Impact of viewing angle variation in the calculated (g-r) color values for five example systems spanning a range of stellar masses. For each galaxy, Model A and Model B results are shown by single circles or squares, respectively. The PDF of colors for Model C is then shown for each using high N side parameters of 4 (dashed) or 8 (solid lines), corresponding to 192 and 768 projections per galaxy, respectively. The most massive galaxy has a predominantly spherical distribution of both stars and gas and thus exhibits little color variation due to projection-dependent dust effects. The less massive galaxies show asymmetric color distributions skewed to the red, although the total widths are typically less than 0.1 mag for all but the most dusty systems.  \n<!-- image -->  \nconvergence behaviors. Given the low-statistics at the high-mass end for TNG100-1, the stellar ages appear to be well converged, better than most galaxy properties. Stellar metallicities increase, at fixed M ? , with increasing resolution, by roughly 0 : 05 GLYPH<0> 0 : 1 dex between the TNG300-1 (intermediate) and TNG100-1 (high) resolution levels. We caution that, as in any case when showing M ? on the x-axis, any intrinsic shift in a specific galaxy property is convolved with the fact that TNG stellar masses are not perfectly converged, and so galaxies will also shift along the x-axis, typically towards higher M ? at fixed halo mass with increasing resolution.  \nFigure A3 presents the resolution convergence of the galaxy color PDFs directly, in the same three stellar mass bins as before. Keeping in mind that galaxies will shift between bins due to changing stellar masses, the general result is that with progressively lower resolution we lose low-mass blue galaxies. As a result, there are too many low mass red galaxies and the separation between the two populations becomes less well defined, and there is also an increase in the number of high-mass blue galaxies.  \nAnother culprit in the apparent lack of resolution convergence for the bluest systems arises directly from the finite mass resolution of the simulations and the stochastic sampling of the physical process of star formation. Namely, a low-mass galaxy on the star formation main sequence with a SFR of GLYPH<24> 1 M GLYPH<12> /yr will on average create only one 10 6 M GLYPH<12> stellar particle every million years at TNG100-1 resolution, and only one 10 7 M GLYPH<12> stellar particle every ten million years at TNG100-2 resolution. This discreteness implies that at lower resolutions there may be no very young star particles present even though a galaxy is actively star forming. Motivated by this exact problem, Trayford et al. (2017) (using the EAGLE simulation) decided to implement a re-sampling technique in postprocessing to effectively generate the expected populations of very young star particles. Indeed, even with perfect resolution convergence of the physical properties of a galaxy, its color would still become redder with worsening mass resolution due to such a stochas-  \nTable A1. Table of physical and numerical parameters for the three resolution levels of each of the two simulations presented in this paper, TNG100 and TNG300. The physical parameters are: the box volume, the box side-length, the initial number of gas cells, dark matter particles, and Monte Carlo tracer particles. The target baryon mass, the dark matter particle mass, the z =0 Plummer equivalent gravitational softening of the collisionless component, the same value in comoving units, the maximum softening applied to blackholes, and the minimum comoving value of the adaptive gas gravitational softenings. Additional characterizations of the gas resolution, measured at redshift zero: the minimum physical gas cell radius, the median gas cell radius, the mean radius of star-forming gas cells, the mean hydrogen number density of star-forming gas cells, and the maximum hydrogen gas density. The numerical parameters are: the number of high time frequency subbox volumes, the number of snapshots in each subbox, the total number of timesteps to z =0, the total run time including substructure identification in millions of CPU core hours, and the number of compute cores used.Table A2. Table of physical and numerical parameters as in Table A1, except for the six dark-matter only analog runs.  \nTables A1 and A2 give the most relevant physical and numerical parameters for the baryonic and dark-matter only runs, respectively. Each includes the highest resolution realization of the two boxes presented herein, TNG100-1 and TNG300-1, as well as the two lower resolution versions of each. Among other values we include the mean radius of star-forming gas cells ¯ r gas ; SF in physical parsecs, and the mean number density of star-forming gas cells  \ntic sampling. This effect provides a caution for accurate colors in lower resolution simulations, and likely drives much of the resolution trend seen in this Figure.  \n¯ n gas ; SF in particles per cm 3 . These provide a rough idea of the spatial scales and gas densities at which the physical process of star formation takes place in our numerical simulations, namely, GLYPH<24> 350 pc and GLYPH<24> 1 cm GLYPH<0> 3 in the case of TNG100-1.  \n5  \nFigure A2. Building up the dust model used in this paper, in the middle three stellar mass bins from 10 9 : 5 < M ? /Msun < 10 11 : 0 . First we show Model A (brown lines), which includes the stellar population synthesis modeling with no dust treatment. Next, Model B (purple lines) includes only the unresolved presence of dense birth clouds. Model C includes the resolved dust column calculations with one random line of sight per galaxy (red lines) or all twelve projections shown simultaneously (green lines). Model D uses a different dust geometry assumption (see text; blue lines).  \n<!-- image -->  \n5  \nFigure A3. The resolution convergence of the color distributions (Model C), shown in the central three stellar mass bins explored in this work, ranging from 10 9 : 5 < M ? /Msun < 10 11 : 0 . Although the resolution convergence is fair, it is far from perfect. In general, with increasing mass and spatial resolution, galaxies are bluer, primarily as a result of larger stellar metallicities and larger stellar masses. Note however that galaxies also shift between bins due to changing stellar masses.  \n<!-- image -->"}