Dust extinction, dust emission, and dust temperature in galaxies at z>=5: a view from the FIRE-2 simulations

We present a suite of 34 high-resolution cosmological zoom-in simulations consisting of thousands of halos up to M_halo~10^12 M_sun (M_star~10^10.5 M_sun) at z>=5 from the Feedback in Realistic Environments project. We post-process our simulations with a three-dimensional Monte Carlo dust radiative transfer code to study dust extinction, dust emission, and dust temperature within these simulated z>=5 galaxies. Our sample forms a tight correlation between infrared excess (IRX=F_IR/F_UV) and ultraviolet (UV)-continuum slope (beta_UV), despite the patchy, clumpy dust geometry shown in our simulations. We find that the IRX-beta_UV relation is mainly determined by the shape of the extinction curve and is independent of its normalization (set by the dust-to-gas ratio). The bolometric IR luminosity (L_IR) correlates with the intrinsic UV luminosity and the star formation rate (SFR) averaged over the past 10 Myr. We predict that at a given L_IR, the peak wavelength of the dust spectral energy distributions for z>=5 galaxies is smaller by a factor of 2 (due to higher dust temperatures on average) than at z=0. The higher dust temperatures are driven by higher specific SFRs and SFR surface densities with increasing redshift. We derive the galaxy UV luminosity functions (LFs) at z=5-10 from our simulations and confirm that a heavy attenuation is required to reproduce the observed bright-end UVLFs. We also predict the IRLFs and UV luminosity densities at z=5-10. We discuss the implications of our results on current and future observations probing dust attenuation and emission in z>=5 galaxies.


