A redshift-dependent IRX-$\beta$ dust attenuation relation for TNG50 galaxies

We study the relation between the UV-slope, $\beta$, and the ratio between the IR and UV luminosities of galaxies, IRX, using TNG50, the latest installment of the IllustrisTNG galaxy formation simulations. We select 7280 TNG50 star-forming galaxies with stellar mass above $10^9M_\odot$ at selected redshifts, $0 \leq z \leq 4$, and perform radiative transfer calculations with SKIRT to model the effects of ISM dust on the simulated stellar light. We adopt a Milky Way (MW) type dust and a fixed dust-to-metal ratio of 0.3 throughout and find that TNG50 star-forming main-sequence galaxies agree with the empirically-derived reference IRX-$\beta$ relations at $z \lesssim 1$. There appears to be a redshift-dependent systematic offset with respect to the reference relations, with the TNG50 IRX-$\beta$ median relation steepening towards higher redshifts. This is to first order driven by variations in intrinsic UV-slopes due to different star-formation histories of galaxies selected at different cosmic epochs. Once differences in the intrinsic UV-slope are corrected for, TNG50 galaxies exhibit a scatter of 0.3 dex in IRX at fixed $\beta$ at all studied redshifts. These galaxy-to-galaxy variations are correlated with intrinsic galaxy properties, such that at fixed UV-slope galaxies with higher star--formation rates, star--formation efficiencies, gas metallicities, and stellar masses exhibit larger IRX values. We additionally demonstrate a degeneracy between stellar age and dust type. The combination of young stellar population age with a SMC type of dust for high-redshift ($z=4$) TNG50 galaxies can result in an IRX-$\beta$ relation similar to what is found for low-redshift galaxies with MW dust. This hampers the use of the IRX-$\beta$ relation as a proxy for dust type. We provide a redshift-dependent fitting function for the IRX-$\beta$ relation with MW dust based on our numerical models.

