The TNG50-SKIRT Atlas: Post-processing methodology and first data release

Galaxy morphology is a powerful diagnostic to assess the realism of cosmological hydrodynamical simulations. Determining the morphology of simulated galaxies requires the generation of synthetic images through 3D radiative transfer post-processing that properly accounts for different stellar populations and interstellar dust attenuation. We use the SKIRT code to generate the TNG50-SKIRT Atlas, a synthetic UV to near-infrared broadband image atlas for a complete stellar-mass selected sample of 1154 galaxies extracted from the TNG50 cosmological simulation at z = 0. The images have a high spatial resolution (100 pc) and a wide field of view (160 kpc). In addition to the dust-obscured images, we also release dust-free images and physical parameter property maps with matching characteristics. As a sanity check and preview application we discuss the UVJ diagram of the galaxy sample. We investigate the effect of dust attenuation on the UVJ diagram and find that it affects both the star-forming and the quiescent galaxy populations. The quiescent galaxy region is polluted by younger and star-forming highly inclined galaxies, while dust attenuation induces a separation in inclination of the star-forming galaxy population, with low-inclination galaxies remaining at the blue side of the diagram and high-inclination galaxies systematically moving towards the red side. This image atlas can be used for a variety of other applications, including galaxy morphology studies and the investigation of local scaling relations. We publicly release the images and parameter maps, and we invite the community to use them.