INTRODUCTION
Improving the constraints on the star formation rate density (SFRD) across cosmic time is important for understanding the assembly history of galaxies (see Madau & Dickinson 2014, for a recent review). At z 5 in particular, the cosmic SFRD directly relates to the number of ionizing photons available from star-forming galaxies for cosmic reionization (dependent upon the escape fraction, e.g. Finkelstein et al. 2012;Kuhlen & Faucher-Giguère 2012;Robertson et al. 2013Robertson et al. , 2015, so understanding the SFRD at z 5 is also crucial for constraining the reionization history. E-mail: xchma@berkeley.edu It is well known that at z 3, the cosmic SFRD is dominated by dusty star-forming galaxies (DSFGs; e.g. Magnelli et al. 2011;Gruppioni et al. 2013), which have very high star formation rates (SFRs) but are too heavily obscured to be seen in the UV and optical. On the other hand, the cosmic SFRD at z 5 is only probed in the rest-frame UV up to z ∼ 10 (e.g. Coe et al. 2013;Ellis et al. 2013;Oesch et al. 2013Oesch et al. , 2014Bouwens et al. 2015Bouwens et al. , 2016aFinkelstein et al. 2015). A consensus of the obscured fraction of star formation at these redshifts is not yet in place. The most commonly adopted approach to correct for dust obscuration in UV-selected galaxies at z 5 is to use the empirical relationship between infrared (IR) excess, IRX ≡ FIR/FUV, and the UV-continuum slope, βUV (e.g. Bouwens et al. 2015;Finkelstein et al. 2015). This so-called IRX-βUV relation was first established in local compact starburst galaxies (e.g. Meurer et al. 1999) and has been confirmed to hold up to z ∼ 2-3 (e.g. Reddy et al. 2006Reddy et al. , 2018Álvarez-Márquez et al. 2016;Bourne et al. 2017;Fudamoto et al. 2017;McLure et al. 2018).
However, at z 5, it is yet unclear whether the IRX-βUV relation, which is a reflection of the dust attenuation law, still applies. Capak et al. (2015) and subsequently Barisic et al. (2017) studied a sample of z ∼ 5.5 galaxies and found that they exhibit a large scatter in the IRX-βUV relation (see also Bourne et al. 2017 andFudamoto et al. 2017 for their z ∼ 5 sample; however, Koprowski et al. 2018 reported that z ∼ 5 galaxies are still consistent with the local relation), where some galaxies fall significantly below the IRX-βUV relation derived from a steeper, Small Magellanic Cloud (SMC)-like reddening law. Moreover, Bouwens et al. (2016b) showed that deep 1.2 mm continuum survey with the Atacama Large Millimeter Array (ALMA) in the Hubble Ultra Deep Field (HUDF) detects much fewer high-redshift sources than what inferred from the IRX-βUV relation using the rest-frame UV slopes and luminosities of galaxies in the same field (see also Dunlop et al. 2017), where a T ∼ 35 K modified black-body (MBB) function is assumed for dust emission. This also suggests that the IRX of high-redshift galaxies are well below the SMC IRX-βUV relation.
There are two possible explanations of these results: one physical, one observational. First, it is likely that z 5 galaxies are much less dust rich than low-and intermediate-redshift galaxies, because the universe has not allowed sufficient time for dust to grow substantially. On the other hand, it is possible that dust luminosity is severely underestimated due to the assumed spectral energy distributions (SEDs) of dust emission. As noted by Bouwens et al. (2016b), the tension can be alleviated if dust temperature reaches as high as 45-50 K at z ∼ 6 (see also e.g. Faisst et al. 2017), such that dust is much less luminous at long wavelengths at the same total IR luminosity. Alternatively, even if the cold dust remains ∼ 35 K, a moderate fraction of warm dust in the galaxy can dramatically reshape the dust SEDs and reduce the apparent flux density on the Rayleigh-Jeans (RJ) tail at fixed IR luminosity (Casey et al. 2018b). In either case, fundamentally different dust properties are not necessarily required for high-redshift galaxies. Therefore, it is critical to understand (1) in what conditions the local IRX-βUV relations still hold (or not) and (2) the typical dust temperature and SEDs in galaxies above z ∼ 5 to properly account for obscured star formation from pure rest-frame UV surveys.
Current knowledge on the population of DSFGs at z 5 is still limited. At the extremely luminous, high-SFR end (LIR ∼ 10 13 L ), there is a growing sample of DSFGs at z ∼ 5-7 built over the past few years (e.g. Riechers et al. 2013;Vieira et al. 2013;Weiß et al. 2013;Strandet et al. 2016;Marrone et al. 2018). Dust emission has also been detected in a small number of less extreme systems even at higher redshifts (e.g. Watson et al. 2015, z ∼ 7.5;Laporte et al. 2017, z ∼ 8.38;Hashimoto et al. 2018, z ∼ 7.15;Tamura et al. 2018, z ∼ 8.3), many of which are gravitationally lensed galaxies. It is obvious from these observational facts that dust plays an nonnegligible role in normal star-forming galaxy even in the very early Universe. However, there is still a lack of efficient ways for finding large samples of DSFGs at moderate luminosities (e.g. LIR ∼ 10 11 -10 12 L ) at z ∼ 5 and beyond.
Future observational facilities have been proposed, including the Chajnantor Sub/Millimeter Survey Telescope (CSST; Golwala 2018), the next-generation Very Large Array (ngVLA; McKinnon et al. 2016), the TolTEC camera on the Large Millimeter Tele-scope (LMT; Bryan 2018), and the Origins Space Telescope (OST; Battersby et al. 2018). Together with ALMA, these facilities are expected to advance our knowledge on high-redshift DSFGs in greater detail. To maximize the efficiency of future observations, it is of great importance to make useful predictions on the properties of DSFGs at z 5. For example, what are the IR luminosity functions (LFs; i.e. the number density of galaxies at a given LIR) at these redshifts (e.g. Casey et al. 2018a)? Is it reliable to infer dust luminosities from the UV-continuum slopes (e.g. as a way of finding targets to observe at longer wavelengths)? What wavelengths are most useful to probe to robustly measure LIR and dust temperature as well as to avoid severe contamination from the Cosmic Microwave Background (CMB; da Cunha et al. 2013)? Are there independent observables in the UV that can be combined with IR data to better constrain the dust properties at z 5?
From a theoretical point of view, it has been broadly appreciated that dust attenuation plays a key role in shaping the bright-end UVLFs at z 5. Numerous studies have demonstrated this both in semi-analytic models (SAMs; e.g. Clay et al. 2015;Liu et al. 2016;Cowley et al. 2018;Yung et al. 2018) and in cosmological simulations (e.g. Cullen et al. 2017;Wilkins et al. 2017;Ma et al. 2018b). Many of these previous studies have assumed simple prescriptions for dust extinction. For example, Ma et al. (2018b) calculated the integrated optical depth along a given sightline from every star particle in the simulations and applied an attenuation e −τ for individual particles to obtain the post-extinction UV luminosity (see also e.g. Katz et al. 2018). Such simple treatment does not properly account for certain radiative transfer effects, such as dust scattering, which can be important in the UV and optical (e.g. Barrow et al. 2017). Moreover, dust temperature and SEDs cannot be self-consistently calculated from this approach (e.g. Wilkins et al. 2018), which limits the predictive power of these calculations. To this end, full dust radiative transfer calculations are necessary.
In recent years, post-processing dust radiative transfer calculations have been conducted on large-volume cosmological simulations, cosmological zoom-in simulations, and idealized simulations of disks and mergers to investigate a broad range of questions, including dust extinction and emission in local galaxies (e.g. Camps et al. 2016;Trayford et al. 2017), the physical origin of the IRX-βUV relation (e.g. Safarzadeh et al. 2017;Narayanan et al. 2018b), the effects of dust geometry on the reddening law (e.g. Narayanan et al. 2018a), and the empirical relation between far-IR/submm flux and molecular gas mass (e.g. Liang et al. 2018;Privon et al. 2018). More relevant to high-redshift galaxies, Behrens et al. (2018) carried on dust radiative transfer calculation on one galaxy from a high-resolution cosmological zoom-in simulation and found that they require an extremely low dust-to-metal ratio (0.08, as oppose to the canonical value of 0.4 in the local Universe; Dwek 1998) and high dust temperature (91 ± 23 K) in order to reproduce the SED of the z ∼ 8.38 dusty galaxy detected by Laporte et al. (2017). Their results suggest that galaxies at such high redshifts are likely dust poor and have very high dust temperature. Population-wise, Cen & Kimm (2014) have run dust radiative transfer calculations on a sample of 198 galaxies in a cosmological zoom-in simulation at z ∼ 7 and predicted dust luminosities, SEDs, and IRLF at z ∼ 7. They found that 60-90% of the starlight is re-emitted in the IR, with the peak wavelength of dust SED around 45-60 µm.
Note that large-volume cosmological simulations usually have mass resolution ∼ 10 6 M and spatial resolution ∼ 1 kpc. Even the zoom-in simulations discussed above are only able to resolve down to ∼ 30 pc. It should be noted that dust geometry (clumpiness and covering fraction) and relative distribution between dust and stars have dramatic effects on the effective dust attenuation law, even if dust properties are constant (e.g. Seon & Draine 2016;Narayanan et al. 2018a). Moderate-and low-resolution simulations sometimes adopt sub-resolution models to account for dust distribution on unresolved scales (e.g. Jonsson et al. 2010), especially the heavy obscuration of young stars from their birth clouds (e.g. Charlot & Fall 2000). These models introduce extra free parameters which the results can be sensitive to (e.g. Liang et al. 2018).
In this work, we present a suite of high-resolution cosmological zoom-in simulations of z 5 galaxies as part of the Feedback in Realistic Environments (FIRE) project 1 (Hopkins et al. 2014. These simulations cover a broad range of halo mass up to 10 12 M at z = 5-10. We adopt a mass resolution ∼ 7000 M or better and the typical spatial resolution in dense gas is ∼ 1 pc. They use the FIRE-2 models of the multi-phase interstellar medium (ISM), star formation, and stellar feedback, which explicitly resolve stars forming in birth clouds and feedback disrupting these clouds. Simulations run down to z ∼ 0 using these models have been shown to reproduce a variety of observables, including the properties of giant molecular clouds (GMCs) in the local Universe , and reference therein). In particular, the z 5 simulations are shown to produce reasonable stellar masshalo mass relation, SFR-stellar mass relation, stellar mass functions, UVLFs, and cosmic SFRD that are broadly consistent with most up-to-date observational constraints (Ma et al. 2018b).
By post-processing these simulations with Monte Carlo dust radiative transfer calculations (without the need of sub-resolution dust recipes), we will study dust extinction and emission in z 5 galaxies that would be detectable in wide-field deep surveys in the rest-frame UV. Our work is built upon previous theoretical studies on dusty galaxies at z 5 by expanding the sample size, using more detailed simulations and post-processing methods, and broadening the scope. The paper is organized as follows. We briefly describe our simulation sample and the baryonic physics used in these simulations in Section 2.1. In Section 2.3, we describe the radiative transfer calculations. Section 3 mainly focuses on the UV and IR properties of dusty galaxies at z 5, where we investigate the IRX-βUV relation in Section 3.2, the bolometric luminosity of dust emission and its correlation with star formation activities in Section 3.3, and dust SEDs and dust temperature in Section 3.4. Section 4 focuses on predicting galaxy UV and bolometric IR LFs and cosmic SFRD at z = 5-10. We discuss the strategies for probing dusty z 5 galaxies and the limitations of this work in Section 5. We conclude in Section 6.

The simulations
This work uses a suite of 34 high-resolution cosmological zoomin simulations at z 5. The zoom-in regions are centered around halos randomly selected at desired mass and redshift from a set of dark matter (DM)-only cosmological boxes with periodic boundary conditions. The initial conditions are generated at z = 99 following the method in Oñorbe et al. (2014) using the MUSIC code (Hahn & Abel 2011), which uses the well-developed multi-scale cosmological zoom-in technique (e.g. Katz & White 1993;Bertschinger 2001). We ensure zero contamination from low-resolution particles within 2Rvir of the central halo, and less than 1% contamination in 3Rvir. 22 zoom-in regions are selected from a (30 h −1 Mpc) 3 box run to z = 5 around halos in Mhalo ∼ 10 9.5 -10 12 M , among which 15 are first presented in Ma et al. (2018b) and 7 more are added to improve the statistics at the high-mass end. Another 6 zoomin regions are selected from a (120 h −1 Mpc) 3 box run to z = 7 and the rest 6 from an independent box with the same size run to z = 9. They are centered on relatively more massive halos from Mhalo ∼ 10 11 -10 12 M at z = 7 and z = 9, respectively.
The initial mass for baryonic particles (gas and stars) ranges from m b = 100-7000 M , and high-resolution DM particles from mDM = 650-4 × 10 4 M , increasing with the mass of the central halo. Force softening for gas particles is adaptive, with a minimum Plummer-equivalent force softening length gas = 0.14-0.42 pc. Force softening lengths for star particles and high-resolution DM particles are fixed at star = 5 gas = 0.7-2.1 pc and DM = 10-42 pc, respectively. The softening lengths are in comoving units at z > 9 and in physical units thereafter. In Table 1, we provide the final redshift, mass resolution, force softening lengths, final halo mass, stellar mass, and selected galaxy properties of the central halo for all 34 zoom-in simulations. We explicitly check and confirm that there is no systematic difference between galaxies of similar masses but simulated at different resolution in all of our results in this paper (examples shown in Appendix A).
All simulation are run using an identical version of the code GIZMO 2 (Hopkins 2015) in the meshless finite-mass (MFM) mode with the FIRE-2 models of the multi-phase ISM, star formation, and stellar feedback , which we briefly summarize here. Gas follows an ionized+atomic+molecular cooling curve in 10-10 10 K, including metallicity-dependent fine-structure and molecular cooling at low temperatures and high-temperature metal-line cooling for 11 separately tracked species (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe). At each timestep, the ionization states and cooling rates for H and He are computed following Katz et al. (1996) and cooling rates from heavier elements are calculated from a compilation of CLOUDY runs (Ferland et al. 2013), applying a uniform, redshift-dependent ionizing background from Faucher-Giguère et al. (2009) and an approximate model for H II regions generated by local sources. Gas self-shielding is accounted for with a local Jeans-length approximation.
Star formation is only allowed in dense, molecular, and locally self-gravitating regions with hydrogen number density above nth = 1000 cm −3 (Hopkins et al. 2013). Each star particles is treated as a stellar population with known mass, age, and metallicity assuming a Kroupa (2002) initial mass function (IMF) from 0.1-100 M . The simulations account for the following feedback mechanisms: (1) local and long-range radiation pressure, (2) photoionization and photoelectric heating, and (3) energy, momentum, mass, and metal injection from supernovae (SNe) and stellar winds. The luminosity, mass loss rates, and Type-II SNe rates of each star particle are obtained from STARBURST99 (Leitherer et al. 1999), and Type-Ia SNe rates following Mannucci et al. (2006). Metal yields from Type-II and Ia SNe and AGB winds are taken from Nomoto et al. (2006), Iwamoto et al. (1999), and Izzard et al. (2004), respectively. Parameters describing the initial conditions and final galaxy properties of our simulations: (1) z final : The redshift which the zoom-in region is selected at and the simulation is run to.
(2) M halo : Halo mass of the central halo at z final .
(3) m b and m DM : Initial baryonic and DM particle mass in the high-resolution region. The masses of DM particles are fixed throughout the simulation. The masses of baryonic (gas and stars) particles are allowed to vary within a factor of two owing to mass loss and mass return due to stellar evolution.
(4) gas and DM : Plummer-equivalent force softening lengths for gas and DM particles, in comoving units above z = 9 and physical units thereafter.
Force softening for gas is adaptive ( gas is the minimum softening length). Force softening length for star particles is star = 5 gas.
(5) M * and M dust : Total stellar and dust mass within the virial radius, assuming a constant dust-to-metal ratio of f dust = 0.4 in gas below 10 6 K and no dust in gas hotter than 10 6 K.
(6) L bol and L IR : Bolometric and dust IR luminosity integrated from 0.08-1000 µm, accounting for all light coming out from the virial radius.
(7) SFR 10 and SFR 100 : Star formation rate averaged over the past 10 and 100 Myr, respectively, measured within the viral radius.
All simulations 3 are run with a sub-resolution turbulent metal diffusion algorithm described in Su et al. (2017) and Escala et al. (2018). We do not account for primordial chemistry nor Pop III star formation, but assume an initial metallicity of Z = 10 −4 Z .
We use the Amiga's halo finder (AHF; Knollmann & Knebe 2009) to identify halos and galaxies in the snapshots, applying the redshift-dependent virial parameter from Bryan & Norman (1998). There are more than one halo in each zoom-in region. In this paper, we restrict our analysis to halos that contain more than 10 4 particles and have zero contamination from low-resolution particles within Rvir to ensure good resolution. We also exclude subhalos from our dex in our simulation sample at z = 5, 7, and 9. These halos will be used to derive rest-frame UVLFs in Section 4.1. Dust radiative transfer calculations are conducted in all halos more massive than 10 10 M and a small number of halos above 10 9.5 M .

Basic properties of the simulated galaxies
In Fig. 2, we present scaling relations for our simulated galaxies at integer redshifts from z = 5-12. Each point represents one galaxy snapshot in our sample, color-coded by its redshift. The top-left panel shows the stellar mass-halo mass relation, where we use the total stellar mass within Rvir. The dashed line shows the best-fit linear relation log M * = 1.53 (log Mhalo − 10) + 7.40. The dotted line shows the linear fit from Ma et al. (2018b). Note that they measure stellar mass within Rmax/3 to exclude satellite galaxies and diffuse stars, where Rmax is the halo maximum velocity radius and Rmax/3 is roughly comparable to 0.2 Rvir in these halos. There is thus a 0.2-0.3 dex difference between the two relations. The redshift evolution of the M * -Mhalo relation is not significant (by less than 0.1 dex from z = 5 to z = 12, still within the scatter of the sample).
The top-right panel shows the relation between dust mass and stellar mass, both measured within Rvir. Here we assume a constant dust-to-metal ratio of 0.4 (Dwek 1998) in gas below 10 6 K and no dust in hotter gas. Unsurprisingly, Mdust is proportional to M * (dust mass equals to 0.48% of the stellar mass), as all the dust is produced in stellar evolution processes (SNe and AGB stars) by assumption. The bottom-left panel presents the relation between the intrinsic (unobscured) rest-frame UV luminosity (LUV,intr ≡ λL λ at 1500 Å, measured within Rvir) and halo mass. At a given Mhalo (and M * ), LUV,intr increases by an order of magnitude from z = 5 to z = 12, because galaxies at higher redshifts have higher SFRs (see also Ma et al. 2018b). The bottom-right panel shows the relation between LUV,intr and bolometric luminosity Lbol of stellar continuum. We find a universal relation LUV,intr = 0.68 Lbol for our simulated sample with no discernible scatter, because young (UV-bright) stars also dominate the total luminosity. We will use these scaling relations to interpret our results in the rest of this paper.

Dust radiative transfer
We post-process our simulations with the public three-dimensional Monte Carlo dust radiative transfer code SKIRT 4 (Baes et al. 2011;  There is a linear correlation M dust = 0.0048 M * , as dust is produced by stars following our assumptions. Bottom left: The intrinsic (unobscured) UV luminosity-halo mass relation. At a given halo mass, L UV, intr increases with redshift by an order of magnitude from z = 5 to 12, as galaxies at higher redshift tend to have higher SFRs. Bottom right: Intrinsic UV luminositybolometric luminosity relation. L UV, intr is a good proxy of L bol with little scatter (L UV, intr = 0.68 L bol ), as young stars dominate both the UV luminosity and the total luminosity. We will use these relations to interpret the results in the rest of this paper.
Camps & Baes 2015) to calculate galaxy continuous spectral energy distributions (SEDs) on a 90-point wavelength grid equally spaced in logarithmic scale from 0.08-1000 µm. For each halo, we include all gas and star particles out to Rvir in our calculations and compute galaxy SEDs and mock images at each wavelength along five random lines of sight.
We do not explicitly model dust formation, growth, and destruction in our simulations, but simply assume a constant dust-tometal ratio (Mdust = fdust Mmetal) in gas below 10 6 K and no dust in hotter gas. The dust grid is reconstructed from gas particles using the built-in octree grid in SKIRT, which includes all particles in a cubic domain with a side length of 2Rvir and adaptively refines the high-density region until the following criteria are met: (1) the dust mass in a cell does not exceed 10 −6 of the total dust mass in the domain and (2) the 15 th refinement level has reached (i.e. the cell size is 2 −15 of the domain size). The minimum cell width is less than 3 pc even for the most massive galaxy in our sample. We use 10 6 photon packets at each of the 90 wavelengths. These choices ensure excellent convergence at a reasonable computational cost. We refer to Appendix A for details about the convergence tests.
We adopt the Small Magellanic Cloud (SMC)-type dust grain size distribution from Weingartner & Draine (2001,   hydrocarbon (PAH) included. We also compare Milky Way (MW)type dust in Appendix A. In Fig. 3, we show the dust opacity from UV to far-IR. Dust absorption and scattering are both important in the UV and optical. At long wavelengths (λ > 30 µm), dust absorption dominates the opacity, which scales with wavelength (frequency) roughly as κ ∝ λ −β (ν β ) with dust emissivity spectral index β = 2. Our fiducial dust-to-metal ratio is fdust = 0.4 (Dwek 1998), which gives a gas opacity 5 at 1500 Å κ 1500 Å = 0.73 × 10 3 cm 2 g −1 fdust 0.4 assuming a solar metallicity Z = 0.02. In this work, we will also explore fdust = 0.2-0.8. We note that it is the absorption coefficient α ≡ κdust ρdust = κgas ρgas = κdust fdust Zgas ρgas that enters the radiative transfer equation and sets the dust temperature and emissivity (via Kirchhoff's law). There is a degeneracy between dust opacity and dust-to-metal ratio in the form of κdust fdust in these calculations. Our experiments with different fdust at fixed κdust should be understood as varying the normalization of the gas opacity. 6 That said, if we vary κdust and fix fdust, all radiative transfer results, including the intensity field and dust temperature, must be identical to our experiments here (varying fdust for fixed κdust). Photon packets are first launched from star particles and propagated in the domain until absorbed or escaped. The SEDs of star particles are calculated from the built-in STARBURST99 stellar population models in SKIRT (nearly identical to those used in our simulations), which are compiled by Jonsson et al. (2010) and include both stellar continuum and Balmer continuum from nebular emission. Dust temperature and emissivity are determined from the local intensity field assuming energy balance. Next, photon packets representing dust emission are launched and propagated in the domain. 5 By definition, the relation between gas opacity due to dust extinction and dust grain opacity is simply κgas = κ dust ρ dust /ρgas. 6 Note that the SMC gas opacity at 1500 Å in Pei (1992) is approximately 154 cm 2 g −1 , a factor of two difference from Equation 1 if using 0.1 Z for SMC metallicity (e.g. Pei 1992 suggested that the B-band gas opacities in the MW and SMC follow roughly 10:1). A factor of a few variation is also seen between different lines of sight. Our experiments with f dust = 0.2-0.8 account for the uncertainties of both dust opacity and dust-to-metal ratio.
The local radiation field and dust temperature are then updated to account for dust self-absorption. This step is done iteratively until the dust SED converges within 1%, when the calculation stops and a final solution is reached. In this work, we do not include heating from the CMB, but defer to a future study on its effects (see Section 5.2 for more discussion).
We use non-local thermal equilibrium (NLTE) dust emission enabled in SKIRT by coupling itself to the DUSTEM code (Compiègne et al. 2011). This accounts for emission from small grains that are transiently heated by individual photons. Moreover, grains of different sizes are no longer at a single equilibrium temperature, but follow a temperature distribution. The NLTE dust emission only affects the dust SED at wavelengths shorter than rest-frame 30 µm. Nonetheless, SKIRT still computes the equilibrium dust temperature for each cell (Teq) following where κ abs ν is the dust absorption opacity and Jν is the local radiation intensity at frequency ν. It is worth noting that the right-handside integral scales as T 4+β eq given κ abs ν ∝ ν β at long wavelengths where Bν (Teq) dominates. We will still use Teq to describe dust temperature in each cell, as the long-wavelength dust emission that we mainly focus on in this paper is not affected by NLTE effects.
We note that some previous works adopted the MAPPINGS III starburst SED models (Groves et al. 2008) for star particles younger than 10 Myr to account for unresolved small-scale dust distribution (e.g. Camps et al. 2016;Liang et al. 2018). These models describe the dynamic evolution of spherical H II regions on spatial scales of ∼ 5-800 pc around star clusters of mass 10 3.5 -10 7.5 M , which are born from at least 10 times more massive clouds. They also include dust extinction and emission from the photodissociation regions. Cosmological simulations at 10 5 M mass resolution or worse may use these models to account for sub-resolution dust distribution (see the discussion in section 2 of Jonsson et al. 2010). Our simulations, however, are able to resolve the mass and spatial scales at which the MAPPINGS III models describe. We thus do not use these models in our radiative transfer calculations, but take the dust distribution 'as such' in the simulations as an alternative to the MAPPINGS III models (see also Behrens et al. 2018). All the results in this paper are fairly converged at the resolution of our simulations (we show examples in Appendix A).

RESULTS: DUST EXTINCTION AND EMISSION IN
HIGH-REDSHIFT GALAXIES

Example images
The built-in 'peeling off' method (a next-event estimator) in SKIRT can produce high signal-to-noise images with a relatively small number of photon packets. Thanks to its Monte Carlo nature, SKIRT also allows us to separate light from different origins (e.g. from sources and dust, from direct transmission and scatter, etc.) using the so-called 'smart detectors' (Baes et al. 2011). In Fig. 4, we present mock images for two example galaxies, z5m12b at z = 5 (top) and z9m12a at z = 9 (bottom), which are the most massive galaxy in our sample at each redshift (using fdust = 0.4). The color scale in columns (a)-(d) includes 95% of the light/mass in the field of view. The white dashed circles show the Rmax/3 radius. Columns (a) and (b) show the rest-frame UV (1500 Å) images detected by a 'smart camera', decomposed into (a) light transmitted directly from stars and (b) light scattered by dust at least once . Example images of two galaxies in our sample, z5m12b at z = 5 (top) and z9m12a at z = 9 (bottom). From left to right: (a) rest-frame UV continuum (1500 Å) directly transmitted from stars, (b) UV continuum from dust scattering, (c) rest-frame 350µm dust continuum, (d) dust column density, and (e) dust effective temperature. The white dashed circles show the Rmax/3 radius. Dust distribution is patchy and extended to large radii. Scattered UV light contributes ∼ 1/3 of the post-extinction UV flux but is distributed over larger spatial scales and at lower surface brightness than the UV light direct from stars.
(i.e. adding (a) and (b) together gives the total UV flux as viewed by a regular camera). Most of the UV light is emitted from a compact region (less than 2 kpc in projected radius) in the center of the galaxy. Some part of the galaxy is heavily obscured by opticallythick dust patches along the line of sight, while other part is almost transparent. Scattered UV light is more spatially extended and at lower surface brightness (note the different color scales), making it more difficult to detect observationally than light from direct transmission (e.g. Ma et al. 2018a). For both galaxies, scattered light contributes 30-35% of the total post-extinction UV flux.
Columns (c) and (d) show the rest-frame 350µm dust continuum image and projected dust column density, respectively. Dust distribution is clumpy and patchy and extends to a much larger spatial scale than stars, owing to feedback-driven outflows pushing gas and dust to large radii. For a dust opacity of 446.5 and 38.1 cm 2 g −1 at 30 and 100 µm, respectively (Fig. 3), the mid-and far-IR is optically thin in most part of the galaxy except for the central dense region where the IR optical depth can reach order unity and dust self-absorption is thus important. Column (e) shows the dust effective temperature defined as Teff = ( ρT 4+β eq dl/ ρdl) 1/(4+β) with integration evaluated along the line of sight (see Section 3.4 for motives for this definition). In the very central region that is close to the sources producing most of the UV light, dust can be heated to over 45 K and a small fraction of dust (less than 0.1% in mass) even reaches up to 100 K. The diffuse dust out to 10 kpc is also heated by diffuse starlight to 20-30 K.

The IRX-βUV relation
The relationship between the infrared excess, IRX = FIR/FUV (FUV is the attenuated UV flux here), and the rest-frame UV continuum slope, βUV (to distinguish from the dust emissivity spectral index β above), where F λ ∼ λ β UV , is an empirical relation first established for local galaxies (Meurer et al. 1999) and being confirmed up to z ∼ 2 (e.g. Reddy et al. 2012). This relation is expected from a simple picture where an intrinsically blue source is obscured by a dust screen: as the amount of attenuation increases, the UV slope appears redder and the observed IR-to-UV flux ratio becomes larger. However, differential extinction between young and old stars (e.g. Charlot & Fall 2000) or clumpiness of dust distribution (e.g. Seon & Draine 2016) can dramatically alter the effective extinction curve (even the dust composition is fixed), which may result in large variations in the IRX-βUV relationship (e.g. Howell et al. 2010;Casey et al. 2014;Narayanan et al. 2018b).
At z > 3, it is not yet clear whether the IRX-βUV relation is still consistent with that for local galaxies (e.g. Koprowski et al. 2018;McLure et al. 2018), where a Calzetti-like attenuation curve should apply (Calzetti et al. 2000), or follows a shallower relation that is more consistent with SMC-like attenuation curve (e.g. Reddy et al. 2018), or there is no well-established IRX-βUV relation due to large scatter (e.g. Capak et al. 2015;Bouwens et al. 2016b;Barisic et al. 2017). It is also unclear whether such discrepancies in these observations are due to inconsistent measurements of the UV slope and IR luminosity (e.g. Gómez-Guijarro et al. 2018), selection biases in different samples, or intrinsic scatter driven by large variations of dust properties and complex dust geometry.
In the top panel of Fig. 5, we present the IRX-βUV relation for our simulated sample, color-coded by redshift (points; fdust = 0.4). FIR is the total dust flux integrated over 8-1000 µm and FUV is the neutral flux density, λF λ , at 1500 Å. βUV is measured using the monochromatic flux at two wavelengths 1500 Å and 2300 Å. We only show one sightline for each galaxy snapshot, but highlight all five lines of sight for galaxies z5m12b and z9m12a at z = 5 and 9, respectively (example images shown in Fig. 4), using blue and or- ASPECS-Pilot sample Capak+15 high-z sample 5 6 7 8 10 12 redshift Pattini+98 (SMC) Figure 5. Top: The IRX-β UV relation for our simulations (using f dust = 0.4). Each point represents a galaxy snapshot along a random line of sight, color-coded by redshift. Galaxies at higher redshifts tend to move slightly toward bluer β UV at fixed IRX due to their younger stellar populations. Points with errorbars show the observational datasets for high-redshift galaxies compiled in Casey et al. (2018b, consisting of two different samples). The lines show the empirical relations from the literature. The simulated sample forms a tight IRX-β UV relation broadly agrees with observations and the SMC IRX-β UV relation expected from a simple dust screen model. Our results suggest that patchy, complex dust distribution and nontrivial effect of dust scatter as shown in our simulations do not drive significant scatter in the IRX-β UV relation. Bottom: The IRX-β UV relation for different normalizations of the extinction curve (represented by varying f dust ; the color points). The grey squares show the results if dust scatter is ignored, which should be understood as changing the attenuation law. These results confirm that the IRX-β UV relation is determined by the shape of the extinction curve but independent from its normalization. The red squares show the results if the IR and UV fluxes are measured in a smaller aperture (Rmax/3 instead of R vir ), which has little effect on the results. ange circles to illustrate the variation from different viewing angles. We compare our results with the observational dataset compiled in Casey et al. (2018b), which consists of the ASPECS-Pilot sample from Aravena et al. (2016) and z ∼ 5.5 sample from Capak et al. (2015) with updated measurements by Barisic et al. and Casey et al. We also show the empirical IRX-βUV relation developed from local starburst galaxies in Meurer et al. (1999, dashed) and the aperturecorrected relation in Takeuchi et al. (2012, dotted). The solid lines show the IRX-βUV relation derived from simple dust screen model applying a SMC-like extinction curve (e.g. Pettini et al. 1998).
The simulated galaxies form a tight IRX-βUV relation. There is a small redshift evolution with galaxies moving toward bluer βUV at fixed IRX with increasing redshift, simply because of younger stellar populations at higher redshifts. Line-of-sight variations for the same galaxy follow the sample IRX-βUV relation. Our simulated sample broadly lies within the scatter of recent measurements and agrees well with the SMC IRX-βUV relation derived from a simple dust screen picture. Surprisingly, although our simulations show patchy, complex dust distribution in these galaxies and the radiative transfer calculations reveal non-trivial effects of dust scattering in the rest-frame UV (see Fig. 4), we suggest that they do not drive significant scatter in the IRX-βUV relation.
In the bottom panel of Fig. 5, we compare the IRX-βUV relation for different dust-to-metal ratios (color points). Note that this must be understood as varying the normalization of the extinction curve. We only show a subsample or our galaxy snapshots (at integer redshifts) along the same line of sight as in the top panel. We find that the IRX-βUV relation does not change with fdust. We also examine the situation where the UV and IR fluxes are measured using a smaller aperture (Rmax/3 instead of Rvir; red squares) and find it has no significant effect on the IRX-βUV relation. Finally, as a proof of concept, we show the results if dust scatter is ignored (grey squares), as changing the shape of the extinction curve. The sample is more consistent with the Calzetti-like IRX-βUV relation by coincidence, as the absorption opacity is shallower than the total opacity in the UV (see Fig. 4). Our results suggest that the IRX-βUV relation is mainly determined by the shape of the extinction curve and independent of its normalization, at least in the mass range we probe in our simulations. If the large scatter in the IRX-βUV relation reported in some z 5 galaxy sample is real, it is more likely to be caused by variations in the effective attenuation law, rather than by smaller dust-to-metal ratios in high-redshift galaxies.

Bolometric IR luminosity
In this section, we present the bolometric luminosity of dust emission (LIR, integrated over 8-1000 µm) for our simulated sample and its dependence on various galaxy properties. We include all the light within Rvir to calculate LIR, which is not a bad treatment as instruments probing these wavelengths usually have large beam sizes. These results are not only useful for understanding the physics of dust obscuration and emission in high-redshift galaxies, but also important for empirically modeling the abundances of dusty star forming galaxies at these redshifts.
At the bright end, the bolometric dust luminosity does not change with fdust by more than 0.1 dex, because most of the obscured sightlines are optically thick. At the faint end where the galaxies are optically thin to dust extinction, LIR is proportional to fdust. Galaxies at z = 5-12 lie on the same LIR-LUV,intr relationship. In the right panels of Fig. 6, we explore why there is no redshift dependence on the LIR-LUV,intr relation and what drives its scatter. To this end, we take all simulated galaxies in a narrow range of UV luminosity (within 0.1 dex from LUV,intr = 10 11 L , as labeled by the grey rectangular in the left panel of Fig. 6) and search for secondary dependence of LIR on total dust mass (Mdust, top left), average dust surface density ( Σ dust ≡ Mdust/(Rmax/3) 2 , top right), and average dust density ( ρ dust ≡ Mdust/(Rmax/3) 3 , bottom left). 7 Interestingly, all three quantities are redshift-dependent as expected: at fixed LUV,intr, galaxies at higher redshifts contain less dust mass 8 but show higher average dust column density (equivalent to optical depth) and density. However, none of these quantities correlate with LIR. This suggests that at a given LUV,intr, the dust luminosity 7 Here we adopt Rmax/3 as a characteristic size of dust distribution mainly for illustrative purposes. Other size measures, such as half-mass/half-light radius, are statistically scaled with Rmax up to a constant factor. 8 As shown in Fig. 2, L UV, intr increases with redshift at fixed halo mass and stellar mass from z = 5-12 (see also Ma et al. 2018b). Therefore, at fixed L UV, intr , galaxies at higher redshifts are less massive and thus less dust rich. is primarily determined by the covering fraction of optically-thick sightlines, regardless of total dust mass and dust density in the system.
The bottom right panel of Fig. 6 shows that the scatter of the LIR-LUV,intr relation is driven by the SFR averaged over the past 10 Myr, in other words, the amount of stars younger than 10 Myr in a galaxy 9 . Note that if the SFR is measured over longer timescale (e.g. 100 Myr), the secondary dependence of LIR on SFR at fixed LUV becomes weaker. The physical picture behind this result is that stars younger than 10 Myr are more heavily obscured by their birth cloud than relatively older stars (e.g. 10-100 Myr, which still contribute a significant fraction of the UV light). This is consistent with models where differential obscuration between young stars and older stars is applied by hand (e.g. Charlot & Fall 2000;Jonsson et al. 2010;Katz et al. 2018), albeit our simulations explicitly resolve this with our star formation and feedback models. Fig. 7 shows the LIR-SFR10Myr relation for the entire simulated sample. Each point represents one galaxy snapshot, color-coded by redshift. Again, the LIR-SFR10Myr relationship does not depend on redshift (as well as dust mass and density as we explicitly checked). This relation is best described by a single power-law function where ( The dependence on f dust is shown by the three red lines. Using galaxies at fixed SFR 10 Myr (0.1 dex from 10 M yr −1 , the grey shaded region), the smaller panel shows that the scatter in the L IR -SFR 10 Myr relation is driven by L UV, intr . This means that differential obscuration for young and relatively older stars is important for understanding dust extinction and emission.
the secondary dependence of LIR on LUV,intr at fixed SFR10Myr (within 0.1 dex from 10 M yr −1 , as labeled by the grey rectangular in the main panel). The red dashed line shows the double powerlaw fit in the left panel of Fig. 6. This is because stars older than 10 Myr provide an extra source for dust emission. Combining the results in Figs. 6 and 7 further confirms that differential obscuration between young and relatively old stars is important in understanding dust extinction and emission. Finally, we inspect the galaxies at both fixed LUV and SFR10Myr and find no further dependence of LIR on other dust properties: this is the intrinsic scatter purely due to variations of dust geometry in these galaxies.

Dust SEDs and dust temperature
In this section, we study the dust SEDs and dust temperatures for our simulated sample. Again, we include all the light in Rvir. As we show in Fig. 4, there is a broad distribution of dust temperatures in a single galaxy, with dust close to the young stars being heated up to 100 K and diffuse dust at large radii at much lower temperatures. It is thus non-trivial to parametrize dust SEDs and even define one dust 'temperature'. One of the most commonly adopted forms for modeling dust SEDs is the MBB function for single-temperature dust (e.g. da Cunha et al. 2013;Bouwens et al. 2016b) where the second expression is valid in the optically thin limit and a power-law opacity κν ∝ ν β is applied. In this situation, the peak wavelength of Lν 10 is λpeak = 96.64 µm (30 K/T ) and the total dust luminosity is ∝ T 4+β . A more realistic form is the two-component dust SED model, consisting of a MBB function for old dust and a power-law component for warmer dust . Nevertheless, an optically-thin MBB function at local equilibrium temperature is 10 Note that the peak wavelengths of Lν , L λ , and νLν are different. still a good approximation for the local dust emissivity at rest-frame λ > 30 µm where NLTE effects are negligible (see Section 2.3). We adopt three definitions of dust temperature that we will refer to in the discussion below. First, we define the peak temperature Tpeak = 30 K (96.64 µm/λpeak), which is motivated from the peak wavelength λpeak of Lν for an optically-thin MBB function. Next, we introduce the mass-weighted dust temperature where Teq is the equilibrium dust temperature given by SKIRT under the LTE assumption (see Equation 2). This is the most straightforward one to calculate from dust radiative transfer calculations for simulated galaxies and adopted by various authors in the literature (e.g. Behrens et al. 2018;Liang et al. 2019). It is worth noting that the mass-weighted temperature directly relates to the dust SED at the R-J tail (e.g. Scoville et al. 2016;Liang et al. 2019). Finally, we define the effective dust temperature Note that the frequency-integrated dust emissivity (power per unit volume) is jν dν = αν Bν (T ) dν = κν,dust ρdust Bν (T ) dν ∼ ν β ρdust Bν (T ) dν ∝ T 4+β ρdust ( jν is emissivity and should not be confused with the radiation intensity Jν in Equation 2), so the effective dust temperature is defined such that in the optically thin limit, the bolometric dust luminosity is LIR ∝ MdustT 4+β eff . 11 All three temperatures correlate with each other with large scatter depending on the exact dust temperature distribution in each galaxy. Casey et al. (2018a) suggested a redshift-independent, empirical relation between observed rest-frame peak wavelength λpeak of Lν and bolometric IR luminosity LIR derived from several observed samples from z = 0-6. These include data from the H-ATLAS survey mostly covering 0 < z < 0.5 (Valiante et al. 2016), the sample in the COSMOS field at 0.3 < z < 2 with Herschel detection (Lee et al. 2013), and the South Pole Telescope (SPT)-detected DSFGs sample with average redshift z ∼ 4.3 from Strandet et al. (2016). Both λpeak and LIR were remeasured by fitting the two-component dust SED model to the original data. In Fig. 8, we compile the individual galaxies in the H-ATLAS z < 0.1 sample (grey points), 1σ and 2σ ranges for the COSMOS sample (orange solid and dashed lines), the SPT z ∼ 4.3 DSFG sample (blue squares), and the bestfit power-law model In Fig. 8, we also present the λpeak-LIR relation for all galaxy snapshots at z > 5 from our simulated sample (color points; using fdust = 0.4). For a sanity check, we conduct dust radiative transfer calculations using identical methods on a sample of 12 Milky Way (MW)-mass galaxies from the FIRE simulations at z = 0, including 8 isolated halos and 2 Local Group (LG)-like galaxy pairs at mass resolution m b = 3500-7000 M 12 (comparable or better than those studied in this paper), run with the identical version of GIZMO. The λpeak-LIR relation for the z = 0 FIRE sample consisting of 12 MWmass galaxies is shown by the red triangles in Fig. 8. The z = 0 FIRE sample agrees well with the observed λpeak-LIR relation from the H-ATLAS z < 0.1 sample and lies along the empirical power-law model from Casey et al. (2018a, Equation 8). However, although the z 5 sample also shows an anti-correlation between λpeak and LIR, these galaxies offset from the observational data and the z = 0 simulations, with λpeak moving toward shorter wavelengths by a factor of 2 at a given LIR. At the most luminous end in our sample, we find λpeak ∼ 60-80 µm, in good agreement with previous simulations at similar redshifts post-processed with dust radiative transfer calculations (e.g. Cen & Kimm 2014). This suggests that dust is much warmer in z > 5 galaxies than in lowredshift galaxies. In fact, the effective dust temperature in our z = 0 MW-mass galaxy simulations is ∼ 18 K, whereas it is typically over 35 K in the z > 5 galaxies at similar IR luminosities (LIR = 10 10 -10 11 L ). Given that our z = 0 simulations in good agreement with observation, we argue that this prediction is a physical effect, as the 12 Among these simulations, 6 isolated halos and 2 LG-like pairs (10 galaxies) have been presented in Garrison-Kimmel et al. (2018). The other two isolated halos will be presented in  Figure 9. Rest-frame dust SEDs for our z > 5 simulations (using f dust = 0.4). In each panel, we show all galaxies brighter than L IR = 10 10 L with λ peak falling in a 0.05 dex bin centered on the wavelength marked by the grey arrow. All SEDs are renormalized to L IR = 10 11.5 L . The red lines in each panel show the optically-thin MBB function at T = 35, 45, and 60 K, also normalized to 10 11.5 L . The vertical cyan dashed lines illustrate the observed-frame 1.2 mm (ALMA Band 6) for z = 6 (rest-frame 171 µm). A 35 K MBB function overestimates the flux density at this wavelength by a factor of 3-10. To convert between L IR and observed-frame 1.2 mm flux density for z ∼ 6 galaxies using an optically MBB function, one must adopt a high dust temperature of 45-60 K. z = 0 and z > 5 simulation samples are run with the same code and comparable resolution. The black dashed lines in Fig. 8 show the best-fit power-law function of the λpeak-LIR relation for our sample, λpeak = 78.78 µm [(1 + z)/7] −0.34 (LIR/10 10 L ) −0.084 , at z = 6, 8, and 10 (see Equation 9 and Table 2 for details).
The SPT-detected DSFG sample at z ∼ 4.3 seems to lie on the same λpeak-LIR relation as low-redshift galaxies. This does not necessarily mean that the λpeak-LIR relation is redshift independent out to z 5. First, the SPT-detected galaxies are much more luminous than our simulated galaxies and λpeak is measured in different ways between the two samples. More important, galaxies of similar LIR but shorter λpeak have weaker flux densities at long wavelengths where the observations are conducted, so they tend to be excluded in a flux-limited sample. Therefore, the SPT sample cannot falsify our prediction that most z 5 galaxies have a factor of 2 shorter λpeak than lower-redshift galaxies.
The black arrows in Fig. 8 indicate the amount and direction that faintest and brightest galaxies in our sample move along on the λpeak-LIR plane, respectively, if fdust increases by a factor of 2 (i.e. fdust = 0.8). At the faint end, the total dust luminosity increases by a factor of 2 (see Fig. 6) while the dust mass also doubles. We thus expect the dust effective temperature (so does the peak temperature or λpeak) remains unchanged, so faint galaxies move horizontally to higher LIR by approximate 0.3 dex. At the bright end, the total dust luminosity changes very little, so the effective temperature should decrease by a factor of 2 1/(4+β) = 1.12. Therefore, bright galaxies move vertically toward longer λpeak by nearly 0.05 dex. This brings our z 5 sample closer to the observed λpeak-LIR relation, although a large offset still remains. Following a similar argument, galaxies will move along the opposite direction by the same distance shown by the arrows if fdust decreases by a factor of 2. This is confirmed by our radiative transfer calculations. T ∼ Σ 1/6 SFR 5 6 7 8 10 12 redshift Figure 10. The correlation between dust temperature and bolometric IR luminosity (left), specific star formation rate (middle, averaged over the past 10 Myr), and SFR surface density ( Σ SFR ≡ SFR 10 Myr /(Rmax/3) 2 , right). Each row represents one definition of dust temperature, including peak temperature (top), mass-weighted temperature (middle), and effective temperature (bottom). We only show snapshots at integer redshifts, color-coded by redshift ( f dust = 0.4). The dust temperature-sSFR correlation can be understood as T ∼ (L IR /M dust ) 1/6 ∼ (SFR/M * ) 1/6 ≡ sSFR 1/6 , given that L IR ∼ SFR 1.3 and M dust ∼ M * . At fixed L IR , dust temperature increases with redshift, as galaxies tend to have higher sSFR at higher redshifts. Σ SFR reflects the intensity of interstellar radiation on dust grains (∼ L/R 2 ), which sets the dust temperature by T ∼ Σ 1/6 SFR . The black dashed lines illustrate the 1/6-power scaling relations as argued above, which are in broad agreement with the simulations. The black arrows show how galaxies at the faint/bright end move if f dust increases by a factor of 2. The best-fit redshift-dependent T -L IR relation (Equation 9) for f dust = 0.2, 0.4, and 0.8 and for all dust temperature definitions are given in Table 2. In Fig. 9, we present the rest-frame dust SEDs for our z 5 simulations (using fdust = 0.4). In each panel, we collect all galaxies brighter than LIR = 10 10 L with peak wavelength falling in a 0.05 dex bin around the wavelength marked by the grey arrow (both λpeak and Tpeak are labeled in each panel). The median λpeak differs by 0.075 dex between adjacent panels. The shape of dust SED in mid-and far-IR does not strongly depend on LIR and redshift for galaxies in such narrow bins of λpeak, so we rescale all galaxies to LIR = 10 11.5 L . In each panel, we also show optically-thin MBB functions (Equation 5) at T = 35, 45, and 60 K (all normalized to 10 11.5 L ; red lines) for reference. The vertical cyan dashed lines in Fig. 9 label the observed-frame 1.2 mm (ALMA Band 6) at z = 6 (rest-frame 171 µm). For our z 5 galaxies, their flux densities at this wavelength are comparable to T ∼ 45-60 K MBB emission at the same LIR. For fdust = 0.8, the peak wavelengths increases and characteristic temperatures decrease by a factor of 2 1/6 = 1.12. The mid-IR SED (around λpeak) is shaped by warm dust, which is not accounted for by the single-temperature MBB function. Note that at rest-frame λ = 6-25 µm, the SED is usually dominated by PAH line emission (e.g. Baes et al. 2011), which is not included in the Weingartner & Draine (2001) SMC dust model. Our results at these wavelengths should be used with caution in this regard. Bouwens et al. (2016b) find that the ALMA Band 6 (observedframe 1.2 mm continuum) deep survey in the HUDF detects much fewer ∼ L * galaxies at z > 4 than what inferred from galaxy restframe UV slopes, even assuming the shallower SMC IRX-βUV re- lation, unless dust temperature increases to 44-50 K (as oppose to 35 K) at z > 4 (such that their far-IR fluxes are too weak to detect). They first derive the total IR luminosities of UV-selected galaxies from their UV fluxes and βUV and then convert LIR to observedframe 1.2 mm fluxes assuming the dust SEDs follow optically-thin MBB functions of assumed temperatures. Our simulated galaxies follow the SMC IRX-βUV relation, but a 35 K MBB function overestimates the flux density at ALMA Band 6 wavelength by a factor of 3-10. In other words, to convert between LIR and observed-frame 1.2 mm flux for z ∼ 6 galaxies using an optically-thin MBB function, one must assume a dust temperature of 45-60 K. These results support the hypothesis that the low detection rate of high-redshift galaxies in mm surveys is caused by galaxies falling below the detection limit because of their high dust temperatures (see also Faisst et al. 2017). The existence of an IRX-βUV relation close to the local or the SMC relation in z > 5 galaxies cannot be ruled out. In Fig. 10, we present the correlation between dust temperature and bolometric IR luminosity LIR (left), specific star forma-tion rate (over the past 10 Myr; middle), and average SFR surface density ( Σ SFR ≡ SFR10Myr/(Rmax/3) 2 , right) 13 . Each row shows one definition of dust temperature, with Tpeak in the top, Tmw in the middle, and Teff in the bottom. We only show a subsample of snapshots, color-coded by redshift. All the three temperatures correlate with LIR, consistent with the negative λpeak-LIR correlation shown in Fig. 8. At the same LIR, dust temperature increases with redshift. We fit the T -LIR relation for our simulated sample by the redshiftdependent power-law function We list the best-fit parameters for Tmw, Teff, and Tpeak and for fdust = 0.2, 0.4, and 0.8 in Table 2. Note that these fitting functions should only apply to star-forming galaxies below LIR ∼ 10 12 L at z = 5-12. Dust temperatures in z ∼ 2-4 DSFGs are studied in more detail in Liang et al. (2019) using a separate suite of FIRE simulations.
In contrast, neither the T -sSFR nor the T -Σ SFR relation depends on redshift. Here we want to provide simple, qualitative understanding on these correlations first, so we do not distinguish the three dust temperatures defined above for simplicity, although they are conceptually and physically different (see e.g. Liang et al. 2019 for more details). The T -sSFR relation can be understood, given LIR ∝ SFR 1.3 (Fig. 7) and Mdust ∝ M * (Fig. 2), as (see also Magnelli et al. 2014;Safarzadeh et al. 2016) where the second term is subdominant. As we have shown in Ma et al. (2018b), SFR increases with redshift at fixed stellar mass from z = 5-12 (see also Fig. 2). The redshift-dependence of the T -LIR relation can thus be attributed to the increasing sSFR with redshift (i.e. luminosity per unit dust mass; see also Imara et al. 2018). The T -Σ SFR relation is probably more physically expected. Given that SFR is proportional to the total luminosity from stellar sources and Rmax/3 is a characteristic scale of dust distribution, Σ SFR reflects the intensity of the interstellar radiation field on dust grains (i.e. Σ SFR ∼ L/R 2 ∼ J), which sets the dust temperature via energy balance. Following Equation 2, T should also scale to the 1/6 power of Σ SFR. The black dashed lines in Fig. 10 show the powerlaw scaling relations derived from the simple arguments above, in broad agreement with our simulated sample. Note that at a given luminosity, galaxies tend to be more compact at higher redshift (e.g. Oesch et al. 2010;Shibuya et al. 2015;Ma et al. 2018a), which also explains why dust temperature increases with redshift in the T -LIR relation. Similar to those in Fig. 8, the black arrows show how galaxies at the faint/bright end move if fdust increases by a factor of 2. Dust temperature does not change at the faint end, but decreases by a factor of ∼ 2 1/6 at the bright end.

The bright-end UV luminosity functions
In this section, we construct galaxy rest-frame UVLFs at z 5 using the entire simulation sample. Fig. 1 shows the number of halos that contain at least 10 4 particles and have zero contamination from low-resolution particles in all 34 zoom-in regions in every 0.25 dex bin from log Mhalo = 7.5-12 at selected redshifts. There are 57 snapshots from z = 12 to z = 5 with 15-20 Myr between snapshots. All halos above Mhalo = 10 10 M and central halos above 10 9.5 M are processed with SKIRT, for each of which we calculate mock images and SEDs along five lines of sight. We treat each halo snapshot and each sightline as independent 'galaxies'. We do not include subhalos and satellites in this work following Ma et al. (2018b). The rest-frame UV luminosity of high-redshift galaxies is usually measured in small apertures and only regions with sufficiently high surface brightness can be picked up (e.g. Ma et al. 2018a;Borlaff et al. 2018). To mimic these effects, we only include the light within an aperture of Rmax/3 in projected radius to exclude satellites and diffuse starlight (see e.g. Ma et al. 2018b, and Fig. 4 for examples). 14 For galaxies processed with SKIRT, we measure their UV luminosities directly from the mock image. For other galaxies, we project their star particles along a random sightline to produce an image, where the UV luminosity of each particle is calculated from the same stellar population synthesis models as in SKIRT for consistency. Note that they are all low-mass galaxies where dust attenuation is negligible (less than 0.01 mag seen in halos below Mhalo ∼ 10 10 M ).
At each redshift, we collect all 'galaxies' within a ∆z = ±0.5 interval. We count the number of objects in 36 halo mass bins from log Mhalo = 7.5-12 (i.e. bin width ∆ log Mhalo = 0.125 dex). On the other hand, we obtain the halo mass function (HMF) at this redshift using the public HMFcalc code , which agrees well with that directly extracted from our DM-only cosmological boxes. Every galaxy in the i th mass bin is assigned a weight representing its abundance in the universe, wi = φi∆ log M/Ni, where φi is the HMF evaluated at the bin center (in Mpc −3 dex −1 ), ∆ log M is the bin width (0.125 dex), and Ni is the number of galaxies in this bin. Next, all galaxies in the ∆z = ±0.5 redshift interval are divided in 30 equal-width bins of UV magnitude from MUV = −24-−14. 15 The number density of galaxies in each MUV bin is thus derived by summing over their weights. In Appendix B, we provide a detailed example about how we derive the z = 6 UVLF from our simulated sample for interested readers. We use a Schechter (1976) to fit the UVLFs derived from our simulations. We visually inspect the results to confirm that the best-fit Schechter function is always a good description for our simulated sample. In Fig. 11, we present the bright-end UVLFs from z = 5-10. Each panel shows the results at one redshift. The lines represent the best-fit Schechter functions for the UVLFs derived from our simulated sample, with the dashed lines showing the intrinsic UVLFs (without dust extinction) and the thin and thick dotted lines showing the post-extinction UVLFs for fdust = 0.4 and 0.8, respectively. Again, we remind that our experiments with different fdust here is equivalent to varying the dust opacity at fixed fdust. We compare our results with the most up-to-date observational constraints at these redshifts from wide-field deep surveys (e.g. McLure et al. 2013;14 Note that this is different from the UV luminosities in Section 3, where we include all the light within R vir . 15 The most luminous galaxies in our sample are slightly above L UV = 10 12 L , corresponding to M UV = −24. In this work, we are mainly interested in the bright-end UVLFs, so a lower limit at M UV = −14 is applied. We also assume that more massive halos contribute little to the UVLFs at these magnitudes, because of their low number densities in the universe as well as heavier dust obscuration in these systems. Oesch+14 Bouwens+16 Figure 11. The UVLFs from z = 5-10. Each panel represents one redshift. The lines show the best-fit Schechter functions for the UVLFs derived from our simulated sample, with the dashed line showing the intrinsic UVLFs without dust extinction and the thin and thick dotted lines show the post-extinction UVLFs for f dust = 0.4 and 0.8, respectively. We compare our results with the most up-to-date observational constraints at these redshifts from wide-field deep surveys. In all cases, the faint-end slope α steepens, the break magnitude M * UV increases, and the normalization φ * UV decreases with redshift, which are inherited from the redshift evolution of the HMFs. The faint-end UVLFs are not strongly affected by f dust . The bright-end (M UV < −21) UVLFs are mainly determined by dust attenuation, with M * UV increasing with f dust . At z 7, the UVLFs where f dust = 0.8 agree better with observation than f dust = 0.4, although there is still a small (less than a factor of 2) discrepancy at the bright end. Such discrepancy disappears at z 8, tentatively suggesting that dust properties in z = 9-10 galaxies are different than those in z = 5-6 galaxies. Table 3. Best-fit Schechter function (Equation 11) for the UVLFs derived from our simulated sample shown in Fig. 11. The parameters are (α, M * UV , φ * UV ).  Oesch et al. 2013Oesch et al. , 2014Bouwens et al. 2015Bouwens et al. , 2016aFinkelstein et al. 2015;Laporte et al. 2015;Bowler et al. 2017;Stefanon et al. 2017;Ono et al. 2018, symbols with errorbars). We also provide the best-fit parameters of the Schechter functions for our UVLFs in Table 3 for reference.
In all three cases ( fdust = 0, 0.4, and 0.8), the faint-end slope α becomes steeper, the break magnitude M * UV increases (i.e. becomes fainter), and the normalization φ * UV decreases with increasing redshift. These features are primarily inherited from the redshift evolution of the HMFs. The UVLFs at MUV > −19 are not strongly affected by dust attenuation and the faint-end slope α remains unchanged with fdust at a given redshift. On the other hand, the brightend UVLFs are determined by dust extinction, with the break magnitude M * UV increasing with fdust, consistent with previous results found by different authors (e.g. Cullen et al. 2017;Wilkins et al. 2017;Ma et al. 2018b;Yung et al. 2018). The intrinsic UVLFs are above the observational constraints at the bright end and dust attenuation reduces the number of bright galaxies at any redshift.
At UV magnitude MUV > −19, the UVLFs derived from our simulations agree well with observations regardless of fdust as dust attenuation is always subdominant in this regime. For fdust = 0.4, the bright-end UVLFs still lie above the observational constraints at z 8. The fdust = 0.8 UVLFs agree better with observations, but there is still a small discrepancy (within a factor of 2) below z = 7 at MUV < −21. Interestingly, such discrepancy disappears at z 8 for fdust = 0.8 and even the UVLFs for fdust = 0.4 agree well with observations at z = 9 and 10. Although we note that the UVLFs at z > 8 are still poorly constrained, our results here show a tentative evidence that dust properties in z = 9-10 galaxies may be different than those in z = 5-6 galaxies (e.g. dust fraction by mass is possibly lower at z 9). This is not unreasonable because the cosmic time at z 9 is too short for dust production from asymp- totic giant branch stars in contrast to the local Universe (e.g. Dwek et al. 2014). Nonetheless, we suggest that better constraints of the bright-end UVLFs at z > 8 with ongoing and future wide-field deep surveys can improve our understanding on dust formation and dust properties in the very early Universe in the foreseeable future.

The IR luminosity functions
We also construct the bolometric IRLFs at z = 5-10. Unlike restframe UV, the IR emission is nearly isotropic, so we do not account for line-of-sight variations, but only include each galaxy snapshot once in our analysis below. Again, at each redshift, we collect all galaxy snapshots within ∆z = ±0.5 and assign weights to halos in 36 equal-with mass bins from log Mhalo = 7.5-12 as in Section 4.1. We divide all galaxies in 15 equal-width bins from log LIR = 7-12 and obtain the number density of galaxies in each LIR bin by adding their weights. Note that most galaxies brighter than LIR ∼ 10 7 L have been processed with SKIRT. For those without dust radiative transfer calculations, their LIR are derived from LUV using Equation 3. This has little effect on our results.
In Fig. 12, we show the derived IRLFs at z = 6, 8, and 10 (symbols) for fdust = 0.4 (left) and 0.8 (right). We fit our results at integer redshifts from z = 5-10 all together using a redshift dependent double power-law function (cf. Casey et al. 2018a) where α1, α2, L * IR , and φ * IR are power-law functions of 1 + z. We show the best-fit IRLFs at z = 6, 8, and 10 in Fig. 12  for fdust = 0.8, respectively. The redshift-dependent double powerlaw function describes our results very well. Note that our simulation sample only covers up to LIR ∼ 10 12 L and does not capture rare, most heavily obscured, extremely luminous IR galaxies (e.g. Our results broadly agree with observations. Using f dust = 0.4 underestimates the obscured fraction at z < 8, but a heavy dust attenuation is not required at higher redshifts. LIR ∼ 10 13 L ), so our results should not be extrapolated to higher LIR without caution.

The cosmic star formation rate density
Current observational constraints on the cosmic SFRD at z 5 are converted from rest-frame UV luminosity density using L1500 = 8.0 × 10 27 SFR M yr −1 erg s −1 Hz −1 .
In this section, we calculate the dust (un)obscured UV luminosity densities at z = 5-10 using our simulations. At each redshift, we derive the UVLF following the steps in Section 4.1 and integrate the best-fit Schechter function over MUV < −17 (as most observational studies do) to compute the UV luminosity density at that redshift. In Fig. 13, we present our results for fdust = 0 (unobscured, dashed line), 0.4 (thin dotted line), and 0.8 (thick dotted line) and compare with observational constraints (symbols with errorbars; e.g. Ellis et al. 2013;Oesch et al. 2013Oesch et al. , 2014Bouwens et al. 2015;Finkelstein et al. 2015;McLeod et al. 2016;and CLASH detections from Zheng et al. 2012;Coe et al. 2013;Bouwens et al. 2014). The open (filled) symbols show the (un)obscured results, respectively. Our predicted unobscured UV luminosity density agrees fairly well with current observational constraints. Similar to the case with UVLFs in Section 4.1, we find our fdust = 0.8 results agree better with the observed dust obscured UV luminosity density than those using fdust = 0.4, as the latter underestimate the obscured fraction of the UV light at z < 8. At higher redshifts, the difference between fdust = 0.4 and 0.8 becomes much smaller and thus a heavy dust attenuation is no longer required at z 8.

DISCUSSION
5.1 Strategies for probing dusty galaxies at z 5 In this paper, we present a broad spectrum of predictions on dust extinction and emission in high-redshift galaxies that can be tested . The grey dashed line shows the consensus relation for z ∼ 2-3 galaxies (Reddy et al. 2010;Whitaker et al. 2014;Álvarez-Márquez et al. 2016). The thick segment represents the mass range where current observational constraints are available, while the thin segment shows the interpolation of the relation to lower masses. For individual galaxies, IRX is typically smaller by ∼ 0.3 dex for f dust = 0.4 than f dust = 0.8 (each galaxy appears as a pair of red and black points at the same M * ). This suggests that the IRX-stellar mass relation can be used to constrain f dust . and motivate future observations. First of all, we argue that current data cannot completely rule out the existence of an IRX-βUV relation in z 5 galaxies that is consistent with local relations. The UVcontinuum slope of a galaxy may still be a good indicator of dust attenuation. Where galaxies lie on the IRX-βUV relation is predominantly determined by the shape of the dust extinction curve in the UV, which reflects the dust composition. Better constraints of the IRX-βUV relation at z 5 with the possibility of constraining the attenuation law in the rest-frame optical in the future can help understand dust formation history and the evolution of dust properties across cosmic time.
We predict that dust temperatures can be much higher in z 5 galaxies than in low-redshift galaxies, a consequence of high sSFR (luminosity per unit dust mass) and/or high SFR surface densities (intensity of radiation on dust grains) in high-redshift galaxies. This can be tested using multi-band observations of dust emission from a sample of intrinsically bright, dust obscured galaxies at z 5. It would be interesting to select targets based on their UV-continuum slopes, but the observations must be deeper than current surveys in the (sub)mm, as the flux densities at the R-J tail are reduced by a factor of a few (cf. Fig. 9). Having coverage on at least one band at wavelength shorter than rest-frame λpeak is critical for measuring dust temperature and bolometric dust luminosity, which will be achievable with the OST (wavelength coverage from 5-600 µm). This also helps us understand whether galaxies with red UV slopes but low apparent IRX are real or just because their bolometric IR luminosities are underestimated from single-wavelength data due to the presence of warmer dust.
We find that the IRX-βUV relation does not depend on dust fraction or dust-to-gas ratio. On the other hand, as shown in Section 4.1, the shape of the bright-end UVLFs is sensitive to dust fraction. We therefore propose that the bright-end UVLFs can be combined with IRLFs, the IRX-βUV relation, and other observables at long wavelengths as a new method to infer dust properties in z 5 galaxies in a statistical sense. Ongoing and future observations with the Hubble Space Telescope and the James Webb Space Telescope in the rest-frame UV as well as current and next-generation radio telescopes (e.g. ALMA, CSST, ngVLA, OST) probing dust emission in high-redshift galaxies are very promising to this end.
In Fig. 14, we also examine the IRX-stellar mass relation for our simulated sample, using fdust = 0.4 (black points) and 0.8 (red points), respectively. Each galaxy appears as a pair of red and black points at the same M * . This relation is independent of redshift and we only show one sightline for each galaxy snapshot at integer redshift from z = 5-12. The IRX-stellar mass relation depends on fdust, with IRX decreasing roughly by ∼ 0.3 dex for individual galaxies if fdust drops from 0.8 to 0.4. The grey dashed line shows the consensus z ∼ 2-3 IRX-stellar mass relation (e.g. Reddy et al. 2010;Whitaker et al. 2014;Álvarez-Márquez et al. 2016). The thicker part represents relatively massive galaxies (i.e. M * > 10 9 M ) for which current observational constraints are available at z ∼ 2-3, while the thinner part represents its interpolation to lower masses. Our simulations broadly agree with this relation, in line with the results in Bouwens et al. (2016b) where they find that the inferred IRX-stellar mass relation of typical ∼ L * galaxies at z ∼ 4-10 from 1.2 mm ALMA-HUDF deep survey is consistent with the z ∼ 2-3 relation if dust temperature increases with redshift to ∼ 44-50 K at z ∼ 6. Our results suggest that better constraints on the IRX-stellar mass relation at z 5 can also be used to infer dust fraction in z 5 galaxies in addition to bright-end UVLFs.

The effects of the CMB
In this work, we include starlight as the sole source heating the dust and only study the 'intrinsic' dust emission and dust temperature. We note that heating from the CMB can play a significant role at z 5, when the CMB temperature starts to become comparable to the dust temperature. Following the argument in da Cunha et al. (2013) for single-temperature dust, the CMB first heats the dust to a higher temperature where β = 2 is the dust emissivity index and TCMB is the CMB temperature at the redshift of interest. Second, the CMB serves as a background which the dust emission from high-redshift galaxies is measured against. Subtracting this background reduces the observed flux by a factor of 1 − Bν (TCMB)/Bν (Tdust,withCMB) at a given frequency. In general, galaxies with higher intrinsic dust temperatures are less affected than those with primarily cold dust. The net effect is stronger at longer wavelengths than at shorter wavelengths. We exclude the CMB in this paper on purpose for two reasons. First, empirical models of number counts or luminosity functions of high-redshift galaxies in the IR and (sub-)mm often start from intrinsic dust emission and then convert to observed flux following da Cunha et al. (2013) as summarized above (e.g. Bouwens et al. 2016b;Casey et al. 2018a). Therefore, it is important to understand the dust SEDs and dust temperatures in z 5 galaxies without the CMB. Second, the effects of the CMB are more complicated in reality given the broad distribution of dust temperature in individual galaxies. Warm dust close to young stars is barely affected, while the diffuse dust at much lower temperature mostly becomes invisible. It is not clear which dust temperature is applicable to Equation  14 and by what fraction the observed flux is reduced in different regions of a galaxy. In a future study, we will investigate how CMB heating affects the observed far-IR flux from our simulated galaxies using full radiative transfer calculations where we include the CMB as an extra source that produces a uniform radiation field with a black-body spectrum.

Limitations of this work
Our dust radiative transfer calculations assume a fixed dust composition and dust-to-metal ratio everywhere in a galaxy as well as in all galaxies. In reality, dust composition may vary in different regions of a galaxy, as seen in the MW where the extinction curve varies between lines of sight. Moreover, the dust-to-metal ratio in cold, dense gas is presumably higher than that in warm, diffuse gas, because dust growth is more efficient and the grains are less likely to be destroyed in cold, dense gas. This may further enlarge the discrepancy of dust attenuation between young stars just born in dense clouds and relatively older stars preferentially living in more diffuse gas, leading to a dramatic effect on the galaxy-averaged extinction law. Furthermore, in the local Universe, it has been suggested that the dust-to-metal ratio decreases at low metallicity (below ∼ 0.1 Z , e.g. Rémy-Ruyer et al. 2014). However, it is not clear if it is also the case at high redshifts, given that the ISM conditions are very different from those at lower redshifts. Given that the dust properties in high-redshift galaxies are still poorly constrained because there are only a small set of data available so far, we adopt the simplest treatments in this work to qualitatively predict what to expect for future observations. More sophisticated treatments of dust physics will be necessary in the future if new data suggest so.
Last but not least, our simulation sample only includes normal star-forming galaxies that are typically discovered in current deep surveys in the rest-frame UV. The most massive galaxies in our sample have stellar mass ∼ 10 10.5 M and bolometric IR luminosity ∼ 10 12 L . We do not yet simulate more massive, heavily obscured, and luminous systems (e.g. M * ∼ 10 11 M , LIR ∼ 10 13 L ) at these redshifts at comparably high resolution, like the extremely luminous DSFGs detected by SPT (e.g. Strandet et al. 2016). They are relatively rare objects that may involve major mergers of two massive galaxies or rapidly accreting supermassive black holes. It is not clear where such galaxies lie on the IRX-βUV (stellar mass) relation and whether they have higher dust temperatures than lowerredshift galaxies at similar luminosities. There is no guarantee that our predictions in this paper still hold for more massive and luminous systems.

CONCLUSIONS
In this work, we utilize a suite of 34 cosmological zoom-in simulations that consist of thousands of sufficiently resolved halos spanning a halo mass range Mhalo ∼ 10 8 -10 12 M with stellar mass up to ∼ 10 10.5 M and intrinsic UV luminosity up to ∼ 10 12 L (MUV ∼ −24) at z 5. These simulations use the FIRE-2 models of the multi-phase ISM, star formation, and stellar feedback. With a mass relation of 7000 M or better and typical spatial resolution in dense gas of 1 pc, these simulations explicitly resolve star formation in dense birth clouds and feedback destroying these clouds. We post-processing all halos above Mhalo = 10 10 M and central halos above 10 9.5 M in our sample using the three-dimensional Monte Carlo dust radiative transfer code SKIRT to study dust extinction, dust emission, and dust temperature in high-redshift galaxies. Our calculations assume a SMC-like dust composition from Weingartner & Draine (2001) and a constant dust-to-metal ratio fdust in all gas below 10 6 K (no dust in hotter gas). We fix dust composition and opacity but experiment with different fdust, which accounts for uncertainties of dust opacity and dust-to-metal ratio in a single parameter. We do not adopt any models for sub-resolution dust distribution but instead process the simulations directly. Our main findings include the following.
(i) Dust geometry is clumpy and patchy. The young stars emitting most of the UV photons are usually concentrated in the central region of the galaxy, but dust is distributed on much larger spatial scales. Dust scatters UV light to an extended distribution at relatively low surface brightness, which contributes a non-negligible fraction of the escaped UV flux (Fig. 4).
(ii) Our sample shows a tight relationship between IR excess (IRX) and UV-continuum slope (βUV), consistent with the SMC IRX-βUV relation, despite the patchy dust geometry in our simulations. Galaxies at higher redshifts tend to move slightly to bluer βUV at fixed IRX due to their younger stellar population. Viewing the same galaxy from different sightlines gives the same IRX-βUV relation as the entire sample (Fig. 5, top panel).
(iii) The IRX-βUV relation does not depend on the normalization of the attenuation law (represented by different fdust). However, it does depend on the shape of the extinction curve (Fig. 5, bottom panel), which reflects the dust composition.
(iv) Our simulations produce an IRX-stellar mass relation in broad agreement with the consensus relation established for z ∼ 2-3 galaxies. This relation depends on fdust, with IRX decreasing by ∼ 0.3 dex if fdust drops from 0.8 to 0.4 (Fig. 14).
(v) There is a positive correlation between bolometric IR luminosity LIR and intrinsic UV luminosity LUV,intr, which can be described by a broken power-law function (Equation 3). At the bright end, LIR changes little with fdust (the optically-thick limit), while at the faint end, LIR is proportional to fdust (the optically-thin limit). The LIR-LUV,intr relation does not depend on redshift (Fig. 6, left  panel).
(vi) The scatter in the LIR-LUV,intr relation is not driven by dust mass, average dust column density, nor dust density, although all three quantities are redshift-dependent. This suggests that dust luminosity is mainly determined by dust covering fraction. There is a secondary correlation between LIR and the SFR averaged over the past 10 Myr (i.e. the amount of stars younger than 10 Myr) at a given LUV,intr, because young stars are more heavily obscured than relatively older stars (Fig. 6, right panel).
(vii) The correlation between LIR and SFR10Myr for the entire sample can be well described by a power-law function. The scatter of this relation is driven by LUV,intr. This further confirms the differential obscuration between stars younger and older than 10 Myr (Fig. 7). Note that LUV,intr and SFR10Myr are not fully degenerated.
(viii) Our simulated sample shows an anti-correlation between the peak wavelength λpeak of dust emission (in terms of Lν ) and LIR. However, the λpeak-LIR relation for z 5 galaxies shows a large offset from the observed relation at lower redshifts, with λpeak moving toward shorter wavelengths by a factor of 2 at a given LIR (Fig. 8), suggesting that dust is on average warmer in high-redshift galaxies.
(ix) The dust SEDs are far from an optically-thin MBB function. At z = 6, the flux densities at ALMA Band 6 (observed-frame 1.2 mm) of our simulated galaxies are comparable to MBB spectra with T ∼ 45-60 K at the same LIR (Fig. 9). The low detection rate of dust continuum at z 5 compared to what inferred from the UV slopes is likely due to higher dust temperatures in these galaxies.
(x) We predict that dust temperature correlates positively with both sSFR (approximately dust luminosity per unit mass) and SFR surface density (intensity of the interstellar radiation). Both corre-lations are independent of redshift. At fixed LIR, dust temperature increases with redshift from z = 5-12 (Fig. 10), because galaxies at higher redshifts tend to have higher sSFR and more compact SFR. Dust temperature does not change significantly with fdust at the faint end, but increases by a factor of 2 1/6 = 1.12 if fdust increases by a factor of 2.
(xi) Using the entire simulation sample, we derive the UVLFs from z = 5-10. The bright-end UVLFs are largely determined by dust attenuation. By comparing our results with most up-to-date observational constraints, we find tentative evidence that dust properties are likely evolving from z = 10 to z = 5 (Fig. 11, Table 3). We suggest that better measurements of the bright-end UVLFs at z > 8 with future observations provide a powerful probe of dust physics in the very early Universe.
(xiii) We derive dust (un)obscured UV luminosity density and cosmic SFRD at z = 5-10 from the UVLFs. Our results are broadly consistent with observational constraints in the literature (Fig. 13).  (2001). The red dashed line uses the same parameters as the black line, except that dust self-absorption is ignored. The SEDs using MW-type dust and SMC-type dust models mainly differ (1) in the mid-IR (rest-frame 3-25 µm) due to PAH emission and (2) the near-UV absorption feature because of the 2175 Å bump in the MW-like extinction curve. The symbols show that our default choices of dust grid resolution and number of photon packets ensure excellent resolution.
The mid-IR SED (6-25 µm) differs significantly between MW-type dust and SMC-type dust due to PAH emission. For MW-type dust model, there is a strong absorption feature in the near-UV caused by the 2175 Å bump in the extinction curve, making the rest-frame UV slope βUV always negative, so we do not use MW-type dust as our default choice. Nevertheless, at rest-frame λ > 30 µm, both models give nearly identical results. We note that dust self-absorption is important around the peak wavelength of dust emission.
The symbols show resolution tests using MW-type dust model without dust self-absorption on a 25-point wavelength grid. Orange circles show the results using 10 7 photon packets per wavelength. Blue triangles show the results using a factor of 4 worse resolution for the dust grid (i.e. each cell has a maximum dust mass equals to 4 × 10 −6 of the total dust mass in the domain). Green square show the results using a factor of 5 better resolution for the dust grid (the maximum dust mass in each cell is 2 × 10 −7 of the total dust mass). These calculations agree precisely well with each other, suggesting that our default choices for grid resolution and the number of photon packets (the red dashed line) ensure excellent convergence.
In Fig. A2, we show the same LIR-LUV,intr relation as in the left panel of Fig. 6, but separate galaxies simulated at mass resolution m b ∼ 7000 M and at resolution m b ∼ 900 M or better with grey and red points, respectively. The two subsamples of our simulated galaxies overlap at halo mass Mhalo 10 10.5 M . There is no significant difference on the LIR-LUV,intr relation between the two subsamples. We also explicitly check all the results in this paper and find that none of our conclusions is sensitive to the resolution of our simulations, at least in the mass range where the two subsamples overlap.

APPENDIX B: HOW WE DERIVE UVLFS FROM THE SIMULATED SAMPLE
In Section 4.1, we briefly describe the methods to construct UVLFs using our simulated sample. Here we walk through the steps in detail for interested readers using an example, the z = 6 UVLF. The histogram in the left panel of Fig. B1 shows the number of 'galaxies' in 36 equal-width bins of log Mhalo from 7.5-12 (i.e. bin width 0.125 dex) in z = 5.5-6.5. Note that we treat those in different snapshots and the same galaxy viewed along different sightlines as different objects. The black solid line in the left panel shows the HMF at z = 6 obtained from HMFcalc code. Every galaxy in the same halo mass bin is given the same weight representing its number density in the universe, w = φ ∆ log M/N, where φ is the HMF evaluated at the bin center, ∆ log M = 0.125 dex is the bin width, and N is the number of objects in this bin.
In the right panel, we show the number of galaxies in 30 bins of MUV from −24 mag to −14 mag (orange cross). The blue circles show the total weight of galaxies in each MUV bin, divided by the bin width (i.e. 1/3 mag). This is thus the UVLF derived from our simulated sample. Note that there is small noise in the result, which is caused by the fluctuations in the number of galaxies in each MUV bin. Using a smaller number of bins will reduce the noise, but do not affect our results significantly. The black dashed line shows the best-fit Schechter function of the z = 6 UVLF (Equation 11). We visually inspect and confirm that a Schechter function is always a good description of the UVLF derived from our sample at any redshift. Number of galaxy snapshots # of galaxies in each bin Figure B1. Left: Number of galaxies in 36 bins of halo mass for all galaxies in 5.5 < z < 6.5 in our sample from log M halo = 7.5-12 (histogram). The black solid line shows the z = 6 HMF calculated from the HMFcalc code. Right: Number of galaxies in 30 bins of M UV from −24 to −14 mag (orange symbols). The blue circles show the z = 6 UVLF derived from our simulated sample, by summing the weight over all galaxies in each M UV bin and dividing it by the bin width. Using a smaller number of M UV bins reduces the noise but does not affect our results. The black dashed line shows the best-fit Schechter function of the z = 6 UVLF (Equation 11).