z ≈ 2, declining again as we go towards higher redshifts (see the review by Madau & Dickinson 2014). Measurements of star formation rates (SFRs) at z > 3 are often based on the ultraviolet (UV) emission as a tracer for the abundance of young stars in a galaxy, and therefore as a tracer for SFR (Kennicutt & Evans 2012). The attenuation of stellar light by dust in the interstellar medium (ISM) is a major hindrance at these wavelengths. UV-emission is being scattered and absorbed, changing the measured shape and luminosity of the UV spectrum. The absorbed UV-light is being re-emitted in the infrared (IR) regime, and ideally we can account for the absorbed emission by measuring both the UV spectrum and the IR spectrum of a galaxy (e.g. Reddy et al. 2012). At high redshifts z > 3, reliable IR data is often missing, which is why alternatives for accounting for UV dust attenuation have been explored.
For local UV-bright starburst galaxies, an empirically derived tight relation between the rest-frame UV continuum slope β (where we assume a power law: f λ ∝ λ β ) and the ratio between IR and UV luminosity, called the infrared excess (IRX) has been found (Meurer et al. 1999). This IRX-β relation provides a connection between the shape of the UV continuum (the "color" of the galaxy) and the amount of dust attenuation. If this relation proved to hold for various galaxy types at different redshifts, it could be used as a reliable tool to account for dust attenuation where only UV data is available and IR data is missing (e.g. as in Bouwens et al. 2011, Oesch et al. 2014, McLeod et al. 2015. Recent studies have investigated the IRX-β relation of various galaxy types and at diverse redshifts. Some adjustments to the Meurer relation have been proposed to correct for biases in aperture size and sample selection (Overzier et al. 2010, Takeuchi et al. 2012, Casey et al. 2014. These corrected relations seem to agree qualitatively at z = 0 for galaxies on the star forming main sequence (SFMS), and also agree with observations of main sequence galaxies at z < 4 (e.g. Heinis et al. 2013, Bouwens et al. 2016, Álvarez-Márquez et al. 2016, Bourne et al. 2017, Fudamoto et al. 2017, Koprowski et al. 2018, McLure et al. 2018, Reddy et al. 2012, Fudamoto et al. 2019. However, allowing for a more complete sample of galaxies including IR-bright starbursts and submillimeter galaxies (SMGs) and going towards higher redshifts increases the scatter in the IRX-β plane (Kong et al. 2004, Seibert et al. 2005, Takeuchi et al. 2010, Oteo et al. 2013, Casey et al. 2014. Various groups have ventured to quantify this scatter and find the physical drivers of it. For instance, the broad distribution of stellar population ages is presumed to affect the resulting IRX-β of a galaxy, increasing the scatter and leading to a systematic shift along the β-axis as the median stellar population age varies across galaxies (with β increasing with stellar population age, Kong et al. 2004, Boquien et al. 2009).
Measurements of dusty star forming galaxies (DSFGs) and ultra luminous infrared galaxies (ULIRGs), which are characterized by a central star forming region the size of a few hundred pc, and 90% of their energy output being in the infrared regime, have shown that they have significantly higher IRX than usual galaxies, when compared to normal local galaxies with similar β (Casey et al. 2014, Goldader et al. 2002. A proposed explanation for this is that these highly star forming galaxies have very complex geometries regarding their dust distribution, leading to parts of the UV emission being only weakly attenuated, with other parts being heavily obscured by dust (Calzetti 2001, Casey et al. 2014. This leads to a UV continuum dominated by young stars not obscured by dust, causing low values of β, while the IR continuum is being dominated by regions with heavy dust obscuration, causing high IRX (see also the theoretical explanations in Popping et al. 2017b andNarayanan et al. 2018).
Some galaxies lie significantly below and to the right of the original Meurer relation. It has been proposed that the grain typeand size-distribution of the ISM dust affect the exact location in the IRX-β plane, by influencing the far UV (FUV) extinction curve (Pettini et al. 1998, Witt & Gordon 2000. For example, dust as it is present in the small and large magellanic clouds (SMC and LMC) has a far higher FUV extinction than Milky Way dust, which leads to a stronger dependence of β on the dust opacity. Ultimately, this translates to a less steep IRX-β relation for galaxies with strongly FUV-extincting dust types.
In short, variations in dust properties, dust geometry, as well as stellar population age all affect the location of galaxies in the IRX-β plane, and this leads to an apparent scatter which is present at all redshifts. Other drivers for the observed scatter, such as dust temperature (Narayanan et al. 2018;Faisst et al. 2017), or variations in intrinsic β (the slope of the UV continuum without any dust attenuation, Boquien et al. 2012) have also been suggested, but a definite answer as to which are the main mechanisms behind the observed scatter is yet to be found. The thruth may lie in a combination of all proposed explanations and in the fact that normal star forming galaxies typically consist of many star forming regions disconnected from each other, each with their own stellar population ages, dust types and dust geometries.
Apart from observational investigations, analytical models (e.g., Popping et al. 2017b), employing two component dust models and stellar population synthesis have confirmed the previously mentioned ideas of causes for scatter in the IRX-β plane. However, even though analytical modeling gives us valuable insights into the effects several physical properties have on the IRX-β scatter in a controlled environment, it only provides a simplified view of the subject, and struggles to model the complex interplay of the many variables at hand.
Numerical simulations have succeeded to account for the complex dynamics of galaxies and the interplay between various physical quantities by coupling idealized (Safarzadeh et al. 2017) or zoom-in simulations (Narayanan et al. 2018, Ma et al. 2019) with radiative transfer calculations. These works confirmed a tight IRXβ relation for Milky Way type galaxies even at higher redshifts, as well as higher dust temperatures in high-redshift galaxies partly due to increased star formation rates. The simulations also suggest that certain galaxy types such as starbursts, DSFGs and ULIRGs lie above the Meurer relation, due to complex dust geometries that lead to spatially disconnected UV-and IR-emission. In these cases the attenuated UV-emission is often spatially very extended and dim.
When it comes to numerically simulating galaxies, there have been three main types of simulations: (i) idealized simulations, that simulate single isolated galaxies or galaxy mergers without any cosmological context, (ii) zoom-in simulations, that trace the evolution of small samples of galaxies in a cosmological context through time, with very high resolution (M baryons 10 4 M ), and initial conditions taken from large volume dark matter only simulations, or (iii) full physics large volume cosmological simulations, that start from primordial initial conditions and trace the evolution of a large sets of galaxies in a large cosmic volume through time, with intermediate resolutions (M baryons ≈ 10 6 M ). The reason that other groups have focused on idealized or zoom-in simulations to study the IRX-β relation is simple: complex dust geometries promise to be an explanation for the IRX-β scatter across galaxies, and these geometries had to be resolved well enough in simulations. So far such resolutions have only been achievable with these two simulation types. The downside of using zoom-in simulations is that in order to achieve these high resolutions, they are restricted to small galaxy samples, not necessarily representative of the overall galaxy population. On the other hand, typical large scale simulations have not been an ideal option for IRX-β investigations, as their numerical resolution was too low and did not properly capture the complex distribution of stars and the ISM.
This has changed now, with the new TNG50 simulation (Nelson et al. 2019a, Pillepich et al. 2019. It is the highest resolution variant of the IllustrisTNG simulation suite , Naiman et al. 2018, Nelson et al. 2018, Pillepich et al. 2018b, Springel et al. 2018, and has a mass resolution of M baryons ∼ 8.5 × 10 4 M in a cosmological volume of ∼ 50 cMpc on a side, bridging the gap between zoom-in and large scale cosmological simulations, making it ideal to study a representative sample of galaxies over cosmic time at sufficiently good mass resolution. It is especially well suited to this work, because it combines a high resolution required to accurately model the dust distribution of each galaxy with a sufficiently large sample size of galaxies which is needed to perform a sensible analysis of the problem at hand.
In this work, we will try to quantify the scatter in the IRX-β dust attenuation relation. We will investigate most of the star forming galaxies in the cosmological volume of TNG50 between the current epoch and z = 4. We employ the radiative transfer code SKIRT (Baes et al. 2011, Camps & Baes 2015 to model the attenuation of stellar light by dust and obtain full spectra with high detail in UV and IR. We will calculate IRX and β for all of these galaxies and explore how various physical and observable quantities of those galaxies, which are directly taken or derived from the TNG50 simulation data, affect the location of those galaxies in the IRX-β plane. We will also study the effect different dust properties have on the location of galaxies in the IRX-β plane.
This paper is organized as follows. In Section 2 we describe the numerical framework by introducing the IllustrisTNG simulation suite and by describing the galaxy selection, as well as the SKIRT radiative transfer code. In such section, we also explain how the IRX and β are derived. In Section 3, we show our main results, quantifying the scatter and locus of the IRX-β relation. Section 4 contains a discussion and possible interpretations of our findings. It is followed by an outlook on ideas for future work and finally a summary of our main results in Section 5. Appendix A contains various checks on the dependence of the fiducial radiative transfer model on a number of implementation choices.

The TNG50 Simulation
The IllustrisTNG Project (TNG hereafter: Nelson et al. 2018;Marinacci et al. 2018;Springel et al. 2018;Naiman et al. 2018;Pillepich et al. 2018b) is a suite of magneto-hydrodynamical cosmological simulations for the formation of galaxies employing the moving mesh AREPO (Springel 2010) code. In this work, we use the output of its highest-resolution implementation: the TNG50 simulation (Nelson et al. 2019a;Pillepich et al. 2019).
The TNG project is based on an updated version of the former Illustris galaxy formation model (Vogelsberger et al. 2014a, Vogelsberger et al. 2014b, Torrey et al. 2015 described in Weinberger et al. (2017) and Pillepich et al. (2018a). The updates it received in comparison to the original Illustris model include the addition of magneto-hydrodynamics, a new model of Active Galactic Nuclei (AGN) feedback for low accretion rates (Weinberger et al. 2017), and modifications to the galactic winds, stellar evolution and chemical enrichment schemes (Pillepich et al. 2018a). The simulations trace the evolution in time of dark matter particles, stellar particles and stellar winds, gas cells and black holes, starting from initial conditions at z = 127, ending at z = 0. Identification of haloes and subhaloes was done with the Friends-of-Friends (FoF, Davis et al. 1985) and SUBFIND (Springel et al. 2001b, Dolag et al. 2009) algorithms. The subhaloes within each FoF halo contain all resolution elements that are gravitationally bound to them. Within the TNG50 framework, we consider those subhaloes to be galaxies which have non-vanishing stellar mass and a total dark matter mass fraction of at least 20% (see Nelson et al. 2019b, Section 5.2;or Pillepich et al. 2019, Section 3.1). The IllustrisTNG model roughly matches the observed cosmic SFRD for 0 < z < 8, the galaxy mass function at z = 0, the stellar-to-halo mass relation at z = 0, the halo gas fraction at z = 0 and galaxy sizes at z = 0. IllustrisTNG uses cosmological parameters consistent with the values recently measured by the Planck collaboration (Planck Collaboration et al. 2016).
There are three different flagship realizations of the Illus-trisTNG suite, namely TNG300, TNG100 and TNG50 -named after their simulation box sizes at z = 0 of 50, 100 and 300 comoving Mpc, respectively. They provide different mass resolutions and were designed for different purposes: TNG300 simulates the largest cosmological volume of them all, and thus provides insights on the large scale physics with a large number of simulated galaxies as well as rare objects, like massive galaxy clusters. TNG50, with its smaller volume delivers relatively less statistics, but provides a numerical resolution that is comparable to those of zoom-in simulations, making it possible to study detailed structures of galaxies. TNG100 lies in between the two. All variations also have lowerresolution variants (so that resolution effects can be quantified), and dark matter-only variants. In this work, we use data of the highest resolution version TNG50 (Nelson et al. 2019a, Pillepich et al. 2019, see also Table A3).
TNG50 follows the evolution of 2×2160 3 resolution elements inside a cube measuring 35 cMpc/h on each side, which correpsonds to a volume of 51.7 3 (cMpc) 3 . This translates to a target mass resolution of 0.85 × 10 5 M for the baryonic resolution elements (gas cells and stellar particles) and to 4.5 × 10 5 M for the dark matter resolution elements. Subscale physical mechanisms (below 70-150 parsecs within the star forming regions of galaxies, see Figure 1 of Pillepich et al. 2019) are included via subgrid modeling, this being the case for example for star formation, supernovae and AGN feedback (Springel & Hernquist 2003, Pillepich et al. 2018a, Weinberger et al. 2017. In this simulation, galaxies with stellar mass M * = 10 9 M and above will contain at least about 10,000 stellar particles, so their resolution is in between those of typical cosmological simulations and those of zoom-in simulations. We will make use of a number of TNG50 data output, including various physical quantities of stellar particles and gas cells as well as integrated properties of the halos and subhaloes identified by the SUBFIND algorithm. We will later use these data to identify drivers of scatter in the IRX-β relation.

TNG50 Galaxy Selection
In this work, we focus on star-forming galaxies, meaning those that lie on the SFMS and above. The classification of SFMS galaxies is done as in Pillepich et al. (2019), Section 4.
In short, we stack all TNG50 galaxies into 0.2 dex stellar mass bins in the range M * = 10 8 M to 10 10.2 M . Then, we calculate the median specific star formation rate (sSFR) inside twice the 3D stellar half mass radius r 1/2 in each bin. The sSFR of a galaxy is defined as the instantaneous star formation rate (SFR) inside 2r 1/2 divided by the stellar mass inside 2r 1/2 . This sSFR median defines the SFMS ridge, and green-valley and quenched galaxies are "rejected" as those that have a logarithmic distance of −1 < ∆ log 10 sSFR < −0.5 or ∆ log 10 sSFR < −1, respectively. The process is then repeated, without considering green valley and quenched galaxies, until the median converges. A power law is then fitted on the sSFR median to extrapolate the SFMS to higher masses Figure 1. The selection of TNG50 galaxies used in this work. Here we show the TNG50 galaxy population at the edges of the considered redshift range, z = 0 (left) and z = 4 (right), according to the galaxies' star formation activity and stellar mass (where SFR and stellar mass are taken inside twice the 3D stellar half mass radius r 1/2 ). Gray circles denote all TNG50 galaxies, red filled circles represent those we select for this work: star-forming and well resolved i.e. with M * > 10 9 M ). Notice that as we increase redshift, fewer galaxies are found in the selected mass range (see also Table 1). M * > 10 10.2 M . The galaxies with ∆ log 10 sSFR > −0.5 are then the star-forming galaxies of the SFMS.
We make no distinction between central and satellite galaxies. We do not inspect low-SFR galaxies, because their dust content is expected to be comparably low, making them less prone to dust attenuation effects. We note that this transition is gradual, as for example green valley galaxies' stellar light can be notably attenuated by dust, however we stick to the well-defined sample of SFMS galaxies in our investigations. We select star-forming galaxies across different redshifts (z = 0, 0.5, 1, 2, 3, 4), so as to investigate a possible redshift-dependent evolution in the IRX-β relation. We limit ourselves to redshifts z 4, because at higher redshifts the effect of the cosmic microwave background (CMB) on the infrared emission of galaxies, which is not accounted for in the radiative transfer modeling employed in this work, becomes nonnegligible.
We investigate galaxies that are made of at least 10,000 stellar particles to ensure that they are well enough resolved to reproduce complex dust geometries. This translates to a minimum stellar mass of M * 10 9 M for TNG50. Notice that as we increase redshift, the median stellar mass of galaxies is decreasing, so the higher we go in redshift, the smaller the sample size will be given our selection cut of M * > 10 9 M (see Table 1). In total, we end up with 7280 well resolved SFMS galaxies in the redshift range 0 z 4.

Obtaining TNG50 Galaxy Spectra with SKIRT
The IllustrisTNG simulation suite does not track radiative processes. Hence, in order to accurately model attenuation of TNG50 galaxy spectra by dust in the ISM, we employ radiative transfer calculations in post-processing. We make use of the publicly available radiative transfer code SKIRT 1 (Baes et al. 2011, Camps & Baes 2015, which performs a Monte-Carlo tracing of the paths of  Table 1. A summary of the number of TNG50 galaxies used in this work. From left to right: redshift, total number of galaxies with mass M * > 10 9 M , number of selected SFMS galaxies given our selection criteria at the given redshift. At high redshifts, we have in total fewer galaxies, but a higher percentage of them belongs to the SFMS. photons emitted by a stellar distribution as they are scattered and absorbed by a dust distribution. These photons are captured by a virtual observing instrument which is capable of returning the total spectral energy distribution (SED) of the captured photons, as well as a spatially resolved spectrum (a data cube containing a 2D image of the observed fluxes at each specified wavelength). The stellar distribution used in our SKIRT calculations is predicted by the TNG model and thus directly imported from the TNG50 simulation output. We note that ISM dust is not tracked in the IllustrisTNG simulations. We will therefore assume that the simulations' metal distribution is a tracer for the ISM dust distribution, following the approach of similar studies (e.g. Camps et al. 2016, Trayford et al. 2017, Ma et al. 2019, Rodriguez-Gomez et al. 2019, Vogelsberger et al. 2019a, see Section 2.3.2 for more details).

Stellar Sources
Each stellar particle in TNG50 is a point-like coeval stellar population. For the treatment with SKIRT, we define an adaptive smoothing scale to each stellar particle equal to the 3D distance to the N-th nearest neighbouring stellar particle, which has been described in Torrey et al. (2015). This is to allow SKIRT to calculate internally a smoothed photon source distribution function. Torrey et al. (2015) use N = 16, however we orient ourselves on Rodriguez-Gomez et al. (2019), who have chosen N = 32 and have also found that varying N in the range of 4 to 64 has no particular effect at least on morphological measurements. We tested if this is the case for IRXβ measurements, too, and found no noticeable difference between choosing N = 16, 32 or 64. All smoothing scale calculations are carried out by the default smoothed particle hydrodynamics (SPH) spline kernel (Monaghan & Lattanzio 1985), redefined over the interval [0, h sml ] (Springel et al. 2001a), where h sml is the smoothing scale described before: where q ≡ r/h sml , r is the radial distance, n is the number of dimensions, and σ is a normalization constant with the values 2/3, 10/7π, 1/π, in one, two and three dimensions, respectively. The SEDs of stellar particles older than 10 Myr are modeled with the GALAXEV population synthesis code (Bruzual & Charlot 2003), which is implemented internally in SKIRT. This code uses simple stellar population (SSP) models, which were computed using Padova 1994 evolutionary tracks and a Chabrier (2003) initial mass function (IMF). These models provide the rest-frame luminosity per unit wavelength of an SSP, L λ (λ, t, Z) as a function of wavelength λ, age t, and metallicity Z (the luminosity is normalized to solar mass M ). The luminosity is provided by an interpolation of a grid which is sampled at 221 unevenly spaced ages between 0 and 20 Gyr, 7 metallicities between 10 −4 and 0.5 and 1221 wavelengths between 91 Å and 160 µm. To generate a stellar spectrum, the GALAXEV code requires as an input the initial mass of the stellar population (neglecting mass loss due to stellar evolution), its metallicity, and stellar age.
Young stellar populations present in starbursting regions are formed and embedded in dense and cold molecular dust, which we call "birth clouds". As the lifetime of these molecular birth clouds is about 10 Myr (Murray & Rahman 2010), we assume that stellar populations younger than that are still surrounded by these clouds. These clouds are not resolved by TNG50, as their sizes are on pc scales. To account for attenuation in birth clouds, all stars younger than 10 Myr will be treated as starbursting regions in this work, and their SEDs are modeled using the MAPPINGS-III photoionization code (Groves et al. 2008), which includes emission from HII-regions and their surrounding photodissociation-regions (PDRs), as well as absorption by gas and dust in the birth clouds. The MAPPINGS-III models require the following five input parameters for each star forming region: (i) The SFR, assumed to be constant over the last 10 Myr, (ii) the metallicity, (iii) the compactness parameter C, (iv) the ISM pressure P ISM and finally (v) the PDR covering fraction f PDR . For every young stellar particle, we assume that the SFR is given by its initial mass divided by 10 Myr (to ensure mass conservation), and use its nominal metallicity inherited from its parent gas cell. For the compactness C and the ISM pressure P ISM , we use typical values found in literature (Groves et al. 2008) of log 10 C = 5 and log 10 The exact values of these two parameters only have a noticeable influence on the far-IR SEDs, we tested what effect a variation of these values has on the result of IRX-β (see Appendix A1.3). While it makes a difference whether we include MAPPINGS-III or not, variations of C and P ISM have only marginal effects on the resulting total IRX-β of a galaxy. Following Jonsson et al. (2009) and Rodriguez-Gomez et al. (2019), we adopt a value of 0.2 for the covering fraction f PDR for our fiducial model. Variation of this value also has no significant qualitative impact on the results (see Appendix A1.3).

Dust Modeling
SKIRT offers the option of performing the radiative transfer calculations directly on a three-dimensional Voronoi mesh (Camps et al. 2013), which makes it particularly well suited to the IllustrisTNG simulation suite, which is based on the AREPO code (Springel 2010). We can directly reconstruct the gas distributions exactly as they are implemented in the hydrodynamic solver (except for cell gradients) in order to evolve the system. In practice, we define a cubical region around each inspected galaxy, with a side length of 15r 1/2 (with r 1/2 being the half stellar mass radius of the galaxy), to make sure that all of the parts of the dust distribution relevant to radiative transfer are captured. This ensures that even the contribution of very spatially extended dust is included in our simulations. The coordinates of the gas cells, which are in actuality mesh generating points, are used by SKIRT to reconstruct the Voronoi mesh inside this volume using the VORO++ open source library for computing Voronoi tesselations. Together with the density values for each gas cell, this fully describes a gas density distribution.
We assume that the diffuse dust content of a galaxy is traced by the gas-phase metal distribution assigned to this galaxy. This means that we only include the gas cells which are gravitationally bound to the subhalo the galaxy is located in, all other gas cells' densities are set to zero beforehand. We add the condition that dust can only be present in a gas cell if the gas cell is either star forming (SFR > 0), or if its temperature is less than a typical threshold value (T < T max , T max = 75000 K). This means that we effectively exclude all gas cells that are too hot and have zero SFR as determined by TNG50. Cells that are star forming or that are below the temperature threshold are always considered to contain dust. In the TNG50 model, gas cells with a density above a threshold of n H ∝ 0.1 cm 3 are considered star-forming and are stochastically converted into stellar particles (Springel & Hernquist 2003, Pillepich et al. 2018a. We motivate these aforementioned choices by the fact that dust is rapidly destroyed in hot gas through thermal sputtering (Guhathakurta & Draine 1989). Unfortunately, since TNG50 does not model dust directly, we cannot properly constrain T max using a physically motivated procedure. Other groups working with the EAGLE simulations have chosen a value of T max = 8000 K (Camps et al. 2016, Trayford et al. 2017). We choose a much higher value so as to only eliminate the very hottest gas cells from potentially containing dust. We accomplish this by manually setting the metallicity of all gas cells that will contain no dust to zero, which in turn sets the dust density of these cells to zero.
Only a certain fraction of the metal content of a gas cell will be locked up in the form of dust. This fraction, the dust-to-metal ratio, is treated as a free parameter in this work. Orienting ourselves on other studies, we assume a constant fiducial value of this dust to metal ratio f = 0.3 (Camps et al. 2016). Observational work (Rémy-Ruyer et al. 2014, De Vis et al. 2019) and theoretical models (e.g., McKinnon et al. 2017, Popping et al. 2017a suggest that this fraction evolves as a function of gas-phase metallicity, but the exact shape of this relation is still uncertain. It was also suggested that the dust-to-metal ratio is evolving with redshift, which is why . SKIRT spectral flux density spectra of exemplary high mass galaxies of TNG50 at z = 0.5 (see their images in Figures 3 and 4). For these spectra, we reran SKIRT with a higher wavelength detail of 300 logarithmically distributed wavelength grid points and a higher photon package number of 10 7 instead of 10 6 . The orange lines show the unattenuated spectrum, the blue lines show the attenuated spectrum, assuming Milky Way type dust properties. The grey shaded area marks the region in which we perform the fit to the attenuated (unattenuated) UV spectrum to attain its slope β (β 0 ) (see section 2.4), which is indicated by the solid black transparent lines. The red shaded region shows the range over which we integrate to get F IR . The green dashed line shows the flux density spectrum F λ = λ f λ vs λ, from which we get F UV at λ = 0.16 µm (indicated by the black dashed line). Top row: Galaxies with relatively high IRX and high β. They are strongly attenuated, their UV slope changes significantly, and their UV spectral flux density is greatly decreased. Bottom row: Galaxies with relatively low IRX and low β. They are hardly attenuated at all, and show no significant change in UV slope, nor in UV spectral flux density. by comparing with observational results. For simplicity, we will assume a fixed value in our fiducial model. We will vary this parameter later in this study to see the effect it has on the IRX-β distribution of the galaxies (see section A1.1).
The dust density distribution is therefore derived from the TNG50 gas density distribution the following way: where ρ dust is the dust density, f is the dust-to-metal ratio, ρ gas is the gas density, Z is the gas metallicity, T gas is the gas temperature and SFR gas is the instantaneous SFR of the gas resolution elements. Note that the gas metallicity Z is a factor in this equation -higher metallicity systems will hence have a higher dust-to-gas ratio via this definition. SKIRT offers various options to model the dust grain composition. As mentioned in Section 1, it is expected that the choice here will affect the resulting IRX-β of our simulated galaxies. For our fiducial model, following Camps et al. (2016) and Trayford et al. (2017), we choose the Zubko et al. (2004) multi-component dust mix, which models a composition of graphite, silicate and poly-cyclic aromatic hydrocarbon (PAH) grains, with various grain size bins for each grain type. The size distributions and the relative amount of the dust grains are chosen so as to recreate the properties of Milky Way type dust. We also later vary the dust type to recreate the dust of the Large Magellanic Cloud (LMC) and Small Magellanic Cloud (SMC), by employing a Weingartner & Draine (2001) grain size distribution for the LMC and SMC dust respectively.
The dust of the molecular birth clouds mentioned in section 2.3.1 is treated separately from the diffuse ISM dust, and is already included by the usage of the MAPPINGS-III models for starforming regions before any radiative transfer calculations are being made. We assume that stellar particles older than 10 Myr will not be surrounded by birth clouds anymore, so their emission is only being absorbed and scattered by the ambient ISM, which is modelled using SKIRT. We also repeat the SKIRT run for each galaxy without the presence of any ISM dust and MAPPINGS-III birth clouds for the calculation of the intrinsic (unattenuated) UV-slopes β 0 .

Virtual Instrument Setup
Each galaxy is observed from an angle perpendicular to the xyplane of the simulation box, from a distance of 10 Mpc. This leads to a random orientation of the galaxies within the sample with respect to our viewing angle. We choose a field of view equal to the box size of the dust distribution which is 15r 1/2 , to again make sure that all of the parts relevant to radiative transfer in the galaxy are captured in the output spectra. This ensures that even very spatially extended UV emission, which might cause some IRX-β deviations, is captured by our instrument. For the spatially resolved spectra, we choose a resolution of 100 by 100 pixels. These images will be used to investigate the spatial distribution of UV-and IR-emission. This resolution is high enough to differentiate the galaxies' structures by eye.

Performing the SKIRT Runs
SKIRT is a Monte-Carlo radiative transfer code, meaning that the radiation field is represented as a discrete, large number of photon packages ("rays"). Typical numbers of rays per wavelength go from 10 6 to 10 8 . In our fiducial model, we choose a value of 10 6 rays per wavelength, which is a number that is capable of producing converged spectra. Our fiducial wavelength grid is defined by a custom wavelength grid file, which contains a total of 59 wavelength grid points, with 30 evenly spaced grid points in the relevant UV-range 0.123 µm < λ < 0.32 µm, and 20 evenly spaced grid points in the IR-range 8 µm < λ < 1000 µm, and 9 grid points that are evenly spaced outside of these ranges, so that the whole grid spans the range 0.1 µm < λ < 1001 µm. This resolution is sufficient both for calculating the UV-slope β and integrating the total IR-emission, even though very fine emission lines from the MAPPINGS-III models are not resolved. We performed a test run with a narrower spaced wavelength grid and found no significant deviations from the results produced by the fiducial model.

SKIRT Output
The output of the SKIRT runs are (i) full SEDs of each galaxy (given as flux density F λ = λ f λ in W/m 2 vs wavelength λ in µm) (see Figure 2), and (ii) data cubes containing resolved 100 × 100 pixel intensity maps for each wavelength of the specified wavelength grid, in units of W/m 2 /arcsec 2 (see Figures 3 and 4) which can then be combined, e.g. in order to create RGB images (see Figure 5). One can clearly see the effect that dust has on the appearance of galaxies at different wavelengths.
Summarizing, SKIRT gives us the freedom to setup the physical models and the numerical specifications of the Monte-Carlo simulation. For the physical models, we have to specify a number of input parameters, which are either taken from the TNG50 output, from the literature, or they can be calibrated, or taken as a free parameter. For a summary of our SKIRT setup, and a comparison to setups of other works coupling the IllustrisTNG model to SKIRT, see Table A1. For an overview of the adopted SKIRT input parameters and their origin please see Table A2.

Measuring IRX and Beta
The infrared excess IRX is defined as the ratio between the IR l'uminosity and the UV luminosity and hence quantifies the amount of dust obscured emission from the galaxy: where the IR luminosity L IR (the flux density F IR ) is attained from integrating the spectral luminosity L λ (the spectral flux density f λ ) over the wavelength range 8 µm < λ < 1000 µm, and the UV luminosity L UV (flux density F UV ) is the luminosity (flux density) at 1600 Å. In this work we take the UV luminosity (flux density) at λ = 0.1601 µm, the value in our wavelength grid that lies closest to that. Larger amounts of dust imply larger values of IRX. β is defined as the rest frame UV spectral slope of a galaxy (where we assume that the spectrum follows a power law): where f λ is the spectral flux density in the in the UV range. β is fitted to the spectrum in the wavelength range 0.123 µm < λ < 0.32 µm, which corresponds to UV radiation from O and B stars and ensures a broad wavelength coverage around the 2175 Å dust feature (if present, as for example in Milky Way dust). Typically the spectrum of a galaxy is such that β is negative: younger stellar populations produce a spectrum with a more negative (or bluer) β than older stellar populations. Dust affects a galaxy's spectrum by reddening it and thus by making the UV slope less negative.

RESULTS
In the following section, we investigate the IRX-β relation of galaxies as predicted by the TNG50 simulation. We compare our sample to a set of reference relations: (i) the original Meurer et al. (1999) relation for local starbursts, (ii) recent calibrations by Overzier et al. (2010) for local galaxies, (iii) the Casey et al. (2014) relation which is a fit to low-redshift galaxies, and includes a higher dynamic SFR range and also aperture-corrected data of heterogeneous samples, and Pettini et al. (1998) data for SMC-like dust attenuation curves.
3.1 IRX-Beta of TNG50 Galaxies Figure 6 shows a scatter plot of the IRX-β plane for all galaxies in our sample, color coded by redshift. The bulk of the sample lies very close to the reference relations, and seems to be best approximated by the Overzier et al. (2010) relation (we demonstrate this in practice in the Discussion section). A typical galaxy in our sample has a β between −2 and 0, and a log 10 IRX between 0 and 1.5. We find that higher redshift galaxies lie towards lower β at the same IRX: at z = 4, the median β is −1.80, while the median β at z = 0 is −1.36 (see also Figure 7, top panel). The median log 10 IRX is close to 0.5 for all redshifts. High-redshift galaxies seem to follow a steeper IRX-β relation than low-redshift galaxies, with less scatter in β: the standard deviation of β at z = 0 is σ β (z = 0) = 0.43, while at z = 4 it is only σ β (z = 4) = 0.31. When compared to the scatter that has been found in observations (e.g. Casey et al. 2014), the total scatter in our sample is smaller: at its worst, there is only a scatter of about 0.5 dex in IRX and a scatter of 1 in β (compared to a scatter of ca. 1 dex in IRX and 1 in β in Casey et al. 2014). These deviations can have multiple causes: e.g. inclusion or not of observational measurement uncertainties, differences in the sample selections, variations in stellar population age and variations in dust properties and geometry: see section 4.1.1 for an elaboration on this. In the following, we will present results that  . We show selected slices of the SKIRT output data cubes for two high mass galaxies of TNG50 at z = 0.5 with high IRX and high β (≈ 10 10−11 M ): see their spectra in Figure 2, top row). For these images we reran SKIRT with 300 logarithmically distributed wavelength grid points, a higher photon package number of 10 7 instead of the fiducial 10 6 and a 500 × 500 pixel image resolution. For each galaxy, the respective top row shows the result when dust is included in the radiative transfer calculations, the respective bottom row shows the result when dust is excluded. From left to right, we show the spatial distribution of the UV flux (0.16 µm), the V-band flux (0.43 µm), and the IR flux (102 µm) as detected by a virtual instrument. The white bar in each image denotes a spatial scale of twice the stellar half mass radius r 1/2 of the depicted galaxy. For these strongly attenuated galaxies, the spatial distribution of their flux does not seem to change significantly when including dust attenuation. Noticeable is the strong change in UV and IR flux.  . We show selected slices of the SKIRT output data cubes for two high mass galaxies of TNG50 at z = 0.5 with low IRX and low β (see their spectra in Figure 2, bottom row). Annotations are as in Figure 3. For these weakly attenuated galaxies, too, the spatial distribution of their flux does not change significantly when including dust attenuation. There is less decrease in UV flux when compared to the high IRX galaxies.
help pinpoint those possible causes for the scatter and the redshift dependent trends in IRX-β predicted by our model. Figure 7, top panel, shows explicitly that, according to TNG50 and assuming a Milky Way like dust and a fixed dust-to-metal ratio throughout, the IRX-β steepens towards higher redshifts.

Accounting for the Intrinsic UV-Slope
Differences in stellar population age can be a plausible source for scatter in the IRX-β relation and for the deviations from a local TNG50 z = 0.5: ID = 258162 2r 1/2 8.6pkpc log 10 M * = 11.17 log 10 M without dust without dust TNG50 z = 0.5: ID = 211991 2r 1/2 11.8pkpc log 10 M * = 11.35 log 10 M without dust without dust TNG50 z = 0.5: ID = 251705 2r 1/2 13.8pkpc log 10 M * = 11.11 log 10 M without dust without dust TNG50 z = 0.5: ID = 160812 2r 1/2 17.2pkpc log 10 M * = 11.26 log 10 M without dust without dust with dust TNG50 z = 0.5: ID = 258162 2r 1/2 8.6pkpc log 10 M * = 11.17 log 10 M with dust with dust TNG50 z = 0.5: ID = 211991 2r 1/2 11.8pkpc log 10 M * = 11.35 log 10 M with dust with dust TNG50 z = 0.5: ID = 251705 2r 1/2 13.8pkpc log 10 M * = 11.11 log 10 M with dust with dust TNG50 z = 0.5: ID = 160812 2r 1/2 17.2pkpc log 10 M * = 11.26 log 10 M with dust Figure 5. RGB images of the four high mass TNG50 galaxies of Figures 3 and 4 at z = 0.5 produced from wavelength slices at 0.675µm (red), 0.528µm (green) and 0.452µm (blue). The effect of dust on the perceived flux is clearly visible -many regions with high concentrations of young stars are covered by dust clouds when accounting for dust by application of radiative transfer models, affecting the resulting UV spectrum (see also Figure 2). Note also the reddening of some individual stellar particles due to the inclusion of MAPPINGS-III spectra for very young stars, caused by effectively creating birth clouds in these star forming regions.
shape. As the stellar population age increases, the intrinsic UVslope (β 0 ) of the stellar population becomes redder (i.e. a less negative UV slope). This naturally creates scatter in β, without even invoking dust absorption.
For all the galaxies in our sample, their intrinsic UV-slopes β 0 can be obtained by running completely dust-free radiative transfer simulations (also excluding birth clouds) on them, and then measuring the slopes of the resulting unattenuated UV-spectra. In our model, these intrinsic UV-slopes depend solely on the galaxies' stellar metallicities and stellar population ages.
In the middle panel of Figure 7, we plot the β 0 of all galaxies in our sample against their mass-weighted average stellar population ages t age . Higher redshift galaxies tend to have younger stellar populations, and this correlates with median intrinsic UVspectra more "blue" (a more negative UV slope) when compared to low redshift UV-spectra. The median β 0 at z = 4 is, for example β 0, median (z = 4) = −2.30, while at z = 0, β 0, median (z = 0) = −2.07. At low redshifts, the spread in β 0 is also higher than at higher redshifts. At z = 4, the standard deviation of β 0 is σ β 0 (z = 4) = 0.10 and at z = 0, it is σ β 0 (z = 0) = 0.14.
In the bottom panel of Figure 7 we now account for differences in stellar population age and stellar metallicities by plotting IRX against β − β 0 . The colored contours show the distributions of the galaxies in this IRX-(β − β 0 ) plane at each redshift. The value of β − β 0 is now a more direct measure of the amount of UVattenuation in a galaxy, because it takes into account the natural fluctuations in the intrinsic UV-slopes β 0 . The lower the value of β − β 0 , the less "reddening" due to dust a galaxy has experienced. The redshift trend is now significantly reduced. Still, high-redshift galaxies have on average lower β − β 0 than low-redshift galaxies.
The IRX-(β − β 0 ) distribution is tighter than the IRX-β dis-tribution. The remaining scatter in IRX of about ∆ log 10 IRX ≈ 0.5 in the IRX-(β − β 0 ) relation can not be caused by dust composition differences (all galaxies are assigned the same Milky Way dust type in our fiducial model), nor due to differences in their intrinsic UVslope (we corrected for β 0 ) and must therefore be caused by variations in other galaxy properties or dust geometry. In the following, we will investigate how different galaxy properties correlate with the location of galaxies in the IRX-β plane.

Examining Correlations in the IRX-Beta Plane
We investigate several physical quantities to see if they correlate with the galaxies' IRX-β distributions. These include the star formation rates SFR, the average gas metallicity Z Gas , the specific star formation rate sSFR, the stellar mass M * , the star formation efficiency SF eff , and the gas mass M Gas . All these quantities are measured within two times the stellar half mass radius of the modeled galaxies.
The panels in Figure 8 show the IRX-β distribution of TNG50 SFMS galaxies at z = 1 with M * > 10 9 M . These plots are divided into 100 by 100 grid cells, the color of each grid cell corresponds to the median value of the physical property of interest of all the galaxies that lie inside of it in the IRX-β plane. By plotting the IRX-β distribution at fixed redshift, possible redshift trends that might be present in the whole sample across all redshifts are eliminated. We show z = 1 as a representative redshift and we have checked that the findings below hold at all individual cosmic epochs between z = 0 and z = 4.
• Star Formation Rate (SFR). Except for a handful of sources with β > −1 that fall below the locally derived reference relations, . The IRX-β distribution of the complete sample of selected TNG50 star-forming galaxies at all investigated redshifts. The lines show the various reference relations we compare our result to; colored circles represent galaxies from our selection in the z = 0−4 redshift range, with N gal denoting the selected and depicted number of objects. The bulk of the galaxies agrees with the reference relations, the best agreement being with the Overzier curve (Overzier et al. 2010). We uncover a redshift trend: in general, at a given IRX, high-redshift galaxies have lower β than low-redshift galaxies (see also Figure 7, top left panel). This indicates that high-redshift galaxies follow a different IRX-β relation from low-redshift galaxies, under the assumption that the type of dust and dust-to-metal ratio do not change over cosmic time. TNG50 also predicts a non-negligible intrinsic scatter of the IRX-β relation, for example at redshift 4, where most galaxies have a β close to -2, their log 10 IRX ranges from 0 to 1. the IRX-β relation depends on SFR (Figure 8, top left panel): at fixed IRX, higher SFR implies lower β values. In fact, the layering due to different SFRs is somewhat parallel to the median relation: in other words, at fixed β, higher SFR also implies larger IRX values. The overall effect would be accentuated if galaxy populations across cosmic epochs were considered, as the galaxy populations exhibit lower SFRs towards lower redshift at fixed mass.
• Specific Star Formation Rate (sSFR). The left panel in the second row of Figure 8 shows how the sSFR correlates with the IRXβ distribution. When looking at z = 1, the correlation with sSFR is somewhat weaker, albeit still present, in comparison to e.g. the effects of SFR, and more prominent for galaxies with low IRX values: higher sSFR galaxies tend to lie towards lower β in the plane. This makes sense intuitively, as high sSFR galaxies contain on average larger amounts of younger stellar populations. Also in this case, the trends would be stronger if galaxies across many cosmic epochs were considered at once, as higher redshift galaxies have systematically higher sSFR and younger stellar populations.
• Stellar Mass We discern an almost one-to-one correlation between mass and IRX when showing the IRX-β plane color coded by stellar mass (Figure 8, second row, right panel), with IRX increasing as a function of stellar mass. The trend at z = 1 with stel- . Top: redshift trend in IRX-β from TNG50 galaxies. The distribution of our complete galaxy sample is shown in grey datapoints, the colored contours highlight the galaxy distributions and locations at each redshift (contours include 95% of the population). At z < 2, galaxies generally follow the reference relations, while at z > 2, given IRX, galaxies are typically found at lower β than what is predicted by the local reference relations. Caption continues in next column. Figure 7. Middle: galaxies' intrinsic UV-slopes β 0 against their mass weighted mean stellar population ages t age . The intrinsic slope β 0 correlates with the stellar population age: the younger the stellar population, the lower is β 0 . High-redshift galaxies tend to have younger stellar populations than low-redshift galaxies and hence they have lower β 0 . Bottom: Effects of accounting for differences in stellar population ages and stellar metallicities across redshifts: we plot IRX against β − β 0 , which is now a direct measure for the amount of UV-attenuation. This serves as a first order correction to redshift trends in IRX-β. The overall IRX-(β − β 0 ) distribution is tighter, and the redshift trend has been reduced. There is still a slight shift in locus depending on redshift -this is probably because highredshift galaxies have lower gas metallicities, which lowers the amount of dust they contain following Equation 2, and therefore lowers the amount of attenuation β − β 0 . lar mass largely resembles the trend with SFR, which is a natural cause of looking at main-sequence galaxies at fixed redshift.
• Gas Mass and Star Formation Efficiency The star formation efficiency is defined as the star formation rate inside 2r 1/2 , divided by the galaxy's gas mass inside 2r 1/2 (this includes a contribution from atomic, molecular, and ionized gas): The higher this value, the more efficiently the gas of a galaxy is converted into stars, conversely a lower star formation efficiency corresponds to a galaxy that converts only little of its gas into stars. The correlation of the IRX-β distribution with star formation efficiency is very strong (see Figure 8, third row, left). The color gradient appears to be perpendicular to the reference relations. In general, a higher star formation efficiency leads to a shift to higher IRX and lower β with respect to the reference relations. We can attempt to dissect if this trend is driven by changes in gas mass or SFR by looking at these quantities independently. We already saw there is a clear gradient in the IRX-β plane when color coding galaxies as a function of their SFR. There is a correlation present with gas mass, but it is less strong (see Figure 8, third row, right panel). The IRX of galaxies on average increases with their gas mass along the IRX-β reference relation (i.e., they also have higher β). The gradient when color coding galaxies as a function of their SF eff resembles the gradient with SFR more than the gradient with gas mass, suggesting the SFR is the main parameter driving the gradient in SF eff .
• Gas Metallicity The average gas metallicity of a galaxy is an indicator of how much dust each of its gas cells contain via Equation 2. We find that IRX increases as a function of gas-phase metallicity, whereas β correlates with gas-phase metallicity less strongly (Figure 8, top right panel). As can be expected, an increased dust abundance results in a higher fraction of the emission being absorbed from the UV and re-emitted in the IR.
Note that the dependencies on star formation rate, gas metallicity, stellar mass, and star formation efficiency found in Figure 8 remain similarly strong also in the IRX-(β − β 0 ) plane (i.e. removing the effects of young stellar populations) while that with specific star formation rate is much weakened.

DISCUSSION
In this work, we obtained the IRX-β distribution of simulated galaxies of the TNG50 simulation, the highest resolution flagship run of the IllustrisTNG project. We performed radiative transfer  Figure 8. Correlations of the IRX-β distribution of the SFMS TNG50 galaxies at z = 1 and M * > 10 9 M with various physical galaxy properties, color coded. The IRX-β plane is divided into 100 by 100 pixels, each containing a different number of galaxies (at least 1). In each pixel, we measure the median value of the physical quantity we investigate, and color code the pixel by this value. We study, from top to bottom, left to right: star formation rate; average gas metallicity; specific star formation rate; stellar mass; star formation efficiency; and gas mass, all within two times the stellar half mass radius. calculations using SKIRT to realistically model the effects of dust attenuation in these galaxies. We have focused on main-sequence galaxies at z = 0, 0.5, 1, 2, 3 and 4, with stellar masses M * > 10 9 M .

The Intrinsic Scatter in IRX-Beta
The TNG50 galaxy population investigated in this study broadly agrees with the IRX-β relations that have been observationally established for local galaxies. We see similar trends with SFR as Casey et al. (2014), even though the latter report larger scatter than we find in TNG50, at low redshift. It should be noted from the onset that in this work we do not account for observational uncertainties that may be responsible for inflating the scatter in the observational data (e.g., when inferring the total infrared luminosities of galaxies). Nevertheless, we can speculate that a different amount of scatter can be due to sample selections: these may affect the observed data sets, but also the simulated one. For example, it could be that, because of the still limited cosmic volume covered by TNG50, we do not sample in TNG50 the most extremely star-forming galaxies (starburst, ultra luminous infrared galaxies, and dusty star-forming galaxies), which are rare. The lack of sources with high IRX in our simulations may also be the direct result of our choice for the dustto-metal fraction. An increase in the dust-to-metal fraction from our fiducial value 0.3 will result in a higher β and more so IRX (we find an increase by 0.5 dex when adopting the extreme scenario for local galaxies of a dust-to-metal fraction of 0.9, see Appendix A1.1), naturally increasing the range in the predicted IRX of galaxies.
At z > 1 we find little scatter in the IRX-β relation of simulated galaxies, similar to the scatter found in the observed IRXβ relation for galaxies on the main-sequence. The predictions by our model at z 2 are in agreement with the results by Bourne et al. (2017), Fudamoto et al. (2017), McLure et al. (2018), and Fudamoto et al. (2019 who find an IRX-β relation shifted bluewards of the z = 0 reference relations. Other works, on the other hand, (e.g., Bouwens et al. 2016;Reddy et al. 2012;Álvarez-Márquez et al. 2016) find that the IRX-β of normal star-forming galaxies is similar to the z = 0 reference relations, or even further towards the red.
The general agreement between our model predictions and observations is very encouraging. A more quantitative comparison will require additional steps to mimic the observational methodology, for instance by selecting galaxies in the exact same way as the aforementioned works. Additional scatter and uncertainty in the observations is introduced by assumptions about the dust temperature and other measurement errors (see e.g., the discussion on measuring β in Popping et al. 2017b and the discussion on dust temperature in Narayanan et al. 2018). In other cases the results are based on stacking (e.g., Bouwens et al. 2016) rather than individual detections. These assumptions and approaches are not accounted for in our analysis. Furthermore, it should be kept in mind that we assigned Milky Way type dust to all simulated galaxies globally. In reality there might be a continuum of dust grain types varying throughout galaxies and over time. This would lead to a more spread out observed IRX-β distribution with more scatter in observations. We will get back to this in Section 4.3.

Systematic Redshift Trends
We discovered a systematic shift toward lower β with increasing redshift. We found that the main driver behind this trend is the redshift evolution of the galaxies' intrinsic UV-slopes β 0 , which evolves with redshift due to the following reason: galaxies at higher redshifts tend to have younger stellar populations. This causes their β 0 to be more negative, or more "blue" (see Figure 7) than those of low-redshift galaxies.
After correcting for β 0 , there are still residual non-negligible trends in redshift. Since we assume a uniform dust type, these residual trends could be driven either by changes in the dust-abundance or changes in the dust geometry across galaxies, or both.
We noticed that high-redshift galaxies in our sample also have a smaller scatter in β than the low redshift ones. This can be explained firstly by the more limited sample size at high redshifts, which encompasses progressively smaller ranges in stellar mass. Another explanation can be given by connecting the scatter in sSFR of our selected galaxy population to the scatter in β 0 : as seen in Figure 1, the higher we go in redshift, the less scatter there is in sSFR and, because of its connection with star formation history, the less scatter there is in β 0 and thus in β. At lower redshifts, there is more variety in sSFR and by the same arguments also in β; additionally, there is more diversity in average gas metallicities of the modeled galaxies and hence different effective galaxy-wide dustto-gas fraction, causing the IRX-β distribution to be even more spread out horizontally.
We conclude that it is to be expected that the IRX-β relation of Milky Way type SFMS galaxies with masses above 10 9 M is different at each redshift, with a trend toward lower β as we go to higher redshifts. This is driven by systematic differences in stellar population ages and star formation history, stellar metallicities, and possibly dust geometry. When inferring SFR from such galaxies with the relations derived at z ≈ 0, neglecting this systematic shift would lead to an underestimation of their IRX, which would lead to underestimations of star formation rates (see Section 4.2).

Correlations between Galaxy Properties and IRX-Beta
Since in the fiducial model the dust type has been fixed to Milky Way type dust for all galaxies, and we account for age related systematic shifts in β by correcting for β 0 , the remaining scatter most certainly reflects the galaxy-to-galaxy property variations within the same population at given cosmic epoch, that in turn may determine differences in the geometry of the dust distribution around the light sources and differences in the dust abundance. Indeed, the correlations found in this work are in place both for the entire sample of galaxies across redshifts, as well as at individual redshift snapshots.
For the single redshift samples, the trend in SFR seems to be dominated by the trend in stellar mass (Figure 8, 2nd row, right)higher mass galaxies on the main sequence will have higher SFRs, and higher IRX. The sSFR accounts for mass trends in single redshift samples, because we divide by the stellar mass, and thus eliminate the inherent mass dependency of the SFR at fixed time. There, we do not see a correlation with IRX any more, but with β, which we suspect to be caused by different intrinsic UV-slopes β 0 . The higher the sSFR, the lower is β 0 and hence β. There is a strong correlation between a galaxy's star formation efficiency and its position in the IRX-β plane, perpendicular to the empirically derived reference relations. We saw that this is mostly driven by a correlation with SFR, rather than gas mass.  Table 2. Results of the fit of the TNG50 data to a modified Overzier-like relation, where we introduce a new parameter β z : see equation (6). The best fit β z are listed in the third column, for TNG50 galaxies at a given redshift (first column) whose number is indicated in the second column.
The number of galaxies that strongly deviate from this modified relation (N ∆IRXz (best fit)>0.5 , last column) is small, below 3 percent at any redshift. This implies that after accounting for a redshift shift, the locally observed relations hold for SFMS galaxies.
The correlation between the average gas metallicity and the IRX-β distribution can be explained as follows. Low metallicity galaxies contain less dust, and are therefore less opaque. This results in lower β (less UV-emission is absorbed and scattered away) and lower IRX. We find that the IRX of galaxies increases as a function of stellar mass. At fixed redshift, β also increases as a function of stellar mass. Higher mass SFMS galaxies generally have higher gas mass and dust-to-gas fractions (because they have higher gas metallicities), causing more dust obscuration and therefore higher IRX.
Putting all of this together, we find that the location of different tracks in the IRX-β plane (above and below the different reference relations) is best determined by the specific star formation rate and the star-formation efficiency of galaxies. The infrared excess IRX on the other hand scales more strongly with stellar mass and gasphase metallicity.

Observational Implications: a new redshift-dependent Fit for the IRX-β Dust Attenuation Relation
The bulk of the galaxies modeled in this work have IRX and β that agree well with the relation presented in Overzier et al. (2010). We show this in Figure 9, where the mean offset between TNG50 galaxies across redshifts and the Overzier et al. (2010) relation for local galaxies read approximately 0.15 dex in IRX. This figure shows histograms of the differences between each galaxy's log 10 IRX predicted by our model and the log 10 IRX Overzier obtained by applying the fit presented in Overzier et al. (2010) for the same β. Analog histograms for the fits presented by Meurer et al. (1999) and Casey et al. (2014) (not shown) would give a worse agreement, with an average offset of 0.3 and 0.33 dex in IRX, respectively, at z = 0. Especially at z 0.5 the modeled IRX as a function of β is close to the IRX one would derive from the relation presented in Overzier et al. (2010). The left hand side of the figure shows the histogram of ∆IRX = log 10 IRX − log 10 IRX Overzier for discrete redshift bins. The peak of the distribution moves toward higher ∆IRX as we increase redshift. While most z = 0 galaxies pile up around ∆IRX = −0.15, the z = 4 galaxies have their peak at ∆IRX = 0.79, implying that there is a redshift evolution in the IRX-β relation (driven by variations in the intrinsic UV-slope as a function of redshift, as discussed in Section 3.2). The fraction of galaxies that deviates from the Overzier et al. fit increases with increasing redshift, up to 50% at z = 3 and even 90% at z = 4.
To account for the redshift trend, we fit a modified version of the Overzier et al. (2010) relation to our data at each redshift. To this end, we add a parameter β z , as follows: log 10 IRX z = log 10 1.68 + log 10 (10 0.4(3.85+1.96(β+β z )) − 1) (6) The best fit β z at each redshift are listed in Table 2. Note that at z = 0 β z is not zero, which means we find a slightly different fit from Overzier et al. for local galaxies. We find that β z can then well be described as a linear function of redshift, such that: In the right panel of Figure 9, we show the distribution of ∆IRX z (best fit) = log 10 IRX − log 10 IRX z (best fit) at each redshift. The histograms are now centered on ∆IRX z (best fit) = 0. After including β z , the fraction of galaxies with ∆IRX z (best fit) > 0.5 has dropped to ∼ 2 − 3% at all redshifts (see Table 2).
The presented fit offers a direct approach for observers to improve the usage of the IRX-β dust attenuation relation. We do note that the presented relation should be further tested observationally at redshifts up to z = 4, especially since in this work we assumed a uniform Milky Way type dust distribution. Beyond that, it is worthwhile to explore the validity of the presented relation at z > 4.

The Impact of Dust Composition
The fiducial model of this work assigns Zubko et al. (2004) Milky Way type galaxy dust to all galaxies. However, in reality it is plausible for dust types to be different across different galaxy types and across redshifts. It has been suggested, for example, that highredshift galaxies have dust more similar to SMC type dust, rather than Milky Way type dust (e.g., Capak et al. 2015). To further explore the effects of different dust types, we assign dust as it is present in the LMC and SMC to a subset of galaxies at z = 0.5 and z = 4 and compare the IRX-β distribution to the same subset of galaxies with Milky Way type dust (Figure 10).
LMC and SMC dust types universally decrease the slope of the IRX-β relation, in comparison to the slope with Milky Way dust. The effect is stronger for a SMC dust composition. This is because of different amounts of far-UV extinction (FUV extinction) in the characteristic extinction curves for the dust types. LMC dust -and even more so SMC dust -have a higher FUV to NUV extinction ratio than Milky Way dust. This leads to a quicker change in UVslope β as the dust opacity is increased. Therefore, the galaxies do not follow the typical Milky Way type dust relations any longer, but have a "redder" β for the same IRX. When employing SMC dust, the bulk of our subset of galaxies at z = 0.5 follow the relation found by Pettini et al. (1998) for SMC galaxies.
We also see a change in the IRX-β relation for SMC dust galaxies at z = 4 when varying the dust type. Unlike for z = 0.5, the relation does not follow the Pettini et al. (1998) trend, but is actually close to the z = 0 reference relation for Milky Way type galaxies. Because the intrinsic UV-slope of the z = 4 population is so blue, this balances out the "extra" reddening due to SMC dust (compared to Milky Way dust) and moves the galaxies towards the z = 0 reference relation for Milky Way type dust, rather than the Pettini et al. (1998) relation for SMC dust. This demonstrates a degeneracy between the effects that stellar population age and different dust types have on the location of galaxies in the IRX-β plane.
Some observations have found that z = 2 − 3 galaxies generally agree with the low redshift reference relations (or lie towards slightly redder β) and found no trend with redshift contrary to what is suggested by our models (Reddy et al. 2012;Bouwens et al. 2016). If the dust composition of these galaxies is indeed different from Milky Way type dust, the change in β due to dust composition Number of Galaxies TNG50 SFMS galaxies M * > 10 9 M MW dust z = 0 z = 0.5 z = 1 z = 2 z = 3 z = 4 Figure 9. Left: A probability density histogram of ∆IRX, which is the difference between log 10 IRX as obtained from our radiative transfer and TNG50 calculations and log 10 IRX Overzier as obtained from the Overzier's relation (Overzier et al. 2010), given β, binned as a function of redshift. The peak of the distribution moves to higher ∆IRX as we increase redshift, indicating that there is a redshift evolution in the IRX-β relation. Right: We fit to our data at fixed redshift a modified version of the Overzier's relation, where we introduce a new parameter β z in Equation 6 in order to account for the systematic redshift trend. The panel shows the histogram of ∆IRX z (best fit), which is now the difference between the predicted log 10 IRX at fixed β, and the best-fit log 10 IRX z . The redshift evolution is now accounted for in the new formula with all distributions peaking at 0.  Casey et al. 2014Meurer et al. 1999Overzier et al. 2011 SMC: Pettini+ 1998 z = 0.5 z = 4 Figure 10. IRX-β for a random subset of galaxies at z = 0.5 and at z = 4 for different assumptions of dust type. From left to right: Milky Way dust (as used throughout as fiducial choice) LMC, and SMC dust. As we vary the dust type to higher FUV absorbing dust, the slope in the IRX-β relation flattens.
could balance the change in β due to a young stellar population age in such a way that the galaxies fall close to the z = 0 reference relation. If this scenario is true, it means that the use of the IRX-β plane to study the dust type in different galaxies is hampered by the age/dust type degeneracy. This degeneracy can be broken by either measuring the stellar population age to account for the intrinsic UV-slope, or by measuring the detailed dust properties of these galaxies.

Comparison with Previous Work
The IRX -β dust attenuation relation has been studied using hydrodynamical models in the literature before. In this subsection we discuss these efforts and contrast them against the results presented in this work. Safarzadeh et al. (2017) analyzed a set of 51 idealized hydrodynamical simulations of disk galaxies and mergers at z = 0 and z = 2 − 3, on which radiative transfer was performed in postprocessing. They found that at z = 0, galaxies are usually found close to the relation of Meurer et al. (1999), while at z = 2 − 3, they lie above it even though they are not necessarily starbursts. They explain these deviations by disassociated UV-and IR-emission in these galaxies. The dust type significantly impacts the slope of the IRX-β relation. Lastly, they find that ULIRGs will deviate strongly toward high IRX. Narayanan et al. (2018) performed zoom-in simulations of a number of galaxies down to z = 2 (and z = 0 in one case) and concluded that Milky Way like galaxies, with relatively young stellar populations and cospatial IR-and UV-emitting regions lie near the standard relations such as the Meurer et al. (1999) relation. They also find causes for substantial deviations: old stellar populations fall below the canonical relations, complex dust geometries lead to a deviation above the relations, and are common in high-redshift heavily star forming galaxies, SMC dust decreases the slope in IRX-β and ULIRGs generally fall significantly above the reference relations. Ma et al. (2019) analyzed 34 zoom-in simulations of z > 5 galaxies. They assumed SMC dust in their radiative transfer calculations, and confirmed a tight relation in agreement with the Pettini et al. (1998) relation, despite the patchy dust geometry they find in their modeled galaxies. They also report that higher redshift galaxies move toward bluer β at fixed UV due to younger stellar populations and also find that dust type affects the slope in the IRX-β distribution as well.
Our results are consistent with those obtained by previous theoretical efforts in the literature. However, this is the first time the IRX-β relation is investigated through theoretical models for a homogeneous and large unbiased sample of many thousands of galaxies that are fully representative of SFMS galaxies from z = 0 to z = 4 in a full cosmological context, albeit volume and mass limited. This strengthens the applicability of our results, particularly for the redshift evolution in the IRX-β relation due to stellar population age.

Limitations
Throughout the analysis, we have fixed the ISM dust type of each galaxy globally, which means that every gas cell of each whole galaxy contains the same type of dust (Milky Way dust). In reality, it is plausible to expect that galaxies consist of a mix of dust types in different regions, which change the resulting IRX-β compared to our results. It is also probable that the dust content does not scale linearly with the gas metallicity, which is assumed in our model. The dust-to-gas fraction may vary across different galaxy regions and over time. For example,Vogelsberger et al. 2019a obtain a strongly declining dust-to-metal ratio as a function of redshift by comparing the TNG dust-attenuated rest-frame UV luminosity functions to existing observational data: their dust-to-metal varies between almost 0.9 at z = 0 to 0.1 at z = 6. It is plausible that a varying dust-to-metal ratio would return a different redshift evolution of the IRX-β plane than the one put forward here (see Appendix A. In recent years, multiple groups have implemented the tracking of dust and the extinction curve of dust in hydrodynamical simulations (McKinnon et al. 2016(McKinnon et al. , 2017(McKinnon et al. , 2018Hou et al. 2019;Li et al. 2019;Vogelsberger et al. 2019b). These approaches are a promising avenue to avoid having to make assumptions about the type and amount of dust in modeled galaxies.
The effects of birth clouds on the attenuation of stellar emission has been modeled by employing MAPPINGS-III spectra for young stellar populations. This crude implementation of unresolved birth cloud attenuation limits the predictive power of our model, especially for high-redshift galaxies and starbursts that contain high fractions of young stars. A better treatment of birth clouds will require higher resolution simulations that resolve the ISM structure in these objects.
Finally, we caution the reader not to over interpret some cases of TNG50 objects that appear as strong outliers from the average IRX-β relation. As documented in Section 5.2 of Nelson et al. 2019b, not all entries in the SUBFIND catalogs should be considered 'galaxies', as they may result from the fragmentation or collapse of gas within already formed galaxies. The criterion adopted in this work (non-vanishing stellar mass and a total dark matter mass fraction of at least 20 per cent, see Section 2.2) might not be sufficient to exclude all these objects from the analysis and it is expected for them to appear as outliers in most relations among galaxy physical properties.

SUMMARY AND OUTLOOK
We have performed radiative transfer calculations with the publicly-available software SKIRT on star-forming galaxies at 0 z 4 taken from the TNG50 simulation. Our fiducial model adopts a universal Milky Way type dust and a dust-to-metal fraction of 0.3 throughout. It considers a volume-limited sample of galaxies above 10 9 M in stellar mass (more than ten thousand stellar particles) at six different redshifts (from z = 0 to z = 4), corresponding to a total of 7280 galaxies. We have therefore quantified the evolution of the IRX-β dust attenuation relation of galaxies and how it correlates with galaxy properties. We summarize the insight we gained from this work as follows: • The bulk of the TNG50 SFMS galaxies at z < 1 follows the reference relations for z ≈ 0 by Meurer et al. (1999), Overzier et al. (2010 and Casey et al. (2014). Only a small fraction ( 5 per cent) of the low-redshift TNG50 SFMS galaxies with stellar mass M * > 10 9 M deviate from the local relation found by Overzier et al. 2010 (with ∆IRX > 0.5). This justifies the use of these relations at these redshifts to account for dust obscuration when IR data is missing (see Figure 9 and Table 2.) • More generally, we find a systematic redshift dependent shift along β in the IRX-β plane, where higher redshift TNG50 SFMS galaxies tend to be shifted toward lower β. This is at first order caused by a lower median intrinsic UV-slope β 0 due to younger stellar population ages, which makes TNG50 galaxies at high redshifts poorly described by an Overzier et al. (2010) relation. We therefore derived a new redshift-dependent version of the relation of Overzier et al. (2010), such that: log 10 IRX z = log 10 1.68+log 10 (10 0.4(3.85+1.96(β+β z (z))) −1), (8) where β z (z) = 0.142z − 0.081. This fit describes the IRX of galaxies at 0 z 4 well, with only a few per cent of strong outliers (see Figs. 7,9 and Table 2).
• Several physical properties determine the location of galaxies on the IRX-β plane. We find that IRX increases with larger stellar mass, gas-phase metallicity, SFR and star formation efficiency. The IRX-β as a whole shifts towards bluer (i.e. more negative) β for increasing sSFR.
• Dust composition plays an important role for both the shape and evolution of the IRX-β relation. The IRX-β relation for LMC and SMC type dust (with a higher FUV-to-NUV extinction ratio than MW type dust) is flatter than it is for Milky Way type dust. In certain cases (we demonstrate this for z = 4), the effects of young stellar age (resulting in a blue intrinsic UV-slope) balance out the effects of SMC type dust (which increases β to redder values compared to a typical Milky Way type of dust), such that a galaxy may appear to follow the reference local IRX-β relation for galaxies with Milky Way type dust, despite having SMC type dust. Without additional knowledge on the stellar population properties of galaxies, this hampers the use of the IRX-β dust attenuation relation as a proxy for dust types in galaxies.
Our results, especially the redshift dependent IRX-β dust attenuation relation and the degeneracy between stellar age and dust type, have clear consequences for the interpretation of observations of z > 1 galaxies where IR information is lacking. We look forward to future observations testing the redshift dependent trend we found in our simulations and hopefully being able to discriminate between different dust types, ultimately testing the applicability of the IRX-β dust attenuation relation in z > 0 SFMS galaxies.

APPENDIX A: CHECKS ON THE FIDUCIAL CHOICES AND FIDUCIAL RESULTS
In this Appendix, we quantify how different choices on the underlying radiative transfer and dust models may affect our results. Tables A1 and A2 summarize our SKIRT fiducial setup and the adopted SKIRT input parameters in our fiducial model. Here, we comment upon or explicitly show the effects of changing selected and relevant choices from tests performed on a subset of 440 randomly selected galaxies at z = 0.5. We distinguish between physical and SKIRT implementation choices and comment on the possible effects of numerical resolution.

A1.1 Varying the Gas Dust Fraction
In the fiducial model, we set the dust-to-metal fraction of the gas to be f = 0.3 (see Eq. 2). Vogelsberger et al. (2019a) performed a calibration with the TNG simulations to infer f at each redshift via comparison to observed UV luminosity functions. They find f to vary between about 0.9 at z ∼ 0 to 0.1 at z ∼ 6. In the following, we test what effect a change of this parameter has on the IRX-β distribution of the TNG50 SFMS galaxies. In principle, higher values of f imply increased gas cell opacity in each galaxy, because each cell contains more dust. For a homogeneous dust screen, a higher opacity would lead simply to a shift along the respective local IRX-β reference relations, toward higher IRX, and higher β. Only if the dust geometry is very complex (dust and young stellar populations are not cospatial) would a change in f result in a more complex track of galaxies in the IRX-β plane (e.g., see the explanation in Casey et al. 2014 andPopping et al. 2017b). We vary the dust-tometal fraction for a random subset of 440 galaxies at z = 0.5 in our sample and see how this affects their IRX and β values ( Figure A1, top row). The difference between f = 0.3 and f = 0.1 (red symbols) is obviously smaller than between f = 0.3 and f = 0.9 (blue symbols): as the dust-to-metal fraction increases, log 10 IRX and β increase as well, with an upward shift of about 0.5 and 0.3 on average, respectively. Still, we verified that the TNG50 galaxies are still roughly consistent with the z ∼ 0 observed IRX-β dust attenuation relations.

A1.2 Self absorption
During the SKIRT setup, it is possible to turn on dust self absorption. When this option is turned on, the Monte-Carlo simulation will take into account absorption and re-emission of dust that has previously been re-emitted by dust itself. This effect only plays a role where very thick layers of dust are present, but significantly increases computation time. For this reason, in our fiducial model we have left this option turned off. The test with a subset of z = 0.5 galaxies (not shown) confirms the bounty of our choice, as with or without dust self absorption results on log 10 IRX and β are uncertain to the percent level and below. Only for a fraction of massive galaxies the value of log 10 IRX can differ by up to 5 per cent or so.

A1.3 MAPPINGS-III Parameters
In our fiducial setup, young stellar populations of ages younger than 10 Myr are assigned MAPPINGS-III spectra. A number of choices are needed to implement this, but they are overall all less relevant than the effect of using MAPPINGS-III spectra itself.
In Figure A1, middle row, we show the effects of turning off MAPPINGS-III completely and of only assigning Bruzual-Charlot spectra to the stellar particles regardless of their age. The impact is not negligible, at least at z = 0.5 as shown. If we exclude MAPPINGS-III from our radiative transfer calculations, we end up with on average lower IRX (a shift of about 0.1-0.2 in the log) and lower β (an average shift of 0.3-0.35), because the effects of birth clouds are neglected. The shift in IRX is more relevant (up to 0.25 dex) for lower-mass galaxies. We keep MAPPINGS-III spectra included in our fiducial model, to be more realistic and comparable to other works that employed SKIRT radiative transfer calculations to IllustrisTNG.
Among the physical choices needed to implement the MAPPINGS-III emission, the interstellar medium pressure has no significant effect on the results (smaller than per cent level shifts: not shown). The compactness parameter C is directly influenced by P ISM via the equation: where M cl is the star cluster mass. As explained in section 2, stellar particles have masses of around 10 5 M in the highest resolution variant of TNG50. Therefore the compactness derives itself directly from the interstellar medium pressure via: However, these choices are not impactful in any manner that is relevant for the study at hand. We therefore choose as fiducial the value specified in the literature of X = log 10 P ISM /k B cm −3 K = 5 for our fiducial model.
On the other hand, an increase of the PDR covering fraction to f PDR = 0.9 would not be negligible. In Figure A1 this effect is quantified by comparing the fiducial choice of f PDR = 0.2 with f PDR = 0.01 (red symbols) and 0.9 (blue symbols) for a subset of test galaxies at z = 0.5. A higher covering fraction of birth clouds would increase the absorption of UV-emission by the young stellar populations (effectively increasing β) and in turn increase the infrared emission. This only happens, however, when we choose extreme values such as f PDR = 0.9, and therefore we maintain that the qualitative analysis of our results is not affected by this.
A2 SKIRT implementation choices SKIRT offers various choices for e.g. the selection of the wavelength grid which is used for sampling the spectra and the number of photon packages for the radiative transfer.
Throughout, we chose a custom wavelength grid with 59 wavelengths, with high detail in the UV and IR ranges. To test if the results given by this choice are converged, we increase the number of wavelengths to 93 wavelengths and see if the IRX and β results change. Differences are minimal across the mass range (less than a few percent for log 10 IRX and less than 10 percent for β, respectively), giving us confidence that the IRX-β results are converged already with 59 wavelengths.
To test if our choice of number of photon packages grants convergence, we increase the photon package count tenfold, and rerun the radiative transfer on a subset of 267 galaxies at z = 0.5. Differences are smaller than at the per cent level, implying that our fiducial choice is converged. In all panels, we show the shift between the test variations and the fiducial choice, and here we show only the choices that we have found to have a non negligible impact in the quantitative results. Top: Effects of varying the dust-to-metal ratio from the fiducial value of 0.2 to 0.1 (red) and 0.9 (blue data points). Middle: Effects of excluding the MAPPINGS-III spectra for young stellar populations in the radiative transfer model from the fiducial setup. Bottom: Within the MAPPINGS-III implementation, effects of varying the dust covering fraction of the PDR regions from 0.2 (fiducial choice) to 0.01 (red) or 0.9 (blue data points). Turned off Turned off Turned off Wavelength grid 59 points in 0.1 µm < λ < 1001 µm 1168 points in 0.05 µm < λ < 5 µm 50 points in 0.35 µm < λ < 0.95 µm Photon package number 10 6 equal to the total number of bound stellar particles, with 10 2 < Np < 10 5 50 per data cube pixel Observing instrument Viewing distance 10 Mpc 1M pc 10 Mpc Field of view 15r 1/2 variable 15r 1/2 Viewing angle Along z-axis Along z-axis Along z-axis Output SED and 100 × 100 cata cube SED and data cube SED and data cube with pixel scale 0.25arcsec/pixel   Table A3. A summary of the most important numerical specifications of the highest resolution TNG50 variant (TNG50-1), which is used in the main body of this work, and of a low resolution TNG50 variant used in the resolution test in Appendix A3 (TNG50-3). L z=0 box denotes the proper length of the simulation box at z = 0. N gas is the number of gas cells at the start of the simulation run, which can change over the course of the simulation due to star formation or refinement of grid cells. N DM is the number of dark matter particles, which stays constant over the whole simulation. m baryon is the mean gas cell mass, and also roughly the initial mass of stellar particles, while m DM is the mass of the dark matter particles (all DM-particles have the same mass). The gravitational softening length for dark matter particles at z = 0 is given by z=0 DM, stars . Softening lengths for gas cells are adaptive, with a minimum value of min gas (see e.g. Fig.1 of Pillepich et al. 2019).
Finally, in order to perform the SKIRT simulation, we need to calculate adaptive smoothing scales for the stellar particles. In our fiducial setup we have chosen the smoothing scale of each particle to be the distance to the 32nd closest neighboring stellar particle. We tested if the results change when opting for smoothing scales that correspond to the distance to the 16th and 64th closest neighbors. We find no noticeable effects on the derived IRX and β values (at per cent levels).

A3 Resolution Effects
To uncover possible effects of numerical resolution, we qualitatively compare the IRX-β distributions of our flagship simulation, TNG50, with one with much worse resolution. For this purposes, we use a low-resolution variant of TNG50, namely TNG50-3 (see Table A3): this has only 2×540 resolution elements inside the same comoving volume. The mass resolution of the baryonic resolution elements is therefore 54 × 10 5 M and for dark matter resolution elements it is 290 × 10 5 M : 64 times worse than TNG50. In Figure A2, we compare the IRX-β distribution of all selected z = 0.5 galaxies of the high resolution variant TNG50 (aka TNG50-1) with the low resolution variant TNG50-3. A large number of galaxies in TNG50-3 fall below the reference relations. This is due to the fact that galaxies with masses of 10 9 M are not well resolved, as they only consist of a few hundreds of stellar particles in total. On the right hand side of the figure, we therefore experiment by excluding TNG50-3 galaxies with stellar masses below 10 10 M , so as to only include reasonably resolved galaxies. When including only galaxies of TNG50-3 that are similarly resolved as in TNG50 (with a few thousand stellar particles each), the reference relation is recovered, as in TNG50. This underlines the importance of using well resolved galaxies (both in mass and spatial numerical resolution) but, conversely, also shows that within TNG50 and above the adopted minimum mass of 10 9 M (i.e. more than 10 4 stellar particles each) the results provided throughout are trustworthy.