Introduction
Cosmological hydrodynamical simulations are powerful tools to investigate the origin and evolution of galaxies.These simulations use numerical methods to emulate the behaviour of gas, stars, black holes, and dark matter in a virtual universe.They take into account the effects of gravity, hydrodynamics, star formation, feedback from supernovae and active galactic nuclei, and other physical processes that are known or expected to influence the formation and evolution of galaxies.For recent reviews, we refer to Somerville & Davé (2015), Vogelsberger et al. (2020a), and Crain & van de Voort (2023).
In general, our confidence in the conclusions drawn from cosmological hydrodynamical simulations increases with the level of agreement between simulations and observations.Vice versa, detailed comparisons are required to calibrate the subgrid physics recipes in the simulations and to guide the improvements for the next-generation simulations.Only a decade ago, the overall agreement between hydrodynamical simulations and observations of galaxies was relatively poor (e.g.Navarro & Benz 1991;Navarro & Steinmetz 1997, 2000;Sommer-Larsen et al. 1999;Scannapieco et al. 2012).In recent years the agreement has improved significantly, mainly due to an increase in resolution and better implementations of subgrid physics recipes.Large-volume cosmological hydrodynamical simulations such as Illustris and IllustrisTNG (Vogelsberger et al. 2014;Pillepich et al. 2018a), EAGLE (Schaye et al. 2015;Crain et al. 2015), Magneticum (Dolag et al. 2016), Horizon-AGN and NewHorizon (Dubois et al. 2016(Dubois et al. , 2021)), or SIMBA (Davé et al. 2019) succeed in reproducing many observed global galaxy properties to a fair degree, including stellar and gas mass functions, the colour bimodality, the star-forming main sequence, the massmetallicity relation, and many other scaling relations (e.g.De Rossi et al. 2017;Bottrell et al. 2017;Torrey et al. 2019;Davé et al. 2020;Akins et al. 2022;de Graaff et al. 2022;Zenocratti et al. 2022, to name just a few examples) Galaxies are no structureless point sources but complex ecosystems.The comparison of simulated galaxies to observations on local rather than global scales provides a powerful and meaningful test for cosmological hydrodynamical simulations.With many large imaging and integral-field spectroscopic surveys just finished, ongoing, or about to start, such as the SAMI Galaxy Survey (Bryant et al. 2015;Croom et al. 2021), MaNGa (Bundy et al. 2015;Abdurro'uf et al. 2022a), the DESI Legacy Imaging Surveys (Dey et al. 2019), WEAVE-StePS (Iovino et al. 2023), or the Euclid Wide Survey (Euclid Collaboration: Scaramella et al. 2022), the required observational data are nowadays available.This offers the prospect of testing the simulation's fidelity in much more detail, and to identify limitations due to both particle resolution or incorrect physics implemented in the models.
Galaxy morphology is a particularly interesting characteristic to assess the realism of cosmological hydrodynamical simulations.Morphology has served as the main basis for galaxy classification since the early days of extragalactic astronomy (Sandage 2005, and references therein).However, galaxy morphology is important beyond classification alone.It is a powerful diagnostic of the physical processes driving galaxy evolution (e.g.bars, mergers, spiral structure, tidal features).Morphological parameters are found to correlate with several fundamental galaxy properties such as stellar mass, colour, star formation history, merger history, and local galactic environment (Dressler 1980;Gómez et al. 2003;Kauffmann et al. 2003;Conselice 2003Conselice , 2014;;Lotz et al. 2008;van der Wel 2008;Blanton & Moustakas 2009;Bluck et al. 2014).
In order to determine the morphology of simulated galaxies, we have to apply a forward modelling or post-processing approach on them.Since we know the 3D position and characteristics of all stellar particles in a cosmological hydrodynamical simulation, this seems a fairly straightforward projection operation.However, the generation of realistic synthetic images can be quite involved, especially due to the effects of absorption and scattering by dust grains in the interstellar medium.Indeed, dust attenuates a third to half of all the starlight in typical spiral galaxies in the local Universe (Popescu & Tuffs 2002;Viaene et al. 2016;Bianchi et al. 2018), and detailed radiative calculations have shown that the effects of dust attenuation are often complex (Witt et al. 1992;Byun et al. 1994;Pierini et al. 2004;Möllenhoff et al. 2006;Gadotti et al. 2010).The only way to generate synthetic images that properly account for dust absorption and scattering in a realistic geometry is by 3D radiative transfer modelling, a demanding and challenging numerical problem (Steinacker et al. 2013).
In the past few years, there have been several efforts to generate synthetic images for large sets of galaxies extracted from cosmological hydrodynamical simulations.Torrey et al. (2015) used the SUNRISE radiative transfer code (Jonsson 2006) to generate a broadband image atlas for about 7000 z = 0 galaxies extracted from the Illustris simulation (without including dust in their models).Trayford et al. (2015) used the SKIRT code (Camps & Baes 2015, 2020) to produce an image atlas for more than 30,000 galaxies from the EAGLE simulation at z = 0.1, fully accounting for dust attenuation.Similarly, Rodriguez-Gomez et al. ( 2019) generated realistic synthetic Pan-STARRS-like images for more than 12,000 galaxies from the TNG100 simulation using SKIRT.The power of these large-volume-simulation-based image databases is that they cover the entire galaxy population, from massive ellipticals to less massive star-forming galaxies.
The drawback is that the physical resolution of the underlying simulations, and hence also of the mock images, is limited to about 1 kpc, a limiting factor in both realism and range of applications.
This disadvantage can be addressed by shifting from largevolume to higher-resolution zoom-in simulations, such as FIRE (Hopkins et al. 2014), NIHAO (Wang et al. 2015), APOSTLE (Sawala et al. 2016), Latte (Wetzel et al. 2016), Auriga (Grand et al. 2017), RomulusC (Tremmel et al. 2019), or ARTEMIS (Font et al. 2020).Kapoor et al. (2021) and Camps et al. (2022) generated a suite of synthetic high-resolution images for the present-day galaxies from the Auriga and ARTEMIS simulations, respectively, with a combined image database that covers the entire UV-submm wavelength range.Faucher et al. (2023) also used SKIRT to generate UV-submm images for a set of present-day galaxies from the NIHAO simulation.The main limitation of these studies is that they only cover a relatively small number of galaxies (30 from Auriga, 45 from ARTEMIS, 65 from NIHAO) and the limited dynamical range in galaxy properties.
The TNG50 simulation (Pillepich et al. 2019;Nelson et al. 2019b) offers the prospect of combining the advantages of both approaches as it combines a high spatial and mass resolution typical for zoom-in simulations with a large number and diversity of galaxies typical for large-volume simulations.Very recently, Guzmán-Ortega et al. (2023) presented a database of simulated KiDS-like gri images for more than 5000 galaxies from the TNG50 simulation at z = 0.034 to investigate the connection between galaxy mergers and optical morphology in the local Universe.Their approach is similar to Rodriguez-Gomez et al. ( 2019) and fully accounts for dust attenuation using SKIRT.At the same time, several teams have generated synthetic optical data cubes for samples of galaxies extracted from the TNG50 simulation, with different levels of refinement in the way dust effects are approximated (Bottrell & Hani 2022;Nanni et al. 2022Nanni et al. , 2023;;Sarmiento et al. 2023).
In this paper we present and release the TNG50-SKIRT Atlas (hereafter TSA), a synthetic image database for 1154 galaxies at z = 0 extracted from the TNG50 simulation.We cover a wide range of galaxy properties, we generate mock images in 18 broadband filters in the UV to near-infrared (NIR) wavelength range with high spatial resolution and a wide field of view, while fully accounting for dust absorption and scattering.Moreover, we also generate images that are free of dust attenuation, as well as synthetic maps of intrinsic physical parameters such that connections between the observed and the intrinsic properties can easily be made.
This paper is organised as follows.In Sect. 2 we briefly present the TNG50 cosmological hydrodynamical simulation and the SKIRT radiative transfer code, and we discuss our methodology to generate synthetic images for the TNG50 galaxies.In Sect. 3 we present the multi-wavelength image database.In Sect. 4 we discuss the UVJ diagram of the galaxy sample as an illustration of the range of applications that is possible with these data.We present a summary and our conclusions in Sect. 5.

The TNG50 simulation
The TNG50 simulation (Pillepich et al. 2019;Nelson et al. 2019b) is the highest-resolution version of the IllustrisTNG cosmological magneto-hydrodynamical simulations suite (Marinacci et al. 2018;Naiman et al. 2018;Nelson et al. 2018;  1. Relation between stellar mass, half-mass radius, and SFR for the TNG50 galaxies in our sample.In the bottom panel, the dashed line indicates the separation between quiescent and star-forming galaxies (sSFR = 10 −11 yr −1 ).Galaxies with no ongoing star-formation are plotted at SFR = 2 × 10 −4 M ⊙ yr −1 .Pillepich et al. 2018b;Springel et al. 2018).It follows the evolution of a cubic volume of 51.7 comoving Mpc on the side, with cosmological parameters based on the Planck 2015 results, namely Ω m = 0.3089, Ω b = 0.0486, Ω Λ = 0.6911, and H 0 = 67.74km s −1 Mpc −1 (Planck Collaboration et al. 2016).TNG50 uses the moving-mesh code AREPO (Springel 2010) as its hydrodynamics solver.The galaxy formation model (Weinberger et al. 2017;Pillepich et al. 2018a) in the TNG50 simulation is identical for all three TNG simulations.The physical processes accounted for in the simulation include gas cooling and heating, stochastic star formation, stellar evolution, chemical enrichment of the ISM, feedback from supernovae, seeding and growth of supermassive black holes, AGN feedback, and magnetic fields.TNG50 reaches a baryonic mass resolution of 8.5 × 10 4 M ⊙ ; the average cell size in the star-forming regions of galaxies is 70-140 pc.For a full description of the TNG50 simulation we refer to Pillepich et al. (2019) and Nelson et al. (2019b).

Sample selection
To build the TSA, we selected all TNG50 galaxies at z = 0 with total stellar mass between 10 9.8 and 10 12 M ⊙ from the TNG50 public database (Nelson et al. 2019a), resulting in a sample of 1154 individual galaxies.The characteristics of the sample are illustrated in Fig. 1, which shows the correlations between stellar mass, star formation rate (SFR), and half-mass radius.As evident from the bottom panel, most galaxies in the sample are starforming galaxies lying in the main sequence, whereas there is also a population of more quiescent galaxies.Adopting a specific star formation rate (sSFR) of sSFR = 10 −11 yr −1 as the boundary between star-forming and quiescent galaxies (e.g.Brinchmann et al. 2004;Fontanot et al. 2009;Donnari et al. 2019;Paspaliaris et al. 2023), our sample contains 869 star-forming and 285 quiescent galaxies.

The SKIRT radiative transfer code
SKIRT (Camps & Baes 2015, 2020) is a three-dimensional Monte Carlo radiative transfer code.Originally set up as a dust radiative transfer code and designed to model the effects of dust in galaxies (Baes et al. 2003(Baes et al. , 2011)), it has transformed into a more generic Monte Carlo radiative transfer tool.The physical ingredients within SKIRT include absorption, multiple scattering, and stochastic emission (Camps et al. 2015); polarisation caused by scattering or by emission from aligned nonspherical dust grains (Peest et al. 2017;Vandenbroucke et al. 2021); Lyα resonant line scattering (Camps et al. 2021); absorption and emission at rotational or electronic transitions for selected ions, atoms and molecules (Gebek et al. 2023;Matsumoto et al. 2023); photo-absorption, fluorescence, and scattering at X-ray wavelengths (Vander Meulen et al. 2023).SKIRT is equipped with a library of flexible input models, routines to import the output from various kinds of hydrodynamical simulations, and a selection of advanced spatial grids for discretising the medium (Baes & Camps 2015;Saftly et al. 2013Saftly et al. , 2014;;Camps et al. 2013).Many Monte Carlo radiative transfer optimisation mechanisms and a hybrid parallelisation strategy are implemented to maximise the efficiency of the code (Baes 2008;Steinacker et al. 2013;Baes et al. 2016Baes et al. , 2022;;Verstocken et al. 2017).

The SKIRT setup to generate high-resolution images
Each SKIRT simulation is completely determined by specifying the characteristics of the primary radiation sources, the transfer medium, a suite of instruments to capture the emerging radiation field, probes to measure other aspects of the simulation, and a number of technical or numerical simulation parameters.The recipe for the SKIRT simulations we adopted in this work largely follows the prescriptions outlined by Trčka et al. (2022), which were inspired by earlier works by Camps et al. (2016Camps et al. ( , 2018Camps et al. ( , 2022) ) and Kapoor et al. (2021).We only give a brief overview, focusing on the particular aspects of the current work, and refer to these papers for more details.

Primary radiation sources
For each galaxy, the primary radiation sources are the stellar particles that belong to the galaxy, which were extracted from the TNG50 database.To each stellar particle older than 10 Myr, we assigned a simple stellar population (SSP) SED from the Bruzual & Charlot (2003) SSP family with a Chabrier ( 2003) initial mass function.The mass, metallicity, and age of the SSP are directly inherited from the TNG50 particle data.Each particle was given a smoothing length corresponding to the distance to the 32nd nearest neighbour, with a maximum smoothing length of 800 pc.
Stellar particles younger than 10 Myr are assumed to be still partly embedded in their birth cloud.They were assigned an SED from a library of Hii region templates based on the MAP-PINGS III library (Groves et al. 2008).Each template in this library is characterised by five free parameters, two of which (metallicity and SFR) can be directly obtained from the particle data.Since the specific choice of the ISM pressure value hardly affects the broadband SED shape (Groves et al. 2008, Fig. 4), we used a fixed value, P = 1.38 × 10 −12 Pa.For the compactness, which essentially determines the dust temperature and thus the shape of the far-infrared spectrum, we sampled a value from a lognormal distribution (Kapoor et al. 2021;Trčka et al. 2022) with parameters calibrated on the dust temperature distribution in observed and simulated star forming regions (Utomo et al. 2019;Kannan et al. 2020).Finally, the PDR covering factor was set to f PDF = e −t/τ with t the particle age and τ the PDR clearing timescale, a free parameter in the post-processing recipe.The value τ = 3 Myr was determined by Trčka et al. (2022) from a calibration of the integrated TNG50 fluxes against observational data from the DustPedia nearby galaxy sample (Davies et al. 2017;Clark et al. 2018).

The dusty medium
The properties of the transfer medium, in our case the dusty interstellar medium, are based on the characteristics of the Voronoi gas cells in the hydrodynamical simulation, which were again downloaded from the TNG50 database.The density of the dust at every location in the SKIRT simulation volume is based on the assumption that a fixed fraction f dust of the metals in the ISM gas is locked up in dust grains, that is with Z gas the metallicity and ρ gas the density of the gas.To determine which gas cells qualify as ISM gas cells, we applied the prescription by Torrey et al. (2012Torrey et al. ( , 2019) ) that separates the hot circumgalactic medium from the cooler ISM gas.In this framework, gas is considered as belonging to the ISM if its temperature T gas satisfies the condition log T gas K < 6 + 0.25 log ρ gas 10 10 h 2 kpc −3 . (2) The density, metallicity and temperature of each gas cell were directly taken from the cell data.The only free parameter left was the dust-to-metal fraction f dust , for which Trčka et al. (2022) determined f dust = 0.2 as the best value in combination with the PDR clearing timescale of τ = 3 Myr (see previous section).This value is comparable to observational estimates based on different nearby galaxy samples (e.g.De Vis et al. 2019;Galliano et al. 2021;Zabel et al. 2021).
We assumed a fixed dust grain model at every location in the galaxy.As our dust grain model, we used the diffuse ISM THEMIS model (Jones et al. 2017).It consists of a distribution of carbonaceous and silicate grains, with optical properties based on laboratory data where possible, and can reproduce many observed properties of the dust in the Milky Way, including the UV to infrared extinction curve, extinction correlations, the thermal dust emission spectrum, and red and blue luminescence.For the present paper we focused on the UV to near-infrared wavelength range, and we did not include diffuse dust emission in the radiative transfer procedure (ISM dust emission is only important The SKIRT code requires a grid structure on which the medium is discretised.Since the TNG50 simulation uses the AREPO moving mesh technique to solve the hydrodynamics, we could have opted to use the native Voronoi grid for the SKIRT post-processing as well.SKIRT is equipped with a mechanism to efficiently traverse photon packets through a Voronoi grid (Camps et al. 2013), and this approach has been adopted to post-process moving-mesh hydrodynamical simulations (e.g.Rodriguez-Gomez et al. 2019;Schulz et al. 2020;Camps et al. 2021;Popping et al. 2022;Guzmán-Ortega et al. 2023).However, like in many other works (e.g.Vogelsberger et al. 2020b;Kapoor et al. 2021;Trčka et al. 2022;Hsu et al. 2023;Barrientos Acevedo et al. 2023) we chose to re-grid the dust density distribution onto a hierarchical octree grid to boost the speed of the radiative transfer process.We used up to 12 levels of refinement to ensure that we fully resolve the dust density distribution.

Instruments and probes
We defined a set of broadband imagers as instruments in our SKIRT simulation.Contrary to Kapoor et al. (2021) and Camps et al. (2022), who positioned the observers at specific sight lines relative to orientation of each galaxy (face-on, edge-on, and at specific inclination angles), we aimed at arbitrary viewing points.We used fixed viewing points relative to the simulation box, which has no connection to the orientation of each individual galaxy.In order to limit the computation time and the data volume, we settled on N obs = 5 viewing positions per galaxy.These positions were spread on the unit sphere in optimal arrangement, that is, so as to maximise the angular distance between them.This problem is generally known in geometry as the Tammes problem (Tammes 1930).Exact solutions for low values of N obs are available in the literature (e.g.Fejes Tóth 1943; Schütte & van der Waerden 1953;Erber & Hockney 1991), as well as approximate numerical solutions for larger values of N obs (e.g.Kottwitz 1991;Hardin et al. 1994;Steinacker et al. 1996).For the case N obs = 5, two of the observer positions are antipodal, and this semi-redundancy provides an opportunity to test the accuracy of any analysis method on the images.The details of the observer positions are provided in Table 1.
For the instruments, a pixel scale of 100 pc provides a nice match to the spatial resolution of the TNG50 simulation (Pillepich et al. 2019;Nelson et al. 2019b).The actual size of the detectors, or equivalently, the field-of-view, does not significantly affect the simulation run time, but it does impact the data volume.We chose detectors with 1600×1600 pixels, correspond- ing to a field of view of 160 kpc on the side, in order to cover the outer regions of the most extended galaxies.SKIRT offers the opportunity to define broadband data cubes, that is, synthetic instruments that contain multiple broadband images automatically convolved with the correct transmission curves (Camps & Baes 2020, Sect. 4.5).We generated images in the GALEX FUV-and NUV-bands, the Johnson UBVRIbands, the LSST ugrizy-bands, the 2MASS JHK s -bands, and the WISE W1-and W2-bands.Together, these bands cover the wavelength range between 0.1 to 5 µm, the domain of dominance of stellar emission.
Apart from instruments that record the radiation escaping from the simulation volume, SKIRT offers the opportunity to install probes, which sample numerical or physical quantities internal to the simulated model.Probes thus provide relevant diagnostics on the simulation setup and offer an opportunity to investigate properties of the simulated model that could never be directly observed from the outside, such as for example the radiation field (Camps & Baes 2020, Sect. 3.5).Specifically for this project, a set of probes was implemented in SKIRT to generate intrinsic physical parameter maps for an arbitrary observer's position.These maps are generated by projecting the 3D physical fields such as the stellar mass density on the observer's plane of the sky.Concretely, we generated stellar mass surface density maps, stellar-mass-weighted metallicity and age maps, and dust mass surface density maps for each of the five observers, with the same pixel scale and field-of-view as the instruments described above.
In a similar way, the emissivity of the stellar particles within a given broadband can be projected on an observer's plane of the sky, resulting in dust-free broadband images.The main advantage of this approach is that it does not involve any Monte Carlo radiative transfer calculation, and hence generates no Monte Carlo noise.We generated a corresponding dust-free image for each observer and each broadband filter.

Numerical parameter settings
Apart from the definition of the sources, the dust medium, the instruments, and the probes, each SKIRT simulation requires a number of numerical parameters related to the different wavelength grids used internally or to the different Monte Carlo optimisation techniques.There was no reason to deviate from the default values, which are optimised to assure the best performance (see also Trčka et al. 2022).
A crucial parameter that directly affects the quality of the output data and the simulation run time is the number of launched photon packets.The optimal value was determined empirically using a limited number of galaxies of different stellar masses and sizes.We re-simulated the same galaxies with increasing number of photon packets until we found convergence on a pixel by pixel basis.Convergence was defined using the reliability statistics introduced in SKIRT by Camps & Baes (2018, 2020), based on work in the field of nuclear particle transport simulations (X-5 Monte Carlo Team 2003).In order to obtain convergence in the optical r-band out to at least twice the halfmass radius of the galaxies, we found that at least 10 9 photon packets are required.For the FUV-and NUV-bands, the corresponding convergence area is typically smaller due to the lower intrinsic emissivity and the larger dust attenuation.We fixed the number of photon packets to 10 9 for all galaxies.For a typical TNG50 galaxy the corresponding SKIRT run time on a dedicated 64-core machine was 3 to 4 hours, though there was a large variety depending on the number of stellar particles and the number of grid cells in the hierarchical octree grid.

Atlas characteristics and availability
For each of the 1154 TNG50 galaxies in our selection and for each of the 5 observing positions, the data set consists of dustobscured and dust-free images in 18 broadbands, and 4 physical parameter maps.In total, this adds up to exactly 200 images per galaxy, or 230 800 images in total.Each image or map is stored as an individual 1600 × 1600 pixel FITS file.The data and user guidelines on how to easily access and use them are available on the SKIRT website.1More details on the TSA characteristics are listed in Table 2.
The broadband images have units of MJy sr −1 and the parameter maps have natural units (M ⊙ pc −2 for the stellar and dust surface density maps, dimensionless units for the mean stellar metallicity maps, and Gyr for the mean stellar age maps).We note that the images are not convolved with a PSF and are noisefree, except for the presence of the Monte Carlo noise in the dust-obscured images.Users who wish to generate synthetic images matching a particular instrument or survey can convolve the images with the appropriate PSF and add background noise (see e.g.Rodriguez-Gomez et al. 2019;de Graaff et al. 2022;Guzmán-Ortega et al. 2023).

Example galaxies
Fig. 2 shows RGB images for 36 randomly selected galaxies, using the LSST u-, g-, and z-band images for observer position O1.The galaxies are sorted by increasing stellar mass.This gallery illustrates the wide variety in optical morphology in the galaxy population: some galaxies are regular and almost featureless, others have a clear bulge-disc structure, and some have a disturbed morphology.In a number of galaxies, the effects of dust attenuation are clearly visible.Also the optical colours and SFRs of the galaxies vary widely, both globally and on spatially resolved scales.
Figs. 3 and 4 show a more detailed view on two handpicked galaxies from the data set.Fig. 3 shows TNG000008, a spiral galaxy with stellar mass M ⋆ = 3.73 × 10 10 M ⊙ and SFR = 3.98 M ⊙ yr −1 , here seen at an intermediate inclination.
The stellar mass distribution is characterised by a small but dense central component of metal-rich stars, whereas the disc is on average less metal-rich.The mass-averaged age of the stellar population is around 7 Gyr for the central bulge component and only about 3 Gyr for the populations in the spiral arms.The dust distribution is strongly concentrated in the central regions and the spiral arms and has a filamentary morphology.The broadband images show the signature as expected based on the physical parameter maps.The u-band image is dominated by the young stars and has a very clumpy appearance.It also shows clear signs of dust attenuation at the locations of the most prominent dust mass surface density enhancements.Moving to longer wavelengths, the images become less clumpy, the apparent signatures of attenuation gradually disappear, and the bulge gradually becomes more prominent.
Fig. 4 shows a very different galaxy, TNG096764, a massive early-type galaxy with M ⋆ = 1.05×10 11 M ⊙ and no ongoing star formation.The stellar mass surface density map shows a smooth distribution without obvious substructure, except a small bar in the central regions (viewed from this single observer's position, the central feature could be either a bar or an edge-on disc, but the combination of the different viewing angles shows that it is a bar rather than a disc).The stellar population shows a strong metallicity gradient, and the bar seen in the stellar mass surface density map stands out by its high metallicity.The same structure is also seen in the age map, where the feature stands out due to a slightly younger mean stellar age.On average, the stellar populations are relatively old and the age gradients modest.Furthermore, the galaxy displays a ring-like structure with a radius of about 5 kpc that is slightly younger (∼6 Gyr) than the overall stellar population (∼9 Gyr).Also interesting is the dust mass surface density map, which shows a conspicuous, albeit very low surface density, spiral-arm-like feature.Since there is no ongoing star formation and hardly any dust attenuation, the morphology of all broadband images is similar.It is very smooth without much substructure, except the inner bar and a faint shelllike structure in the bottom right quadrant.The u-band image of this galaxy is smoother than the u-band image of the spiral galaxy shown in Fig. 3.

Comparison of TNG50, TNG100, and EAGLE
The UVJ diagram was first presented by Wuyts et al. (2007) and Williams et al. (2009) as a way to distinguish between quiescent and star-forming galaxies.A single optical colour is not enough to separate these two classes as different combinations of dust attenuation and SFR can give rise to similar optical colours.However, the combination of an optical and an optical-NIR colour proves successful at making that distinction.The UVJ method has been widely applied to select quiescent galaxies from observed galaxy samples (Whitaker et al. 2013;Straatman et al. 2014;Barro et al. 2014;Papovich et al. 2015Papovich et al. , 2018;;Fang et al. 2018;Tan et al. 2022;Valentino et al. 2023) and from cosmological hydrodynamical simulations (Davé et al. 2017;Trayford et al. 2017;Donnari et al. 2019;Akins et al. 2022;Kurinchi-Vendhan et al. 2023).
The first panel of Fig. 5 shows the UVJ diagram for the galaxies in our TSA.We calculated the rest-frame U-, V-, and J-magnitudes by integrating the surface brightness of our images over the entire field-of-view.The colour scale in the UVJ diagram represents the sSFR, which is taken directly from the TNG50 database.The other panels show equivalent UVJ diagrams based on other simulated galaxy data sets, all at z = 0 and within the same stellar mass range we used.The second panel contains the same set of galaxies, but used the 'random' viewing position flux values generated by Trčka et al. (2022).The third panel is based on SKIRT-generated fluxes presented by Camps et al. (2018) for the flagship EAGLE simulation.The final panel corresponds to the TNG100 simulation, for which Gebek et al. ( 2024) generated broadband fluxes using the same methodology as applied by Trčka et al. (2022) for the TNG50 simulation.
In all panels, the data show a clear correlation between the UVJ colours and the sSFR.The dotted line, taken from Donnari et al. (2019), separates the quiescent galaxy population from the star-forming galaxies.For all four cases, the separation line corresponds roughly to a fixed sSFR ∼ 10 −10.7 yr −1 .The general agreement between the different simulations is a good sanity check for the SKIRT calculations executed in this work.
The differences between the different panels are at least as interesting.The main difference is the coverage of the UVJ diagram, which can primarily be attributed to the sample size.The TNG100 fluxes generated by Gebek et al. (2024) reach lower U − V colours for a fixed V − J colour than the other diagrams.The galaxies populating these blue regions at the bottom of the UVJ diagram are among the most actively star-forming galaxies.A remarkable difference is that the distribution of star-forming galaxies in our new TNG50 UVJ diagram extends towards very red colours in the diagonal direction, beyond the vertical separation line at V − J = 1.6, whereas it ends rather abruptly before this line in the other three panels (apart from some scattered galaxies in the TNG100 diagram).The galaxies populating this part of the UVJ diagram are very actively star-forming with sSFR > 10 −10 yr −1 .At the same time they have very red colours, which implies that they must be heavily attenuated.It is, at first sight, surprising that these galaxies are lacking in the UVJ diagram based on the fluxes calculated by Trčka et al. (2022) since the galaxy samples used for the two TNG50 diagrams are exactly the same.The reason is that we consider four different random viewing positions for each galaxy, compared to a single random flux for the Trčka et al. (2022) catalogue.As one can visually infer from Fig. 2, the level of dust attenuation in the TNG50 galaxies can differ substantially depending on the viewing angle, which can move galaxies across the UVJ plane for different

Dust attenuation and physical properties in the UVJ diagram
While the UVJ diagram is widely used to discriminate between star-forming and quiescent galaxies, its origin and reliability is being investigated in more detail (Belli et al. 2019;Leja et al. 2019;Díaz-García et al. 2019;Wu et al. 2020;Nagaraj et al. 2022;Antwi-Danso et al. 2023).More specifically, it has been investigated to which degree the UVJ colours can be used to constrain other physical properties of galaxies, and how dust attenuation affects the UVJ colours of galaxies.The general effect of dust attenuation is that galaxies move upwards in the UVJ diagram, roughly parallel to the diagonal separation line.The exact direction depends on the shape of the attenuation curve, which can vary significantly among galaxies (e.g.Salmon et al. 2016;Leja et al. 2017;Salim et al. 2018;Narayanan et al. 2018;Qin et al. 2022;Zhang et al. 2023).The attenuation in our models is calculated in full 3D and takes into account absorption and scattering, which can sometimes lead to counterintuitive effects (Witt et al. 1992;Byun et al. 1994).Moreover, our models have more complex, and hopefully more realistic, star formation histories than the parametric models often assumed.
In the top row panels of Fig. 6, we show the UVJ diagram based on the dust-obscured images for the galaxies in our sample colour-coded by six different physical parameters.The bottom row contains the same information, but now for UVJ colours based on the dust-free images.The correlation between UVJ colours and sSFR (see also Fig. 5) is obvious for both the dustobscured and dust-free diagrams, but it is not the only systematic trend.Also mean stellar age, mean metallicity, V-band mass-to-light ratio, and dust-to-stellar mass ratio show a clear trend with the UVJ colours.These results are well in agreement with the results obtained by Leja et al. (2019).By means of a Bayesian inference method applied to synthetic SED models, they demonstrated that the mass-to-light ratio is well constrained by UVJ colours alone, whereas the trends with age and metallicity are induced by galaxy scaling relations.
Looking at the quiescent galaxies region of the second panel on the top row, one can see that the mean stellar ages of the galaxies systematically increase when moving upward in the diagonal direction.To investigate this in more detail, we used the rotated coordinate system on the UVJ diagram introduced by Fang et al. (2018).In this system, the rotated axes S Q and C Q , defined as run parallel and perpendicular to the boundary of the quiescent box, respectively.In Fig. 7 we show the relation between the colour index parallel to the separation line and the logarithm of the mean stellar age of each galaxy in the quiescent region.We recovered a reasonably strong linear correlation (R 2 = 0.44), which means that the ages of quiescent galaxies can in principle be estimated from UVJ colours alone.This supports the conclusions by Belli et al. (2019), who derived a similar linear relation between the logarithmic median stellar age and the S Q index.
Comparing the panels on the top row with the corresponding ones on the bottom row, a number of interesting effects of dust  attenuation can be discerned.One can immediately note the difference in coverage of the UVJ plane, which is most obvious for the star-forming galaxy population.The dust-free star-forming galaxies predominantly occupy a region with 0.5 ≲ V − J ≲ 1.1 which moves into the quiescent galaxies region as soon as U − V ≳ 1.7.The dust-obscured star-forming galaxies, on the other hand, occupy a much more extended region that runs parallel to the separation line up to V − J > 1.5.The galaxies in this new territory have high sSFR values, and so they moved from the bottom left corner in the dust-free UVJ diagram.Interestingly, a large number of these galaxies are oriented almost edge-on, in line with the findings by Patel et al. (2012).This segregation by inclination is clearly visible when comparing the panels in the fifth column.While the galaxies in the bottom panel have a uniform distribution in inclination with a median value of about 60 • in every pixel of the parameter space, the distribution is clearly separated in the top panel, in particular for the star-forming galaxies: low-inclination galaxies remain at the blue side of the diagram, whereas high-inclination galaxies move towards the red side.At the extreme end of the star-forming distribution we find edge-on galaxies with high sSFR values and very high dust-to-stellar-mass ratios.
There is also a difference in the dust-free and the dustobscured UVJ diagram for the quiescent galaxy population.In the dust-free diagram it forms a well-defined sequence out to V − J ≈ 1.5, populated exclusively by galaxies with very low sSFR values, old ages, high metallicities, large mass-to-light ratios, and very small dust-to-stellar-mass ratios.When dust atten-Article number, page 9 of 13 A&A proofs: manuscript no.TSA uation is turned on, however, the population extends to redder colours, but most importantly, the region in the UVJ space below the original tight sequence is populated, across the boundary line into the star-forming galaxies regimes.The galaxies in this new region are, again, oriented almost edge-on and heavily attenuated.We can quantify the nature of these additional galaxies by calculating some global statistics of the two different populations (as defined by the dotted line region in the plots) with and without dust attenuation.In the dust-free case, the quiescent galaxy region contains just 36% of the total stellar mass budget, and the mean sSFR and stellar age are ⟨sSFR⟩ = 10 −11.52 yr −1 and ⟨t⟩ = 8.84 Gyr.In the dust-obscured case, the pollution by star-forming galaxies increases the fraction of the stellar mass budget of the quiescent galaxy region to 48%.The mean sSFR increases by 0.4 dex to ⟨sSFR⟩ = 10 −11.12 yr −1 , and the mean stellar age decreases to ⟨t⟩ = 8.55 Gyr.In summary: dust attenuation pollutes the quiescent galaxy region with younger and more actively star-forming highly inclined galaxies.

Summary
The ambition of this work was to generate, present, and publicly release a synthetic UV-NIR broadband image atlas for a complete stellar-mass selected sample of 1154 galaxies extracted from the z = 0 snapshot of the TNG50 cosmological simulation (Pillepich et al. 2019;Nelson et al. 2019b).The images were generated with the SKIRT radiative transfer code (Camps & Baes 2015, 2020) and account for different stellar populations and absorption and scattering by interstellar dust in a realistic 3D setting.
For each galaxy, we generated a suite of 100 pc resolution images in 18 broadband filters, for five different observer positions.In addition to the dust-obscured images, we also released synthetic images without dust attenuation, and stellar mass surface density, mean stellar age, mean stellar metallicity, and dust mass surface density maps.While other teams have already released synthetic image datasets for cosmological hydrodynamical simulations (e.g.Torrey et al. 2015;Trayford et al. 2015;Rodriguez-Gomez et al. 2019;Kapoor et al. 2021;Camps et al. 2022;Guzmán-Ortega et al. 2023), we argue that the present atlas is unique in its kind due the realism of both the underlying simulation and the radiative transfer treatment, the large sample size, the high spatial resolution, the number of filters, and the combination of matching images and physical parameter maps.
As a sanity check and a first application of our image atlas, we investigated the UVJ diagram.Comparing our UVJ diagram with the one based on the fluxes generated by Trčka et al. (2022) for the same TNG50 galaxy sample, and with the EAGLE and TNG100 UVJ diagrams we found excellent agreement in terms of the relation between sSFR and UVJ colours.The diagrams also show some interesting differences, mainly in the coverage of the UVJ diagram.In particular, the distribution of star-forming galaxies in our new TNG50 UVJ diagram extends towards very red colours in the diagonal direction, whereas the other diagrams lack these very actively star-forming and heavily obscured galaxies.These differences can be interpreted as a result of different sample size and the difference in spatial resolution.
We also investigated the trends of galaxy physical parameters in the UVJ diagram.We found that, apart from the strong correlation with sSFR, the UVJ colours also show systematic trends with mean stellar age, mean stellar metallicity, V-band mass-to-light ratio, and dust-to-stellar-mass ratio.We found a reasonably strong positive correlation between the mean stellar age and the UVJ colours for the quiescent galaxy population, which suggests that the ages of quiescent galaxies can be well constrained by UVJ colours alone (Belli et al. 2019).Finally, we investigated the effect of dust attenuation on the distribution of the galaxy population in the UVJ diagram.As expected, dust attenuation spreads the galaxy population towards redder colours parallel to the separation line between the quiescent and starforming galaxy populations.In the dust-free UVJ diagram the quiescent galaxy population forms a well-defined sequence populated exclusively by galaxies with very low sSFR values and old ages, but dust attenuation pollutes this quiescent galaxy region with younger and more actively star-forming galaxies.Dust attenuation generates a clear separation in inclination of the starforming galaxies: low-inclination galaxies remain at the blue side of the diagram, whereas high-inclination galaxies move towards the red side.The reddest star-forming galaxies are edgeon, dusty galaxies with high sSFR values.

Possible applications
The prime goal of this paper was to present and publicly release the TSA.We demonstrated its usefulness by investigating the UVJ diagram.We hope that this image atlas can be used for many more applications.
A prime application is the connection between the morphology of galaxies, their fundamental physical properties, and the environment in which they reside (e.g.Conselice 2003Conselice , 2014;;Blanton & Moustakas 2009;Holwerda 2021).The availability of multi-band dust-obscured and dust-free images allows for a systematic investigation of the wavelength dependence of galaxy morphology (Muñoz-Mateos et al. 2009;Kelvin et al. 2012;Vulcani et al. 2014;Baes et al. 2020;Nersesian et al. 2023;Martorano et al. 2023) and the effects of dust attenuation on photometric and morphological parameters (Gadotti et al. 2010;Pastrav et al. 2013a,b;Savchenko et al. 2023).
In a companion paper to this release paper (Baes et al. 2024) we used the TSA to investigate the wavelength dependence of the effective radii of galaxies.In the near future we plan to employ single or multiple-component Sérsic fitting and nonparametric morphological indices to quantify morphology.All global physical properties (also including dark matter properties), the intrinsic particle and cell data, and the entire history for all galaxies can be readily accessed through the TNG public database (Nelson et al. 2019a), which allows for a thorough investigation on what drives galaxy morphology.
In the past few years, several well-known global galaxy scaling relations have been investigated on local, sub-kpc scales as well.Examples include the local dust scaling relations (Viaene et al. 2014;Casasola et al. 2022), the spatially resolved star-forming main sequence (Cano-Díaz et al. 2016;Hsieh et al. 2017;Morselli et al. 2020;Pessa et al. 2022;Abdurro'uf et al. 2022b;Baker et al. 2022), the resolved stellar mass gas metallicity relation (Rosales-Ortega et al. 2012;Barrera-Ballesteros et al. 2016;Gao et al. 2018;Boardman et al. 2023;Baker et al. 2023), or the resolved stellar mass stellar metallicity relation (González Delgado et al. 2014;Zibetti et al. 2020;Neumann et al. 2021;Zibetti & Gallazzi 2022;Pessa et al. 2023).Our image and physical parameter map database can be used to investigate the universality of these local scaling relations and to determine the physical scales at which they potentially break down.
We plan to address several of the questions raised above.However, we also publicly release these data and warmly invite the community to use them in any way they see fit.

Caveats and future work
While we believe that the current image atlas is sufficiently rich and realistic to allow for a range of applications, we are aware that it also has its caveats and limitations.A first important aspect is that our data are built on simulated galaxies extracted from the TNG50 cosmological hydrodynamical simulation.While this simulation is generally considered as one of the most powerful large-volume simulations, it comes with its own caveats and limitations.One of them is that the TNG model was designed and calibrated at the resolution of the original Illustris simulation, while the TNG50 simulation has a roughly 16 times better mass resolution.This improved resolution results in somewhat larger galaxy masses and SFRs (Pillepich et al. 2018a(Pillepich et al. ,b, 2019;;Donnari et al. 2019;Trčka et al. 2022).
A second limitation of our atlas is that it is limited to the UV-NIR range, that is, the range where stellar emission dominates the SED.This limitation is not inherent to SKIRT code: the code has been used to generate synthetic SEDs and images for galaxies that cover the entire UV-mm wavelength range.The main reason is computational: to simulate the UV-NIR range only dust attenuation and no dust emission is required.SKIRT simulations with dust emission are computationally much more demanding (up to an order of magnitude, depending on the details of the simulation) compared to attenuation-only simulations.
In previous post-processing work of cosmological hydrodynamical simulations, we encountered difficulties to reproduce the UV and MIR fluxes and colours of observed galaxies: we typically found insufficient UV attenuation and too much emission in the MIR (Baes et al. 2019;Trčka et al. 2020Trčka et al. , 2022;;Kapoor et al. 2021;Camps et al. 2022).We argued that the MAP-PINGS III templates we use for the young stellar particles are at least partly responsible for these discrepancies.Kapoor et al. (2023) recently generated a new template library, called TOD-DLERS, to be used in SKIRT with the aim of addressing this problem.The first tests of this new library are very promising (Kapoor et al., in prep.).Our ambition is to rerun our image library with this new emission library in the near future, extending the range to mm wavelengths.This new emission library is expected to not affect the optical images, but will probably improve the UV images.We advise users to take this caveat into account when they use our data.
Looking forward, we see this UV-NIR broadband image atlas as an intermediate step in an effort to generate increasingly more realistic synthetic data products for simulated galaxies.As discussed above, our next ambition is to extend this image database to the UV-mm wavelength range, incorporating the new TODDLERS library.The step beyond that could be a transition from broadband imaging to full spectral resolution.Such an effort could be similar to the works of Bottrell & Hani (2022), Nanni et al. (2022, 2023), and Sarmiento et al. (2023), but with a completely self-consistent treatment of dust attenuation and emission, and ideally covering the entire UV-mm wavelength range.Several intermediate steps towards realising this ambition have already been taken or are under development: the TOD-DLERS templates have full spectral resolution and contain a detailed treatment of nebular emission lines, SKIRT has full support for gas and stellar kinematics (Camps & Baes 2020;Barrientos Acevedo et al. 2023), and we are currently working on a sub-grid model for the multi-phase ISM, inspired by Olsen et al. (2021) and Ramos Padilla et al. (2021Padilla et al. ( , 2023)).
Finally, applying the SKIRT radiative transfer postprocessing recipe for different redshift snapshots would open up the possibility to directly investigate and test the cosmic evolution of galaxy morphology and many of the scaling relations mentioned above.
Fig.1.Relation between stellar mass, half-mass radius, and SFR for the TNG50 galaxies in our sample.In the bottom panel, the dashed line indicates the separation between quiescent and star-forming galaxies (sSFR = 10 −11 yr −1 ).Galaxies with no ongoing star-formation are plotted at SFR = 2 × 10 −4 M ⊙ yr −1 .

Fig. 3 .
Fig. 3. Representative images and physical parameter maps for the spiral galaxy TNG000008 from observer position O2 (i = 47.6 • ).Top row: stellar mass surface density, mean stellar metallicity, mean stellar age, and dust mass surface density.Bottom row: u-, g-, i-, and K s -band images.All images and parameter maps zoom into the central 30 × 30 kpc region.

Fig. 5 .
Fig.5.The UVJ plane of the TNG50, EAGLE, and TNG100 cosmological hydrodynamical simulations.Each panel contains only galaxies at z = 0 within the same stellar mass range (10 9.8 M ⊙ < M ⋆ < 10 12 M ⊙ ).All colours are calculated from SKIRT-generated fluxes or images that take dust attenuation into account.The colour scale represents the median sSFR of all galaxies with UVJ colours within each pixel.The dotted line in each panel indicates the separation between the quiescent and star-forming galaxy populations and is taken fromDonnari et al. (2019).

Fig. 6 .Fig. 7 .
Fig. 6.Physical properties in the UVJ plane, with (top row) and without (bottom row) dust attenuation.The different columns correspond to different physical parameters, indicated at the top of the column.The colour scale represents the median value of the physical parameter within each bin in UVJ colours.

Table 1 .
Polar angle (θ) and azimuth (ϕ) of the observer positions with respect to the simulation volume.

Table 2 .
Main characteristics of the TNG50-SKIRT Atlas.