On the dust temperatures of high redshift galaxies
Lichen Liang Robert Feldmann Dušan Kereš Nick Z. Scoville Christopher C. Hayward Claude-André Faucher-Giguère Corentin Schreiber Xiangcheng Ma Philip F. Hopkins Eliot Quataert


firstpage–lastpage
Tachyonic  δ-Tsallis entropy of  a thermal tachyonic BIon
     Aroonkumar Beesham^1[beeshama@unizulu.ac.za], Alireza Sepehri^1,2[alireza.sepehri3@gmail.com]
    December 30, 2023
==================================================================================================




Dust temperature is an important property of the interstellar medium (ISM) of galaxies. It is required when converting (sub)millimeter broadband flux to total infrared luminosity (L_ IR), and hence star formation rate, in high-redshift galaxies. However, different definitions of dust temperatures have been used in the literature, leading to different physical interpretations of how ISM conditions change with, , redshift and star formation rate. In this paper, we analyse the dust temperatures of massive (M_ star>10^10 M_⊙) z=2-6 galaxies with the help of high-resolution cosmological simulations from the Feedback in Realistic Environments (FIRE) project. At z∼2,  our  simulations  successfully  predict  dust temperatures in good agreement with observations. We find that dust temperatures based on the peak emission wavelength increase  with redshift, in line with the higher star formation activity at higher redshift, and are strongly correlated with the specific star formation rate. In contrast, the mass-weighted dust temperature, which is required to accurately estimate the total dust mass, does not strongly evolve with redshift over z=2-6 at fixed IR luminosity but is tightly correlated with L_ IR at fixed z. We also analyse an `equivalent' dust temperature for converting (sub)millimeter flux density to total IR luminosity, and provide a fitting formula as a function of redshift and dust-to-metal ratio. We find that galaxies of higher equivalent (or higher peak) dust temperature (`warmer dust') do not necessarily have higher mass-weighted temperatures. A `two-phase' picture for interstellar dust can explain the different scaling relations of the various dust temperatures.



galaxies: evolution — galaxies: high-redshift — galaxies: ISM — submillimetre: galaxies




§ INTRODUCTION



Astrophysical dust, originating from the condensation of metals in stellar ejecta, is pervasive in the interstellar medium (ISM) of galaxies in both local and distant Universe <cit.>. Dust scatters and absorbs UV-to-optical light, and therefore strongly impacts the observed flux densities and the detectability of galaxies at these wavelengths <cit.>. Despite that it accounts for no more than a few percent of the total ISM mass <cit.>, dust also plays a key role in star formation process of galaxies <cit.>. Constraining and understanding dust properties of galaxies is therefore essential for proper interpretation of the multi-wavelength data from observations and for facilitating our understanding of galaxy formation and evolution.

Much of the stellar emission of star-forming galaxies is absorbed by dust grains and re-emitted at infrared (IR)-to-millimeter (mm) wavelengths as thermal radiation, encoding important information about dust and galactic properties, such as dust mass, total IR luminosity[In this paper, L_ IR is defined as the luminosity density integrated over the 8-1000 μ m wavelength interval.] (L_ IR) and star formation rate (SFR) <cit.>. The advent of the new facilities in the past two decades, such as the Spitzer Space Telescope <cit.>, Herschel Space Observatory <cit.>, the Submillimetre Common-user Bolometer Array (SCUBA) camera on the James Clerk Maxwell Telescope (JCMT) <cit.>, the AzTEC millimeter camera on the Large Millimeter Telescope (LMT) <cit.>, the South Pole Telescope (SPT) <cit.> and the Atacama Large Millimeter/sub-millimeter Array (ALMA) has triggered significant interests in the study of ISM dust. In particular, observations with the Photodetector Array Camera and Spectrometer <cit.> and the Spectral and Photometric Imaging Receiver <cit.> instruments aboard Herschel made it possible to study the 70-500 μ m wavelength range where most of the Universe's obscured radiation emerges, and many dust-enshrouded, previously unreported objects at distant space have been uncovered through the wide-area extra-galactic surveys <cit.>. Far-infrared (FIR)-to-mm spectral energy distribution (SED) modelling of dust emission has therefore become possible for objects at high redshift <cit.> and various dust properties can be extracted using SED fitting techniques <cit.>. 

The Rayleigh-Jeans (RJ) side ( λ>λ_ peak, where λ_ peak is the wavelength of the peak emission) of the dust SED can empirically be well fit by a single-temperature (T) modified blackbody (MBB) function <cit.>. However, the shape of the Wien side ( λ<λ_ peak) of the SED, which is tied to the warm dust component in vicinity of the young stars and active galactic nuclei (AGN), has a much larger variety <cit.>. Studies have shown that one single-T MBB function cannot well fit both sides of the observed SEDs <cit.>. Motivated by this fact, fitting functions of multi-T components have been adopted <cit.>. Meanwhile, empirical SED templates have been developed based on an assumed distribution of interstellar radiation intensity (U) incident on dust grains <cit.>. Both approaches can produce functional shape of dust SED that better matches the observed photometry of galaxies compared to a single-T MBB function.

At high redshift (z2), however, it is more common that a galaxy has only a few (two or three) reliable photometric data points in its dust continuum so that SED fitting by multi-T functions or more sophisticated SED templates is not possible. Therefore, the widely adopted approach is to simply fit the available data points with one single-T MBB function <cit.>. The T parameter that yields the best-fit is then often referred to as the `dust temperature' of the galaxy in the literature. We specify this definition of dust temperature as the `effective' temperature (T_ eff) in this paper. Another temperature also often used is the `peak' temperature (T_ peak), which is defined based on the emission peak of the best-fit SED assuming Wien's displacement law <cit.>. These observationally-derived dust temperatures (both T_ eff and T_ peak) can depend on the assumed functional form of SED as well as the adopted photometry <cit.>. Despite that it is unclear how well these simplified fitting functions represent the true SED shape of high-redshift galaxies and the physical interpretation of the derived temperatures is not obvious, this approach is frequently used to analyse large statistical samples of data <cit.>.

The scaling relations of T_ peak (T_ eff) and other dust/galaxy properties, including the L_ IR-temperature and specific star formation rate (sSFR)-temperature relations, may be related to the physical conditions of the star-forming regions in distant galaxies and have attracted much attention <cit.>. While observational studies derive dust temperatures in a variety of ways (using different fitting functions and/or different photometry), they generally infer that the temperature increases with L_ IR and sSFR of galaxies. Accurate interpretation of these observed scaling relationships requires a knowledge of how different dust and galaxy properties ( stellar mass, SFR and dust mass) shape the dust SED <cit.>, and hence the derived dust temperatures. Radiative transfer (RT) analyses of galaxy models are important tools for understanding these temperatures since the intrinsic properties of the simulated galaxies are known <cit.>.

One important question is how the derived temperatures are related to the physical, mass-weighted dust temperature (T_ mw). Observations of the local galaxies have shown that the bulk of the ISM dust remains at low temperature <cit.>. The cold dust component determines T_ mw of the galaxy, which sets the shape of the RJ tail. For distant galaxies, it is very challenging to measure T_ mw due to the limit of resolution. However, a good estimate of T_ mw is important for deriving the ISM masses of high-redshift galaxies via the RJ method <cit.>. It is unclear whether, or how, one can infer T_ mw from the observationally-derived temperatures. Alternatively, one can simply adopt a constant value if T_ mw has relatively small variation among different galaxies, given that the mass estimates scale only linearly with T_ mw <cit.>. If that is the case, it can also be one major advantage of the RJ approach because the main difficulty of the traditional CO method is the large uncertainty of the CO-to-H_2 conversion factor <cit.>. RT analyses are useful for understanding the relation between the derived temperatures and T_ mw <cit.>.

Over the past two decades, many ground-based galaxy surveys at (sub)mm wavelengths ( SCUBA, AzTEC, SPT and ALMA) that are complementary to Herschel observations <cit.> have been carried out. Deep (sub)mm surveys are capable of probing less actively star-forming (SFRs100 M_⊙ yr^-1) galaxies at z4  <cit.>. Furthermore, they are effective at uncovering sources at z>4 thanks to the effect of “negative-K correction" <cit.>.
The (sub)mm-detected sources do not necessarily have Herschel counterparts. Deriving L_ IR (and hence SFR) of these sources from a single (sub)mm flux density (S) requires adopting an assumed dust temperature, which we refer to as `equivalent' temperature (T_ eqv) in this paper, along with an assumed (simplified) functional shape of the dust SED <cit.>. T_ eqv is conceptually different from T_ eff introduced above because the former is an assumed quantity for extrapolating L_ IR from a single data point while the latter is a derived quantity through SED fitting to multiple data points.

A good estimate of T_ eqv is important for translating the information ( source number counts) extracted from the ALMA blind surveys to the obscured cosmic star formation density at z4 <cit.>, where currently only reliable data from rest-frame UV measurements are available <cit.>. One common finding of the recent (sub)mm blind surveys is a dearth of faint sources at these early epochs, as noted by <cit.> (see also the references therein). This can suggest that the early Universe is relatively dust-poor and only a small fraction of stellar emission is absorbed and re-emitted by dust <cit.>. Alternatively, it could also be accounted for by a significantly higher T_ eqv at high redshifts, meaning that galaxies of the same L_ IR appear to be fainter in the (sub)mm bands <cit.>.  Hence, understanding how T_ eqv evolves with redshift and how it depends on different galaxy properties are crucial for constraining the total amount of dust and the amount of obscured star formation density in the early Universe <cit.>.

In this paper, we study in detail the observational and the physical (mass-weighted) dust temperatures with the aid of high-resolution cosmological galaxy simulations. In particular, we study a sample of massive (M_ star>10^10 M_⊙) z=2-6 galaxies from the FIRE project[<fire.northwestern.edu>] <cit.> with dust RT modelling. This sample contains galaxies with L_ IR ranging over two orders of magnitude, from 10^10 to 10^12 L_⊙ and few dust-rich, ultra-luminous (L_ IR10^12 L_⊙) galaxies at z∼2 that are candidates for both Herschel- and submm-detected objects. A lot of them have L_ IR∼ a few ×10^11 L_⊙, which is accessible by Herschel using stacking techniques <cit.>. Our sample also contains fainter galaxies at z=2-6 with observed flux densities S_ 870μ m (S_ 1.2 mm)0.1 mJy, which could be potentially detected with ALMA. We calculate and explicitly compare their T_ mw with the observationally-derived temperatures (T_ peak or T_ eff), as well as their scaling relationships with several galaxy properties. We also provide the prediction for T_ eqv that is needed for deriving L_ IR of galaxy from its observed single-band (sub)mm flux.

The paper is structured as follows. In Section <ref>, we introduce the simulation details and the methodology of radiative transfer modelling. In Section <ref>, we provide the various definitions of dust temperature in detail, discuss the impact of dust-temperature on SED shape, and compare the specific predictions of our simulations with observations. In Section <ref>, we focus on the conversion from single-band (sub)mm broadband flux to L_ IR and provide useful fitting formulae. In Section <ref>, we discuss the observational implications of our findings. We summarise and conclude in Section <ref>.

Throughout this paper, we adopt cosmological parameters in agreement with the nine-year data from the Wilkinson Microwave Anisotropy Probe <cit.>, specifically Ω_ m=0.2821, Ω_Λ=0.7179, and H_0=69.7  km s^-1 Mpc^-1.



§ SIMULATION METHODOLOGY



In this section, we introduce our simulation methodology. In Section <ref>, we briefly summarize the details of the cosmological hydrodynamic simulations from which our galaxy sample is extracted. In Section <ref>, we introduce the methodology of our dust RT analysis and present mock images produced with SKIRT.
  


 §.§ Simulation suite and sample



We extract our galaxy sample from the MassiveFIRE cosmological zoom-in suite <cit.>, which is part of the Feedback in Realistic Environments (FIRE) project. 

The initial conditions for the MassiveFIRE suites are generated using the MUSIC code <cit.> within a (100  Mpc/h)^3 comoving periodic box with the WMAP cosmology. From a low-resolution (LR) dark matter (DM)-only run, isolated halos of a variety of halo masses, accretion history and environmental over-densities are selected. 
Initial conditions for the `zoom-in' runs use a convex hull surrounding all particles within 3R_ vir at z=2 of the chosen halo defining the Lagrangian high-resolution (HR) region. The mass resolution of the default HR runs are m_ DM=1.7×10^5 M_⊙ and m_ gas=3.3×10^4 M_⊙, respectively. The initial mass of the star particle is set to be the same as the parent gas particle from which it is spawned in the simulations. 

The simulations are run with the gravity-hydrodynamics code GIZMO[A public version of GIZMO is available at <http://www.tapir.caltech.edu/phopkins/Site/GIZMO.html>] (FIRE-1 version) in the Pressure-energy Smoothed Particle Hydrodynamics (“P-SPH") mode <cit.>, which improves the treatment of fluid mixing instabilities and includes various other improvements to the artificial viscosity, artificial conductivity, higher-order kernels, and time-stepping algorithm designed to reduce the most significant known discrepancies between SPH and grid methods <cit.>. Gas that is locally self-gravitating and has density over 5  cm^-3 is assigned an SFR ρ̇=f_ molρ/t_ ff, where f_ mol is the self-shielding molecular mass fraction. The simulations explicitly incorporate several different stellar feedback channels (but not feedback from supermassive black holes) including 1) local and long-range momentum flux from radiative pressure, 2) energy, momentum, mass and metal injection from supernovae (Types Ia and II), 3) and stellar mass loss (both OB and AGB stars) and 4) photo-ionization and photo-electric heating processes. We refer the reader to <cit.> for details.

In the present study we analyse 18 massive (10^10<M_ star<10^11.3 M_⊙ at z=2) central galaxies (from Series A, B and C in ) and their most massive progenitors (MMP) up to z=6, identified using the Amiga Halo Finder <cit.>. These galaxies were extracted from the halos selected from the LR DM-only run. In order to better probe the dusty, IR-luminous galaxies at the extremely high-redshift (z>4) Universe, we also include another 11 massive (10^10<M_ star<10^11 M_⊙ at z=6) galaxies extracted from a different set of MassiveFIRE simulations that stop at z=6, which are presented here for the first time. The latter were run with the same physics, initial conditions, numerics, and spatial and mass resolution, but were extracted from larger simulation boxes (400 Mpc/h and 762 Mpc/h on a side, respectively).

FIRE simulations successfully reproduce a variety of observed galaxy properties relevant for the present work, such as the stellar-to-halo-mass relation <cit.>,  the sSFRs of galaxies at the cosmic noon (z∼2) <cit.>, the stellar mass – metallicity relation <cit.>, and the submm flux densities at 850 μ m <cit.>.



 §.§ Predicting dust SED with SKIRT





We generate the UV-to-mm spectral energy distribution (SED) using the open source[SKIRT code repository: <https://github.com/skirt>] 3D dust Monte Carlo RT code SKIRT <cit.>. SKIRT accounts for absorption and anisotropic scattering of dust and self-consistently calculates the dust temperature. We follow the approach by <cit.> (see also ) to prepare our galaxy snapshots as RT input models.

Each star particle in the simulation is treated as a `single stellar population' (SSP). The spectrum of a star particle in the simulation is assigned using starburst99 SED libraries. In our default RT model, every star particle is assigned an SED according to the age and metallicity of the particle.

While our simulations have better resolution than many previous simulations modelling infrared and submm emission (e.g., ) and can directly incorporate various important stellar feedback processes, they are still unable to resolve the emission from HII and photo-dissociation regions (PDR) from some of the more compact birth-clouds surrounding star-forming cores. The time-average spatial scale of these HII+PDR regions typically vary from ∼5 pc to ∼800 pc depending on the local physical conditions <cit.>. Hence, in our alternative RT model, star particles are split into two sets based on their age. Star particles formed less than 10 Myrs ago are identified as `young star-forming' particles, while older star particles are treated as above. We follow <cit.> in assigning a source SED from the MappingsIII <cit.> family to young star-forming particles to account for the pre-processing of radiation by birth-clouds. Dust associated with the birth-clouds is removed from the neighbouring gas particles to avoid double-counting <cit.>.

We present in Section <ref> and <ref> the results from our default (`no birth-cloud') model. In Section <ref> we will show that none of our results are qualitatively altered if we adopt the alternative RT model and account for unresolved birth-clouds. 




Our RT analysis uses 10^6 photon packets for each stage. We use an octree for the dust grid and keep subdividing grid cells until the cell contains less than f=3×10^-6 of the total dust mass and the V-band optical depth in each cell is less than unity. The highest grid level corresponds to a cell width of ∼20 pc, , about twice the minimal SPH smoothing length. For all the analysis in this paper, we adopt the <cit.> dust model with Milky-Way size distribution for the case of R_ V=3.1. At FIR, the dust opacity can be well described by a power law, κ_λ∝ 0.05 (λ/ 870μ m)^-β  m^2/kg, where β≈2.0 (see the upper panel of Figure <ref>) is the dust emissivity spectral index <cit.>. Gas hotter than 10^6 K is assumed to be dust-free due to sputtering <cit.>.  We self-consistently calculate the self-absorption of dust emission and include the transient heating function to calculate non-local thermal equilibrium (NLTE) dust emission by transiently heated small grains and PAH molecules <cit.>. Transient heating influences the rest-frame mid-infrared (MIR) emission (80 μ m) but has minor impact on the FIR and (sub)mm emission <cit.>. SKIRT outputs T_ mw for each cell that is obtained by averaging the temperature over grains of different species (composition and size). A galaxy-wide dust temperature is calculated by mass-weighting T_ mw of each cell in the galaxies. At high redshift (z>4), the radiation field from the cosmic microwave background (CMB) starts to affect the temperature of the cold ISM. We account for the CMB by adopting a corrected dust temperature <cit.>


    T^ corr_ dust(z)=[T^4+β_ dust+T^4+β_ CMB(z)-T^4+β_ CMB(z=0) ]^1/(4+β),


where T_ CMB(z) = 2.73 (1 + z) K is the CMB temperature
at z.

For this study, we assume that dust mass traces metal mass in the ISM, and adopt a constant dust-to-metal mass ratio δ_ dzr=0.4 <cit.> for our fiducial analysis. We also try two different cases where δ_ dzr=0.2 and δ_ dzr=0.8, and throughout the paper, we refer to these two dust-poor and dust-rich cases, respectively. In the lower panel of Figure <ref>, we show the galaxy SED for the three models. L_ IR increases when δ_ dzr increases because a higher optical depth leads to more absorption of stellar light and more re-emission at IR.

SKIRT produces spatially resolved, multi-wavelength rest-frame SEDs for each galaxy snapshot observed from multiple viewing angles. For the analysis in this paper, SEDs are calculated on an equally spaced logarithmic wavelength grid ranging from rest-frame 0.005 to 1000 μ m. We convolve the simulated SED output from SKIRT with the transmission functions of the PACS (70, 100, 160 μ m), SPIRE (250, 350, 500 μ m), SCUBA-2 (450, 850 μ m), ALMA band 6 (870 μ m) and 7 (1.2 mm) to yield the broadband flux density for each band.

We show in Figure <ref> the result of running SKIRT on one of our galaxies. In particular we show a composite U, V, J false-color image with and without accounting for dust absorption, scattering, and emission. We also show the image of ALMA 1.2 mm flux density, total IR luminosity, dust surface density and temperature. It can be seen that the 1.2 mm flux density traces the dust mass distribution, while IR luminosity appears to be more localised to the high-temperature region, since it is expected to be sensitive to temperature (L∼ MT^4+β). The local intensity of radiation, the dust temperatures, and the dust density all peak in the central region of the galaxy.





§ UNDERSTANDING DUST TEMPERATURE AND ITS SCALING RELATIONS



In this section, we at first review the different ways of defining galaxy dust temperature that have been used in different observational and theoretical studies (Section <ref>), and compute the different temperatures for the MassiveFIRE sample (Section <ref>). We compare the calculated dust temperature(s) of the simulated galaxies with recent observational data (Section <ref>). Finally, we reproduce several observed scaling relations ( L_ IR vs. temperature, sSFR vs. temperature) with the simulated galaxies and provide physical insights for these relations (Section <ref>). 



 §.§ Defining dust temperature



Dust temperature has been defined in different ways by observational and theoretical studies. Here, we focus on four different possibilities, which we call mass-weighted, peak, effective, and equivalent dust temperature.



  
*Mass-weighted dust temperature T_ mw

T_ mw is the physical, mass-weighted temperature of dust in the ISM.  T_ mw is often explicitly discussed in theoretical studies where dust radiative transfer modelling is applied to the snapshots from the galaxy simulations, and dust temperature is calculated using LTE (for large grains) and non-LTE (for small grains and PAH molecules) approaches <cit.>. 



  
*Peak dust temperature T_ peak

The peak dust temperature is defined based on the wavelength λ_ peak at which the far-infrared spectral flux density reaches a maximum <cit.>


    T_ peak=2.90×10^3 μ m· K/λ_ peak.


The peak wavelength λ_ peak is commonly derived from fitting the SED to a specific functional form, for instance, a modified blackbody (MBB), see below. λ_ peak (and T_ peak) in practice depends on the adopted functional form as well as the broadband photometry used in the fit <cit.>.



  
*Effective dust temperature T_ eff

The effective temperature is obtained by fitting the SED with a parametrised function. The effective temperature is thus a fit parameter, and like T_ peak, depends on both the adopted functional form and the broadband photometry.

For most observed SEDs, the RJ side of the dust continuum can be well described by a generalised modified-blackbody function (G-MBB) of the form <cit.>  


    S_ν_0 (T)    = A (1+z)/d^2_ L ( 1-e^-τ_ν)B_ν(T)
       =  1-e^-τ_ν/τ_ν(1+z)/d^2_ Lκ_ν M_ dust B_ν(T)


where ν_ o is the observer's frequency, ν=ν_ o (1+z) is the rest-frame frequency,  τ_ν is the dust optical depth at ν[Throughout this paper, all ν and λ with no subscript stand for rest-frame quantities, while those with “o" are the observed quantities.], κ_ν is the dust opacity (per unit dust mass) at ν, B_ν(T) is the Planck function, A is the surface area of the emitting source and d_ L is the luminosity distance from the source. τ_ν is often fitted by a power law at FIR wavelengths,  τ_ν=(ν/ν_ 1)^β, where β is the spectral emissivity index and ν_ 1 is the frequency where optical depth is unity. Observational evidence has shown that the value of ν_ 1 can differ between galaxies <cit.>. In principle, ν_1 can be determined from SED fitting given full FIR-to-mm coverage <cit.>. However, in practice, it is often taken to be a constant, ∼1.5-3 THz ( λ_1=c ν^-1_1=100-200 μ m) <cit.>.

The Wien side of the dust emission is expected to be strongly affected by the warm dust component in the vicinity of the star-forming regions, which can significantly boost the luminosity of galaxy with only a small mass fraction <cit.>, knowing L∼ M T^4+β. Observations also show a variety of SED shape at MIR <cit.>. To better account for the emission at MIR, <cit.> introduced a simple (truncated) power-law component to Eq. <ref>, giving rise to a G-MBB with an additional power-law component (GP-MBB)


    S_ν_0 (T) = A (1+z)/d^2_ L [  (1-e^-τ_ν)B_ν(T) + 
    N_ pl ν^-αe^-(ν_ c(T)/ν)^2],


where N_ pl is the normalisation factor, α is the power-law index, and ν_ c is a cutoff frequency where the power-law term turns over and no longer dominates the emission at MIR. We allow N_ pl as a free parameter, fix α=2.5, and adopt the functional form of ν_ c(T) provided by <cit.>. The latter were constrained by fitting the observational data of a sample of local IR-luminous galaxies from the Great-Origins All Sky LIRG Survey <cit.>.

In the optically-thin regime (τ≪1), Eq.<ref> reduces to the optically-thin modified black body function (OT-MBB), <cit.>


    S_ν_0   = (1+z)/d^2_ Lκ_ν M_ dust B_ν(T)
       = (1+z)/d^2_ Lκ_870(ν/ν_870)^β M_ dust B_ν(T)
       =𝒞_ν(z) M_ dust B_ν(T)


where κ_870 is the opacity at 870 μ m (κ_870=0.05  m^2 kg^-1 for the dust model used in this work), ν_870=343  GHz, and 𝒞_ν(z) is a known constant for a given ν, κ_870, β, and z and d_ L.

The long-wavelength (λ200 μ m) RJ tail of the dust emission, where dust optical depth becomes low, can be well fit by the above equation. However,  Eq. <ref> is also frequently adopted to fit the full dust SED, including both the Wien and the RJ sides, especially by the studies in the pre-Herschel era, when not enough data is available to well cover both sides from the SED peak <cit.>. The single-T parameter in Eq. <ref> is then often referred to as the `dust temperature' of the galaxy. However, an effective temperature derived this way should be primarily understood as a fitting parameter and may not correspond to a physical temperature <cit.>. In particular, it differs in general from the mass-weighted temperature of dust in a galaxy. 





  
*Equivalent dust temperature T_ eqv


We define T_ eqv as the temperature that reproduces the actual IR luminosity for a given broadband flux (e.g., at 870 μm) and adopted parametrised functional form of the SED (e.g., OT-MBB). The value of T_ eqv typically depends on both the observing frequency band as well as the SED form (Section <ref>).

In the specific case of optically-thin dust emission, the specific luminosity, can be written as 


    L_ν, OT(T, M_ dust)    = 4π(1+z)^-1d^2_ L S_ν_ o
       =4πκ_ν M_ dust B_ν(T)


By directly integrating the above formula over ν, one obtains the total IR luminosity <cit.>


    L_ IR,  OT (T, M_ dust)    = ∫^∞_0 4π M_ dustκ_ν B_ν(T)  dν
       = 4π M_ dustκ_ν_1ν_1^-β(k_B T/h)^4+β(2h/c^2)
       ==Γ (4+β) ζ(4+β)
       = 𝒟 M_ dustT^(4+β),

where 𝒟(κ_ν_1, ν_1, β) is a constant and Γ and ζ are Riemann functions.

Combining Eq. <ref> and Eq. <ref>, T_ eqv can now be defined as the temperature satisfying

    L_ IR / S_ν_0 = 𝒟 T_ eqv^4+β/𝒞_ν(z) B_ν(T_ eqv).


In the RJ regime, where B_ν(T_ eqv)=2ν^2k_ BT_ eqv/c^2,

    L_ IR / S_ν_0∝ T^3+β_ eqv.



T_ eqv is therefore the temperature that one would need to adopt in order to obtain the correct IR luminosity and match the broadband flux density under the assumption that the SED has the shape of an OT-MBB function. Of course, the latter assumption is often a poor one and the actual SED shape can differ substantially from an OT-MBB curve. In this case, the equivalent temperature will be different from the mass-weighted dust temperature. Furthermore, the dust mass that is derived this way (via Eq. <ref> for a given T_ eqv and S_ν_0) will then differ from the actual physical dust mass.

In this paper, we compute T_ eqv based on Eq. <ref> using the actual integrated IR luminosities and 870 μ m (1.2 mm) flux densities unless explicitly noted otherwise. For equivalent temperatures based on G-MBB or GP-MBB spectral shapes, we numerically integrate Eq. <ref> and Eq. <ref> to obtain the IR luminosity for a given dust temperature and dust mass (analogous to Eq. <ref> for the OT-MBB case).



 §.§ The SEDs of simulated galaxies



In Figure <ref> we show example SEDs of a z=2 galaxy and a z=6 galaxy from the MassiveFIRE sample. We separately discuss z=2 and the z=6 galaxies because the observational strategies for the two epochs are usually different. For z=2, an IR-luminous ( L_ IR 10^12 L_⊙) galaxy may have both Herschel coverage at FIR as well as (sub)mm coverage from ground-based facilities ( SCUBA, ALMA and AzTEC). One can then derive the dust temperature (T_ peak or T_ eff) from the observed FIR-to-mm photometry via SED fitting. At z>4, the sources that have a good coverage of the SED peak (via Herschel surveys) are currently limited to higher IR luminosity ( L_ IR 10^13 L_⊙) and the majority are strongly lensed objects <cit.>. Meanwhile, the unprecedented sensitivity of ALMA has allowed us to detect the dust continuum of a growing population of galaxies at these epochs <cit.> that do not have detected Herschel counterparts. Most of these observations cover only a single band (typically at ALMA band 6 or 7). Physical properties, such as L_ IR and SFR, are thus often derived based on a single data point at (sub)mm <cit.>, by assuming a dust temperature for the object <cit.>. This approach is sensible if the adopted dust temperature is close to T_ eqv of the given galaxy (see section <ref>).



  §.§.§ Example: The SED of a galaxy at z=2


Figure <ref> shows the SED of a selected z=2 MassiveFIRE galaxy (upper left panel). This galaxy has L_ IR=3.0×10^12 L_⊙, SFR=210 M_⊙  yr^-1, M_ dust=5.4×10^8 M_⊙ and M_*=5.3×10^11 M_⊙ [Physical properties of the simulated galaxies reported in this paper are measured using the material within a 30 pkpc kernel around the DM halo centre,  the minimum gravitational potential.].

We calculate the PACS (70, 100 and 160 μ m) + SPIRE (250, 350 and 500 μ m) + SCUBA-2 (450 and 850 μ m) + ALMA (870 μ m and 1.2 mm) broadband flux densities from the simulated SED. We fit its FIR-to-mm photometry — assuming successful detection at every band, as we show in the left panel that the PACS/SPIRE fluxes of this galaxy are above the confusion noise limit (marked by the horizontal ticks) <cit.> and the submm fluxes are above the typical sensitivity limit of SCUBA-2 and ALMA — by a GP-MBB function (with λ_1=100 μ m, β=2.0 and α=2.5) using least-χ^2 method. N_ pl and T are left as two free parameters for the fitting. The best-fitting GP-MBB function is shown by the thick magenta line. The derived T_ eff is 31.4 K, which is similar to its mass-weighted temperature (T_ mw=29.1 K)[How well T_ eff in the best-fitting GP-MBB function approximates T_ mw depends on its parametrisation (see Section <ref>). For instance, increasing λ_1 in Eq. <ref> from 100 to 200 μ m changes T_ eff from 29.1 K to 48.2 K (see also Figure 20 of ).]. From the best-fitting GP-MBB function (and also the simulated SED), T_ peak is found to be 33.9 K. 

For demonstration purpose, we also show with the blue line the exact solution of the OT-MBB function, with T=T_ mw=29.1 K, M_ dust = 5.4×10^8 M_⊙, κ_870=0.05  m^2 kg^-1 and β=2.0. As expected, the OT-MBB function with a mass-weighted temperature is in very good agreement with the galaxy SED at long wavelength. For this galaxy, at λ=100-650 μ m (λ_ o=300 μ m-2  mm), the difference between the flux of the OT-MBB function and the simulated flux is within 10% (illustrated by the lower left panel). At shorter wavelength, the emission is more tied to the dense, warm dust component in the galaxy, which is poorly accounted for by this OT-MBB function with a mass-weighted temperature. Overall, the OT-MBB function accounts for ∼50% of L_ IR of the galaxy, and the discrepancy is largely due to the MIR emission.

We also show the effect of optical depth. In the upper left panel, the green line shows the analytic solution from a G-MBB (Eq. <ref>) function with the same M_ dust and T (T=T_ mw=29.1 K), but with a power-law optical depth that equals unity at rest-frame ν_0=1.5 THz, or λ=100 μ m. While the emission looks identical to the optical-thin case (blue line) at long wavelength (λ_ o>500 μ m), it appears to be lower at shorter wavelength when the effect of optical depth becomes important. The effect of increasing optical depth is that the overall light-to-mass ratio is lower and the emission peak wavelength is longer compared to the optically-thin case <cit.>.



  §.§.§ Example: The SED of a galaxy at z=6



Figure <ref> also shows the SED of a z=6 MassiveFIRE galaxy. This galaxy has lower L_ IR (3×10^11 L_⊙) and M_ dust (8×10^7 M_⊙) compared to the z=2 galaxy, but interestingly, it has similar T_ mw (30.7 K). The calculated flux densities at ALMA band 7 (S_ 870 μ m) and 6 (S_ 1.2 mm) are 0.44 and 0.23 mJy, respectively. Like the z=2 galaxy, an OT-MBB function (blue line) with M_ dust and T=T_ mw can well describe the emission of the z=6 galaxy at long wavelength (for this case, λ_ o>1.2 mm, or rest-frame λ>170 μ m), but it only accounts for ∼30% of L_ IR. A larger fraction of the total emission of this z=6 galaxy origins from the warm dust component.

To estimate L_ IR of a z=6 galaxy from S_ 870 μ m (or S_ 1.2 mm), one often needs an assumed SED function and an assumed T_ eqv for the adopted function. Since it is extremely difficult to constrain the details of SED shape at this high redshift, often a simple OT-MBB or GP-MBB function is used by the observational studies <cit.>. As an example, we fit the OT-MBB function to S_ 1.2 mm of the z=6 MassiveFIRE galaxy with varying T. We show in the right panel of Figure <ref> the OT-MBB function (with fixed β=2.0) that yields the simulated L_ IR with the light blue line. The derived T_ eqv for this function is 49.1 K. This is significantly higher than T_ mw, and as a result, the RJ side of the derived SED of this function appears to be much steeper than the simulated SED. It also poorly fits the simulated SED at wavelength close to λ_ peak. The derived T_ peak is therefore very different from the true T_ peak of the simulated SED. 

We also fit S_ 1.2 mm of this galaxy by a GP-MBB function (λ_1=100 μ m, β=2.0, α=2.5). We show the result for T=30.7 K (purple line), T=65.1 (magenta line) and T=80 K (salmon line). For T=T_ mw=30.7 K, we use the same normalisation of the power-law component as for the z=2 galaxy (upper left panel), so that the SED shape is similar between these two galaxies. For T=65.1 K and T=80.0 K, we use the best-fitting normalisation factor derived based on the local GOALS sample (see Table 1 of C12). We can see that the GP-MBB function appears to better describe the simulated SED shape compared with OT-MBB function, but in order to fit the simulated SED with reasonably good quality, a different choice of N_ pl and λ_1 is needed.  With T=T_ mw=30.7 K, the GP-MBB function under-predicts the simulated L_ IR (3×10^11 L_⊙) by 70%. Using T_ eqv, GP-MBB=65.1 K, this function leads to the right L_ IR. We also show the result for T=80 K, which over-predicts the L_ IR by about a factor two.

In conclusion, we find that an OT-MBB function with a mass-weighted dust temperature well describes the long-wavelength (λ 200 μ m) part of the dust SED, but it does not well account for the Wien side of the SED and leads to significant under-estimate of L_ IR. A GP-MBB function can provide high-quality fitting to the simulated SED with good FIR+(sub)mm photometry of galaxy. Using single-band (sub)mm flux density of z>4 galaxies, T_ eqv is very different from T_ mw of the galaxy. We will discuss T_ eqv for high-redshift galaxies, its evolution with redshift and its dependence on other galaxy properties in more details in Section <ref>.




 §.§ Comparing simulation to observation







Due to the high confusion noise level of the Herschel PACS/SPIRE cameras, most current observational studies on dust temperature at high-redshift are limited to the most IR-luminous galaxies in the Universe. For z=2, the observations are generally limited to L_ IR10^12 L_⊙. Applying the powerful stacking technique to the Herschel images, it is also possible to probe the fainter regime of a few 10^11 L_⊙ at z∼2 <cit.>. Yet another problem with the observational studies is the strong selection bias with flux-limited surveys, meaning that the selected galaxy sample is limited to increasing IR luminosity with redshift. It is therefore non-trivial to disentangle the dependence of dust temperature on redshift and that on other galaxy properties. Using simulated sample, we do not expect to have such problem. 

We start here by comparing the result of the MassiveFIRE sample at z=2 with the observational data from similar redshift. This is where the luminosity range of our simulated galaxies share the largest overlap with the current observational data. The selection methods of the quoted data are summarized in Table <ref>. At higher redshift, the observations are biased to higher L_ IR. In the following section, we will explicitly discuss the redshift evolution of dust temperatures with the MassiveFIRE sample. 

We present the result in Figure <ref>. In the upper panel, we compare the simulations with the observational data of which the (originally effective) dust temperature is derived using SED fitting technique and with MBB functions ( Eq. <ref>-<ref>), while in the lower panels, we show examples where the dust temperature of both the simulated and observation data is derived using the SED template libraries. In order to make fair comparison among different observations and with the simulation data, we convert all different T_ eff presented in the literature to T_ peak in the upper panel. T_ peak of the simulated galaxies are derived from the best-fitting GP-MBB function (Eq. <ref>, with λ_1 = 100 μ m, β=2.0 and α=2.5) to the FIR-to-mm photometry. 

In the upper panel, we show with the blue shaded block the data from the H-ATLAS survey <cit.>, that encompasses the high-redshift (1.5<z<2.0) Herschel-selected galaxies in the COSMOS field. The height of the block represents 1σ dispersion. We also explicitly show the z=1.5-2.5 objects from <cit.> (purple diamonds) and <cit.> (cyan asterisks), which are selected at 850 μ m from the deep SCUBA-2 Cosmology Legacy Survey <cit.> probing the Ultra Deep Survey (UDS) and the Extended Groth Strip (EGS) field, respectively. And finally, we present the stacked result by <cit.> (blue square), which is based on a high-redshift (⟨z⟩=2.23) sample extracted from the High-redshift Emission Line Survey (HiZELS) <cit.>, comprising 388 and 146 H_α-selected star-forming
galaxies in the COSMOS and UDS fields, respectively. And for purpose of reference, we show the binned data from <cit.> by grey filled circles and error bars, which encompasses an IR-selected sample at 0.1<z<2 selected from the COSMOS, GOODS-N and GOODS-S fields. We convert the effective dust temperature T_ eff presented in <cit.>, <cit.>, <cit.> and <cit.> to T_ peak. The relation between T_ peak and T_ eff for the fitting functions that are used by the four studies are plotted in Figure <ref>. <cit.> () adopt an OT-MBB function (Eq. <ref>) with fixed β=1.8 (β=1.5), while <cit.> () use a G-MBB function (Eq. <ref>) with fixed β=1.5 (β=1.6) and λ_1=100 μ m. From Figure <ref>, we can see that T_ eff presented in the four studies is higher than T_ peak. 



In the lower panels, we compare the simulated result with the observational data from <cit.> (left) and <cit.> (right), both of which fit the galaxy photometry to the empirical SED template libraries. In particular, <cit.> adopt the <cit.> SED template library and determine the temperature for each template by fitting their PACS+SPIRE flux densities with an OT-MBB function with fixed β=1.5 and then finding the T_ eff for the best-fitting OT-MBB function. Their sample comprises of near-infrared (NIR)-selected galaxies in GOODS-N, GOODS-S and COSMOS fields with reliable SFR, M_ star and redshift estimates.
The galaxies are binned in the SFR-M_ star-z plane and dust temperatures are inferred using the stacked FIR (100-500 μ m) flux densities of the SFR-M_ star-z bins with least-χ^2 method. We show the stacked result for their 1.7<z<2.3 redshift bin with the black filled dots in the lower left panel. For purpose of reference, we also show with the solid grey line the result of a lower-redshift bin (0.2<z<0.5) in the same panel.

In the lower right panel, we also compare the simulation to the observational data of <cit.>, of which the galaxy catalogue is based on the CANDELS survey <cit.>, a z=2-4 galaxy sample from the ALESS program <cit.>, as well as the local Herschel Reference Survey <cit.>. The temperature is derived by fitting the PACS+SPIRE photometry to the <cit.> SED template library, which is constructed based on the <cit.> library of elementary templates with an assumed power-law distribution of U. The G11 templates are a set of MIR-to-mm spectra emitted by a uniform dust cloud of 1M_⊙ when it is exposed to the <cit.> interstellar radiation field of a range of U. The temperature assigned to each <cit.> template of galaxy SED is the mass-weighted value of the G11 templates being used. We show in the lower right panel the result for the CANDELS sample with the black and grey filled circles. The black circles explicitly represent the objects at z=1.5-2.5. We also show with blue squares the result of the stacked SEDs for z=1.5-2.5 derived based on the PACS/SPIRE photometry in the CANDELS sample. The result of the ALESS sample at higher redshift (z=2-4) is shown with grey crosses. The black curve shows the scaling relation (T/ K)=5.57 (L_ IR/L_⊙)^0.0638 that is derived by <cit.> using the combination of the CANDELS, ALESS and HRS samples. 

For the simulated z=2 galaxies, we fit their PACS/SPIRE photometry to the <cit.> (as ) and <cit.> SED templates using least-χ^2 method and find the temperature associated with the best-fitting template SED as defined in the literature. In other words, the temperature of the MassiveFIRE galaxies is not the same in each of the three panels. The temperature derived following the <cit.> and <cit.> methods are on average 5.2 and 4.2 K higher than T_ peak, respectively. Comparing the simulated with the observational data, we find an encouragingly good agreement over the common range of L_ IR, with either the observational data derived using SED fitting technique (upper panel), or using SED templates (lower panels). And apart from that, T_ peak of the simulated z=2 galaxies appear to show no clear correlation with L_ IR in all three panels, at least at L_ IR 10^11 L_⊙. This is consistent with the recent finding by <cit.> that the mean dust temperature derived from the stacked SEDs of the three L_ IR bins of their z∼2 sample shows almost no correlation over the range of 1.5×10^11-1.5×10^12 L_⊙ (blue squares) and lies systematically above the mean temperature of galaxies at lower redshift (black line). This suggests that high-redshift galaxies do not necessarily follow a single, fundamental L_ IR-T scaling relation, which is typically derived using flux-limited observational data across a range of redshift but without much overlap of L_ IR among different redshift bins. We will also show in Section <ref> that the dust temperatures of our MassiveFIRE sample increase with redshift at fixed L_ IR from z=2 to z=6. <cit.> also report that the same redshift evolution extends to higher redshift (up to z=10) using a different suite of FIRE simulations.



The observational data shows nontrivial scatter, which is particularly clear in the upper and lower right panels. At L_ IR≈3×10^12 L_⊙, for instance, T_ peak (upper panel) is observed to be as low as ∼25 K and as high as ∼45 K. One possible reason is the intrinsic scatter of δ_ dzr.  We show in Figure <ref> the result for the dust-poor (δ_ dzr=0.2) and dust-rich (δ_ dzr=0.8) models in each panel. The former (latter) show ∼3 K increase (decrease) of dust temperature(s) compared with the fiducial model (δ_ dzr). This difference, however, still appears to be relatively smaller compared to the scatter of the observational data. A larger variance of δ_ dzr may lead to a larger scatter of temperature. Apart from that, another reason could be the variance of the conditions of the ISM structure on the unresolved scale ( compactness and obscurity of the birth-clouds embedding the young stars) could also contribute to the scatter. We will discuss more about the impact of sub-grid models later in Section <ref>. And finally, given that the Herschel cameras have fairly high confusion noise level, and it is rare that one galaxy has full reliable detection at every PACS/SPIRE+SCUBA band, we suggest that both factors can cause nontrivial uncertainty of observational result. Future infrared space telescope <cit.> spanning similar wavelength range and with higher sensitivity may help improve the constraint near emission peak and hence the observationally-derived dust temperatures.


We also note that z=2 MassiveFIRE galaxies appear to show higher dust temperature compared to the lower-redshift counterparts in the observed sample, with either the temperature derived using SED fitting (upper panel) technique or SED templates (lower panels). Observationally, how dust temperature evolves at fixed L_ IR (or M_ star) from z=0 to z=2 is still being debated <cit.>. Uncertainties can potentially arise from selection effects (surveys at certain wavelengths preferentially select galaxies of warmer/colder dust) <cit.> and inconsistency in derivation of dust temperature. The dust temperature of galaxies in this redshift regime (z<2) is beyond the scope of this paper.





 §.§ The role of dust temperature in scaling relationships
 


The scaling relationships of dust temperature against other dust/galaxy properties (such as total IR emission, sSFR and etc.) have been extensively studied in the past decade because of the significant boost of the number of detected high-redshift dusty star-forming galaxies by Herschel, SCUBA and ALMA. We now have statistically large sample for revealing and studying the various scaling relationships of dust temperature. Here in this section, we show the result of the MassiveFIRE sample at z=2-6, discuss the physical interpretation of the scaling relations and specifically examine how each scaling relation differs by using different dust temperatures (T_ mw vs. T_ peak). 



  §.§.§ S∝ MT (optically-thin regime)



As mentioned above, the long-wavelength RJ tail can be well described by a single-T OT-MBB function. This is a direct consequence of the rapid power-law decline of the dust opacity with wavelength as well as the fact that the coldest dust dominates the mass budget <cit.>. At very long wavelength, the flux is only linearly dependent on T in the RJ tail, and therefore the overall shape of the SED on the RJ side is largely set by the temperature of the mass-dominating cold dust. Hence, it has been proposed that the flux density originating from the optically-thin part of the RJ tail can be used as an efficient measure for estimating dust and gas mass (by assuming a dust-to-gas ratio) of massive high-redshift galaxies <cit.>. Given the high uncertainties of the traditional CO methods and their long observing time, this approach represents an important alternative strategy for gas estimate <cit.>. 

The RJ approach benefits from the effect of “negative K-correction". Eq. <ref> can be re-written as <cit.>


    S_ν_ o (T)    = (1+z)/d^2_ L2k_ Bκ_ν(ν/c)^2Γ_ RJ(ν_0,T,z) M_ dustT
       = ψ(z) Γ_ RJ M_ dustT


where Γ_ RJ is the RJ correction function that accounts for the departure of the Planck function from RJ approximate solution in the rest frame, and ψ(z) has the unit of mJy M^-1_⊙  K^-1. For given ν_ o, κ_ν(ν/c)^2 scales as (1+z)^4 (β=2.0). On the other hand, (1+z) d^-2_ L and Γ_ RJ decline with redshift. The former term roughly scales as (1+z)^-2, while how Γ_ RJ evolves with redshift depends on both ν_ o and T. The rise of κ_ν(ν/c)^2 with redshift can roughly cancel out or even reverse the decline of the other two components at z1, with typical T of galaxies and (sub)mm bands. For example, with T=25 K and ALMA band 6, ψΓ_ RJ stays about a constant from z=2-6, while with ALMA band 7, ψΓ_ RJ declines only by less than a factor of two over the same redshift range (see Figure 2 of ). (Sub)mm observations are therefore powerful for unveiling high-redshift dusty star-forming galaxies. In the RJ regime (hν≪ k_ BT), Γ_ RJ≈1 and S scales linearly to M_ dustT at a given redshift.

The RJ approach relies on an assumed dust temperature. The proper temperature, T, needed for inferring dust (and gas) masses can be obtained from solving Eq. <ref>, given S_ν_ o, M_ dust and z. This required T value is close to the mass-weighted dust temperature, for galaxies from z=2 to z=6, and with varying δ_ dzr, see Figure <ref>. The difference between these two temperatures is typically less than 0.03 dex. This again confirms that a single-T OT-MBB function well describes the emission from the optically-thin RJ tail. 

However, using T_ peak will lead to a poor constraint on M_ dust and therefore gas mass of galaxy. First of all, it is systematically higher than T_ mw, and therefore can cause systematically underestimate of M_ dust. Secondly, there seems to be no strong correlation between T_ mw and T_ peak by comparing the left and right panels. So even by using T_ peak to infer T_ mw will produce systematic error. We will discuss the discrepancy between T_ peak and T_ mw in more details in the later sections. Using other effective temperatures that have strong correlation with T_ peak will be problematic as well.



  §.§.§ The L_ IR vs. MT^4+β relation



 

The scaling relation L_ IR∝ M_ dustT^(4+β), which is frequently been adopted by many studies to probe and obtain useful physical insights for the star-forming conditions of the IR-luminous sources owing to its simplicity, is derived under the assumption of the optically-thin approximation (Eq. <ref>). 

The temperature in the above scaling relation is a measure of the luminosity per unit dust mass and often viewed as a proxy for the internal radiative intensity. Yet, it is not obvious how this temperature parameter ( ∼(L_ IR/M_ dust)^1/6) is related to the physical, T_ mw or the observationally accessible T_ peak.

We show in Figure <ref> the scaling relation of the light-to-mass ratio, L_ IR/M_ dust against T_ mw (left panel)  as well as T_ peak (right panel) for the MassiveFIRE sample at z=2-6, and we explicitly present the result for the fiducial (filled symbols), dust-poor (unfilled symbols) and dust-rich (semi-transparent symbols) cases. 

In general, galaxy having higher dust temperature (both T_ mw and T_ peak) emits more IR luminosity per unit dust mass. Focusing at first on T_ mw (left panel), we see that L_ IR/M_ dust of the MassiveFIRE galaxies appears to be systematically higher than from a simple single-T OT-MBB function (Eq. <ref>), which is indicated by solid black line in both panels. The offset (∼0.3 dex) between the simulated result and the analytic solution is due to the higher emissivity of the dense, warm dust in vicinity of the star-forming regions (see lower panels of Figure <ref>), which accounts for a small fraction of the total dust mass but has strong emission, and shapes the Wien side of the overall SED of galaxy.

With all the galaxies from z=2 to z=6, we find that L_ IR/M_ dust scales to ≈T^5.4_ mw. This is slightly flatter than the analytic solution derived using a single-temperature, optically-thin MBB function,  L_ IR, OT/M_ dust∝ T^6 (Eq. <ref>, with β=2.0). We understand the shallower slope as an optical depth effect. In the optically-thin regime (τ≪1), L/M∝ (1-e^-τ)/τ≈1, while in the optically-thick regime (τ≫1), L/M∝τ^-1 (Eq. <ref>). In the optically-thick regime, L_ IR/M_ dust therefore decreases with increasing τ. Galaxies of higher T_ mw are more dust-rich (Section <ref>) and their star-forming regions tend to be more optically-thick, resulting in a flattening of the scaling relation.

Comparing the dust-poor (dust-rich) models with the fiducial case, the median of T_ mw is higher (lower) by 0.84 (1.70) K. This is due to the optical depth effect. By reducing the amount of dust, the chance of receiving a short-wavelength photon increases because the optical depth from the emitting sources decreases. Therefore, dust is expected to be heated to higher temperature to balance the increased amount of absorption. Apart from that, δ_ dzr also mildly effects the normalisation of the L_ IR/M_ dust vs. T_ mw relation. The dust-poor (dust-rich) case shows about 0.13 (0.06) dex higher (lower) L_ IR/M_ dust, on the average, than the fiducial case, indicating a high (lower) luminosity emitted per unit dust mass. This is because a larger (reduced) mass fraction of the total dust is heated by (can actually “see") the hard UV photons emitted from the young stars due to the reduced optical depth <cit.>. This dust component can be efficiently heated to a temperature much higher than the mass-weighted average of the bulk <cit.>, and has a much higher L/M ratio than the rest.

T_ peak (right panel) also shows a positive correlation with L_ IR/M_ dust, although the strength of correlation is relatively weaker than that of T_ mw (ρ=0.81 vs. 0.91, where ρ is the Spearman rank correlation coefficient). Beside, T_ peak also shows larger scatter than T_ mw. The 1σ dispersion of L_ IR/M_ dust at fixed T_ peak is 0.21 dex, which is higher than 0.14 dex at fixed T_ mw. This means that T_ peak has relatively lower power to predict the luminosity-to-dust-mass ratio. Furthermore, T_ peak is also more affected by a change of δ_ dzr. The median T_ peak of the dust-poor (dust-rich) case is 2.49 (2.63) K higher (lower) than the fiducial model, which is more than the change of T_ mw with δ_ dzr. This is because T_ peak is more sensitive to the mass fraction of ISM dust that is efficiently heated to high temperature by the hard UV photons emitted from young stars <cit.>.




  §.§.§ L_ IR vs. T relation



The dust temperature vs. total IR luminosity is one most extensively studied scaling relations. We have shown in Section <ref> that our simulations have successfully produced the result at z=2 for galaxies that are in good agreement with the recent observational data at similar luminosity range. Here in this section, we focus on the evolution of dust temperature up to higher redshift. One major problem with the current observational studies on the T-L scaling is the selection effects of the flux-limited FIR samples that have been used to probe such relation. Higher redshift sample is biased towards more luminous systems <cit.>. How dust temperature evolves at fixed luminosity is still being routinely debated <cit.>. We present the result using our sample with L_ IR≈ 10^9-2×10^12 L_⊙ from z=2-6. For z=3-6, there is no current data available that we can make direct comparison to at similar L_ IR of our sample. Future generation of space infrared telescope, such as SPICA, can probe similar regime of IR luminosity at these epochs.

We present the temperature vs. luminosity relation of the MassiveFIRE galaxies at z=2-6 in Figure <ref>. In the upper and lower left panels, we show T_ peak vs. L_ IR and T_ mw vs. L_ IR relation, respectively. 

Focusing at first on T_ peak vs. L_ IR relation (upper left), we find a noticeable increase of T_ peak with redshift at fixed L_ IR, albeit with large scatter at each redshift. Looking at the most luminous galaxy at each redshift, we see that T_ peak increases from about 34 K at z=2 to ∼43 K at z=6 for the fiducial dust model (δ_ dzr=0.4). With all the luminous galaxies with L_ IR>10^11 L_⊙, we fit the evolution of T_ peak with redshift as a power law and obtained


    log (T_ peak(z)/25  K) = (-0.02±0.06)+(0.25±0.09)  log (1+z)


This result is in good quantitative agreement with the recent observational finding by <cit.> and <cit.>, although they use more IR-luminous sample at similar redshift range.

For each redshift, there is also a mild trend of declining T_ peak with decreasing L_ IR over the three orders of magnitude of L_ IR being considered. For instance, T_ peak of the z=6 galaxies at L_ IR=10^10 L_⊙ is about 32 K, which is about 10 K lower than the value at L_ IR=10^12 L_⊙, and is similar to the value of the brightest objects at z=3 and z=4. We find some faint objects at ∼10^10 L_⊙ whose T_ peak is as low as ∼20 K. We also note that the scatter of T_ peak could be very large at the faint end even with the simple fiducial dust model. At z=4, some objects could be as hot as ∼40 K, while some could be as cold as ∼20 K. This large scatter is mainly driven by the difference of sSFR among those galaxies, which we will discuss in more details in the following section.

With such large scatter, the correlation between T_ peak and L_ IR appears to be fairly weak. The Spearman correlation coefficient (ρ) of the T_ peak vs. L_ IR relation at individual redshift ranges from 0.46 to 0.70 at the redshifts being considered. For the z=2 sample, there is no noticeable correlation at L_ IR>10^11 L_⊙.

On the other hand, T_ mw exhibits a tighter correlation with L_ IR (lower left panel) (ρ ranging from 0.88 to 0.97), with an increase of the normalisation of the L_ IR-T_ mw relation with redshift. The increase of T_ mw with redshift at fixed L_ IR is clearly less prominent than T_ peak. At L_ IR≈10^12 L_⊙, for example, T_ mw increases from ∼27 K at z=2 to only ∼32 K at z=6. The CMB heating sets a temperature floor for T_ mw at the low luminosity end. 

The evolution of the T_ mw vs. L_ IR scaling is driven by M_ dust. At fixed L_ IR, galaxies at higher redshift have lower M_ dust. This can be clearly seen from the lower right panel, where we colour the same data as in the lower left panel by M_ dust of galaxy. There is clear sign of anti-correlation between T_ mw and M_ dust at fixed L_ IR <cit.>. Applying multi-variable linear regression analysis to the z=2-6 galaxies, excluding those that are strongly affected by the heating of the CMB background ( T_ mw T_ CMB(z)+5 K), we obtain the scaling relation


    log (L_ IR/10^10 L_⊙)   = (0.81±0.07) + (1.01±0.06)  log (M_ dust/10^7 M_⊙)
       +(5.40±0.36)  log (T_ mw/
        25  K),
     or  L_ IR   ∝ M_ dust T^5.4_ mw  (β=2.0).


It appears to be shallower than the classical L_ IR∝ M_ dustT^(4+β) (β=2.0 for our adopted dust model,  Fig. <ref>) relation derived based on the optically-thin approximation. We will discuss in Section <ref> about using this scaling relation to estimate M_ dust and T_ mw when only single data point is available at FIR-to-mm wavelengths.





  §.§.§ sSFR vs. T relation



The sSFR vs. dust temperature relation is one other frequently studied scaling relation which provide useful physical insights to dust temperature and is complementary to the L_ IR vs. temperature relation. 

In Figure <ref>, we show the relation of dust temperature against sSFR= SFR_ 20 Myrs/M_ star for the MassiveFIRE sample at z=2-6 in the left panels. We present the result for T_ peak and T_ mw in the upper and lower left panels, respectively. 

The dust temperatures are positively correlated with sSFR (ρ=0.68 for the sSFR vs. T_ peak relation and ρ=0.55 for the sSFR vs. T_ mw relation). Galaxies at higher redshift have, on average, higher sSFR, which is a direct consequence of the evolution of the star-formation main sequence. SFR is a proxy for the internal radiative intensity (most UV emission originates from the young stellar populations in the galaxies), and M_ dust is about linearly scaled to M_ star in the MassiveFIRE galaxies, the sSFR (∼SFR/M_ dust) can be viewed as a proxy for the total energy input rate per unit dust mass. It is therefore expected that to first order, dust temperature is positively correlated with sSFR of galaxies. This is indeed what we can see from both of the left panels of Figure <ref>. For instance, the z=2 galaxies (red) have a median sSFR of 3×10^-9  yr^-1 and median T_ mw=20 K (T_ peak=30 K). Both sSFR and dust temperature (both T_ peak and T_ mw), on average, increases with redshift. The z=6 sample (green) have a median sSFR of 2×10^-8  yr^-1 and median T_ mw=26 K (T_ peak=37 K). 

The correlation persists when focusing on each individual redshift.  In the middle panels, we show the result when both temperature and sSFR are normalised by the median value of the whole sample (T_ mw or T_ peak, sSFR) at each different redshift. With T_ peak (upper middle panel), the simulated galaxies, including all objects at z=2-6, exhibit a positive correlation (ρ=0.55) between starburstiness[The median sSFR at z=2, z=3, z=4 and z=6 of the MassiveFIRE sample are 2.1×10^-10, 5.8×10^-10, 8.7×10^-10 and 3.3×10^-9  yr^-1, respectively. SFRs are averaged over the past 20 Myrs.] ( SB=sSFR/sSFR(z)) and normalised T_ peak. The derived scaling relation (solid orange line) is in good qualitative agreement with the recent observations by <cit.> (dotted black line) and <cit.> (solid black line), despite that both studies include samples at lower redshifts (z<2) which our simulations do not probe. We also find that compared to L_ IR, T_ peak is more strongly correlated with sSFR at each given redshift, which is in agreement with the previous finding by <cit.> (see also ).

However, due to the inhomogeneity of dust distribution in galaxies and the complexity in star-dust geometry, the radiative energy emitted from the young stellar populations is not expected to evenly heat the ISM dust in the galaxy. Most of the UV photons are absorbed by the dense dust cloud in vicinity of the young star-forming regions, while the majority of the dust in the ISM is heated by the old stellar populations with more extended distribution, as well as the secondary photons re-emitted from the dust cloud near the young star clusters. For such reason, T_ peak is expected to be more sensitive to the emission from the warm dust component, which is more closely tied to the young star clusters, while T_ mw is determined by the cold dust component and therefore can be relatively less sensitive to the sSFR of galaxy than T_ peak.

This indeed can be seen from comparing the upper and lower middle panels of Figure <ref>. First of all, Δ T_ peak (T_ peak-T_ peak) shows a relatively stronger correlation with SB than Δ T_ mw (T_ mw-T_ mw). With all the z=2-6 MassiveFIRE galaxies, the Spearman correlation coefficient of the Δ T_ peak vs. SB scaling is ρ=0.55, while that of the Δ T_ mw vs. SB scaling is ρ=0.43. Secondly, over about two orders of magnitude of SB (∼0.1-10), the scaling relation with Δ T_ peak is relatively steeper, 


    Δ T_ peak∝ SB^ 8.61± 1.38  vs. Δ T_ mw∝ SB^ 4.03± 0.48.




This is because the UV photons from the young star clusters preferentially heat the dense dust cloud in the neighbourhood to high temperature, which boosts the MIR emission and helps shift the SED peak to shorter wavelength. However, the heating of the bulk of the dust is inefficient. The reason is that once the UV photons get absorbed and re-emit as FIR photons, the chance of them being absorbed by dust again becomes much lower as a consequence of the declining opacity with wavelength (κ_λ∝λ^-2) <cit.>. It is also interesting to note that both T_ peak and T_ mw are less correlated with SB when sSFR is averaged over longer period of time <cit.>. By averaging sSFR over a period of 100 Myrs instead of 20 Myrs, for example, ρ of the T_ peak (T_ mw) vs. SB relation declines from 0.55 (0.43) to 0.13 (0.22). 

We note that comparing to the recent observations of the z=2-4 star-forming galaxies <cit.>, the median sSFR of the MassiveFIRE galaxies is about 0.2-0.3 dex lower, but is still within the lower 1σ limit of the observational data <cit.>. This discrepancy of sSFR is commonly seen in the current cosmological galaxy simulations. A systematic increase of SFR will lead to more heating to the ISM dust and hence higher simulated dust temperatures, but will not affect the slope of the sSFR vs. T relation. The increment of T_ peak is estimated to be about 1.7-2.6 K according to the sSFR vs. T relation (Figure <ref>), which appears to be insignificant compared to the scatter of the observational data (Figure <ref>). It should also be noted that the impact of an increased sSFR on T_ peak can easily be offset by an increase of dust mass, which can potentially be driven by an increased δ_ dzr, dust opacity, or gas metallicities - all these properties are currently uncertain at high redshifts.

Finally, we show the relation between dust temperatures and M_ star in the right panels. Looking at the upper panel, it is clear that T_ peak has very weak correlation with M_ star. This again shows that T_ peak is strongly influenced by the emission from the warm dust that is associated with the recently formed young stars and does not have as strong correlation with the total stellar mass of a galaxy. In contrast, T_ mw is less sensitive to the variance of recent star-forming conditions and therefore shows relatively small scatter at given M_ star at each redshift. The normalisation of the T_ mw vs. M_ star relation increases with redshift, which is driven by the rise of SFR/M_ dust ( energy injection rate per unit dust mass). We also notice a slight increase of T_ mw with M_ star. This is owing to the decrease of M_ star/M_ dust with M_ star of the MassiveFIRE sample. As a result, SFR/M_ dust slightly increases with M_ star ( SFR/M_ dust∝ sSFR(M_ star/M_ dust) ∝ M^0.3_ star) at given redshift.






§ (SUB)MILLIMETRE BROADBAND FLUXES



A major problem for probing the dust properties in the high-redshift (z>4) is that most observations of dust emission at such high redshift are limited to a single broadband flux detected by ALMA band 7 or 6. Deriving infrared luminosities and hence SFRs of these z4 objects is very challenging without FIR constraints and depends highly on the assumed equivalent dust temperature for the flux-to-luminosity conversion. The same problem also applies to many faint ( below a few mJy) submm-selected objects at lower redshift (2<z<4) that do not have Herschel FIR coverage. Therefore, an accurate estimate of T_ eqv of the adopted SED function for different redshifts is critical.

In this section, we will analyse the T_ eqv distribution of galaxies at z=2-6 with the help of the MassiveFIRE sample. Specifically, in Section <ref>, we will examine the redshift evolution of T_ eqv and its dependence on δ_ dzr, offering a `cookbook' for converting between (sub)mm and L_ IR observations. In Section <ref>, we will compare T_ eqv with T_ mw and T_ peak, and provide a physical interpretation of this dust temperature.
 
  



 §.§ The flux-to-luminosity conversion



L_ IR is often extrapolated from a single broadband (sub)mm flux given the lack of additional submm or FIR constraints. A typical approach is to assume that the SED has an OT-MBB (or G-MBB) shape with a chosen value of the dust temperature parameter. However, as we have shown in Figure <ref> and discussed in Section <ref>, choosing a dust temperature parameter that is not compatible with the assumed functional shape of SED can result in significant biases for estimating L_ IR of a galaxy. By definition, this problem is avoided if the adopted dust temperature is chosen to be T_ eqv.

With the MassiveFIRE sample, we are able to predict the full dust SED for the high-redshift (z=2-6) objects covering over two orders of magnitude of IR luminosity (L_ IR≈10^10-10^12 L_⊙). We predict the observed flux densities at ALMA band 7 (S_870μ m) and band 6 (S_1.2 mm) given the SED and redshift as well as L_ IR. Many of these objects have S_ 870μ m (S_ 1.2 mm)0.1 mJy, which are over the 3σ detection limit of ALMA band 6 and 7 using a typical integration time of 1 hour. With the calculated S_870μ m (and S_1.2 mm) of each galaxy, we find the OT-MBB (with β=2.0) and GP-MBB functions (with β=2.0, λ_1=100 μ m, α=2.5 and the suggested value of N_ pl by C12), normalised to match their observed flux densities at both ALMA bands, that can predict their true L_ IR. By adjusting the temperature parameter in the fitting function to match both observed submm flux density and true L_ IR, we obtain T_ eqv, , the value of T necessary for obtaining an accurate estimate of L_ IR from the measured (sub)mm flux densities for each galaxy.

In Figure <ref>, we show the relation of L_ IR against S_870μ m (left panel) and S_1.2 mm (right panel) for the z=2-6 MassiveFIRE galaxies. For each redshift, we also show the expected L_ IR vs. S_ 870μ m (and S_ 1.2 mm) relation using the mean T_ eqv for galaxies above 0.1 mJy. The latter temperature is provided for the two different ALMA bands and for redshifts z=2-6. We present the results for OT-MBB and GP-MBB functional shapes.

There appears to be a clear trend of increasing T_ eqv with redshift, with either forms of fitting function (GP or OT-MBB) and with either ALMA band 6 or 7. This shows that a higher T_ eqv is typically needed for deriving L_ IR of galaxies at higher redshift. Using OT-MBB function, for example, the mean T_ eqv increases from 34.0 K at z=2 (red triangles) to 44.6 K at z=6 (green diamonds) for ALMA band 7. Applying the typical T_ eqv for z=2 to a z=6 galaxy will therefore lead to a significant underestimate of L_ IR.

For the same redshift, the normalisation of the L_ IR vs. S_ 870 μ m (S_ 1.2mm) relation depends on dust mass. We explicitly show in Figure <ref> the result for dust-rich and dust-poor models. At fixed observed broadband flux density, the L_ IR of dust-rich galaxies lies systematically below the fiducial model (vice versa for dust-poor galaxies). This result indicates that a galaxy of given observed (sub)mm flux density tends to have lower (higher) L_ IR if it contains more (less) amount of dust.

This finding can be understood as follows. By increasing the dust mass, both L_ IR and S_ 870μ m (S_ 1.2 mm) increase but the latter changes by a larger degree. Hence, the normalisation of the relation declines. The increase of S_ 870μ m (S_ 1.2 mm) is mainly driven by dust mass, as S_ 870μ m (S_ 1.2 mm) is linearly scaled to M_ dust (Eq. <ref>). On the other hand, the increase of L_ IR is due to enhanced optical depth — a larger fraction of UV photons gets absorbed by dust and re-emitted in the infrared/submm. A lower T_ eqv is therefore needed to account for the decrease of the normalisation of the L_ IR vs. S_ 870μ m (S_ 1.2 mm) relation with increasing dust mass. This anti-correlation of T_ eqv with δ_ dzr is more clearly shown in Figure <ref>.

We therefore provide a two-parameter fit for T_ eqv with δ_ dzr and redshift as predictor variables. Using all the z=2-6 objects with S_ 870 μ m>0.1  mJy, including the data for δ_ dzr=0.2-0.8, we perform a multiple linear regression analysis


    log (T_ eqv/25  K)=a+b  log (δ_ dzr/0.4)+c  log (1+z).


We present the best-fit regression parameters a, b and c for ALMA band 6 and 7, and for OT-MBB and GP-MBB functions in Table <ref>. These derived scaling relations are useful for converting a measured (sub)mm flux density into L_ IR, provided the redshift and dust-to-metal ratio of galaxy can be constrained.

The photometric redshift of the (sub)mm-detected galaxies can be determined when multi-band optical and NIR data are available, and the more accurate spectroscopic redshift can subsequently be determined if several atomic/molecular emission lines ( CO, CII, NII, OIII) are identified <cit.>. In contrast, δ_ dzr is more difficult to constrain from direct observation and is not yet well understood. Recent studies have reported differing results on how δ_ dzr depends on redshift and other galaxy properties <cit.>. We will discuss in more detail about the recent observations of δ_ dzr and the implication of the reduced δ_ dzr at high redshifts in Section <ref>.





 §.§ The equivalent dust temperature







T_ eqv depends on redshift and δ_ dzr in a clear and systematic manner, see Table <ref>. For the OT-MBB functional shape, for example, T_ eqv scales as ∝(1+z)^0.31δ_ dzr^-0.13 for ALMA band 7. This means that by applying a typical T_ eqv for z=2 to a z=6 galaxy would lead to an underestimate of L_ IR by a factor of ∼4 (Eq. <ref>). Also, at a given redshift, an order-of-magnitude increase of δ_ dzr corresponds to a ∼0.13 dex decrease of the best-fitting T_ eqv. This corresponds to a decrease of L_ IR by a factor of ∼4 (Eq. <ref>). Therefore, not taking the correlation of T_ eqv with redshift and δ_ dzr into account can potentially lead to significant biases in the L_ IR (and hence SFR) estimates.

The scaling T_ eqv∝(1+z)^0.31 (for band 7 and OT-MBB) is quantitatively similar to the one for T_ peak (Eq. <ref>), meaning that T_ eqv also evolves more quickly with redshift compared to T_ mw (see left panels of Figure <ref>). A natural question arises — what drives the evolution of T_ eqv with redshift?


To answer this question, we show in Figure <ref> the T_ eqv vs. T_ mw (upper panel) and T_ eqv vs. T_ peak (lower panel) relations of the MassiveFIRE sample at z=2-6. In this figure, T_ eqv  is calculated using an OT-MBB functional form (with fixed β=2.0) given a flux density at ALMA band 6. Using ALMA band 7 or a different form of MBB function results in qualitatively similar results and thus does not affect our conclusions.

It is clear from Figure <ref> that T_ eqv is more strongly correlated with T_ peak than T_ mw, either by looking at the z=2-6 sample as a whole, or each individual redshift. For each redshift, T_ peak scales approximately linearly with T_ eqv, with a high Spearman correlation coefficient ρ0.95. In contrast, the relation between T_ mw and T_ eqv is sub-linear and shows large scatter. As shown in the upper panel, galaxies with similar T_ mw can have very different T_ eqv (Δ T_ eqv>10 K) and thus a large range of L_ IR/S ratios (Eq. <ref>).

To understand the origin of the scatter in T_ eqv and fixed T_ mw, we selected two galaxies from the MassiveFIRE sample with similar T_ mw(≈ 30  K), one from z=2 and the other from z=6, and study their SEDs and their T_ eqv in more detail. The two galaxies are marked in both panels of Figure <ref> by yellow asterisks, and their SEDs are presented in Figure <ref>. The z=6 galaxy has T_ eqv, OT-MBB=49.1 K which is about 14 K higher than the z=2 galaxy.

Figure <ref> shows that the two galaxies have different SED shape at short wavelengths. The z=6 galaxy shows more prominent MIR emission due to its more active recent star formation. Its sSFR (=5.0×10^-9  yr^-1) is about one order of magnitude higher than that of the z=2 galaxy. Young star clusters in this high-redshift galaxy efficiently heat the dense, surrounding dust, which boosts the MIR emission and thus leads to a relatively high T_ peak (=44.6 K) to account for the more prominent MIR emission of this galaxy. Furthermore, the z=6 galaxy is less dust-enriched than the z=2 galaxy (having only 1/7 of dust mass), and its SFR/M_ dust ratio is roughly 4 times higher.

The increased SFR/M_ dust ratio would leave an imprint on the temperature of the diffuse dust if the heat budget of the young stars were evenly distributed in the ISM dust. However, the bulk of the diffuse cold dust is clearly not heated efficiently as the two galaxies have almost the same T_ mw (29.1 K vs. 30.7 K). A number of factors can influence how efficiently the bulk of the dust is heated, such as the spatial distribution of dust in galaxy and the optical depth in vicinity of the star-forming cores <cit.>. These conditions can be significantly different among galaxies and therefore T_ mw is not expected to be well correlated with T_ peak (see the upper right panel of Figure <ref>).  This example strongly indicates that a `two-phase' picture of ISM dust is needed to account for the discrepancy between T_ mw and T_ peak, see Figure <ref>.

Clearly, T_ eqv depends on the exact form of the MBB function and the observing frequency band. As is shown in Table <ref>, T_ eqv is higher at z=2 by 0.07 dex and it increases faster with redshift at z=2-6 when a GP-MBB function is used. Using the same MBB function, T_ eqv also appears to be slightly higher (by ∼0.05 dex at z=2-4 for an OT-MBB function) when a flux density is measured at ALMA band 7 than band 6. As T_ eqv depends both on the specific form of MBB function and the observing wavelength, T_ eqv should not be interpreted as a physical temperature but rather understood as a parametrisation of SED shape.  

It may appear reasonable to use sSFR as a predictor variable instead of (1+z), given that the former depends strongly on redshift (Figure <ref>) and is physically linked to the amount of hot dust in galaxies. The reasons for adopting (1+z) are two-folds. First of all, observationally, redshift of the (sub)mm-selected galaxies can be accurately determined through atomic/molecular emission lines, as discussed in Section <ref>. sSFR estimates, however, are uncertain because SFR derived based on the non-L_ IR indicators ( UV continuum and H_α flux) are uncertain due to the variation of the dust attenuation laws <cit.>. Secondly, the mapping between observed (sub)mm flux and rest-frame SED introduces an explicit redshift dependence on T_ eqv — as is shown in the lower panel of Figure <ref>, the normalisation of the T_ peak vs. T_ eqv relation declines with redshift using the same functional form and the observing frequency, indicating that a higher T_ eqv ( a steeper MBB function) is needed to derive L_ IR when the rest-frame observing wavelength gets closer to the emission peak. Therefore, the (1+z) term in Eq. <ref> accounts both (indirectly) for the cosmic time dependence of the sSFR and (directly) for the redshift of electromagnetic radiation.

Finally, adding M_ star as a predictor variable results in a regression coefficient for the M_ star term being consistent with zero. This means that our obtained fitting functions for T_ eqv do not depend on the selection function of M_ star of the MassiveFIRE sample, which can be different from that of the observations. Replacing the dependence on δ_ dzr by a dependence on M_ star or Z_ gas[Z_ gas is calculated using gas particles with temperature between 7,000-15,000 K and density above 0.5 cm^-3, which represent the nebular gas where the strong nebular emission lines originate <cit.>.] leads to a decreased goodness-of-fit for T_ eqv. 



§ DISCUSSION





 §.§ Deriving M_ dust



Many dust-enshrouded galaxies at high redshift (z>2) have been detected at (sub)mm wavelengths in the past years, thanks to the unprecedented sensitivity of ALMA. These (sub)mm-detected objects often lack a reliable measure of FIR photometry and many are extremely faint at UV/optical wavelengths <cit.>. A reliable estimate of their dust mass from full SED fitting is often not possible <cit.>.
 
In the optically-thin regime, the flux density in the RJ tail has a simple analytic form (Eq. <ref>), and M_ dust can be derived from the flux density given T_ mw (Section <ref>). However, it is difficult to constrain T_ mw of high-redshift galaxies when individual star-forming regions are not resolved.

Fortunately, we find that T_ mw does not strongly vary from galaxy to galaxy. This is noteworthy, given that our sample spans a wide range of cosmic time (z=2-6), stellar mass (M_ star = 10^9-10^12 M_⊙), sSFR (10^-10-10^-8  yr^-1), and IR luminosities (L_ IR=10^9-3×10^12 L_⊙). In particular, 68% ( 1σ) of the galaxies in our sample have mass-weighted dust temperatures T_ mw = 25±5 K,
corresponding to a 20% uncertainty of estimating the dust mass as the mass estimates scale only linearly 
with T_ mw, while 90% of our sample lies within T_ mw = 25±8 K (32% uncertainty of M_ dust). Our findings support the empirical approach of adopting a constant T_ mw = 25 K to estimate the ISM mass of high redshift galaxies via Eq. <ref> and δ_ dgr <cit.>.

While adopting a constant T_ mw is a good assumption to first order, and the only option if the (sub)mm flux density is measured at only a single wavelength, additional constraints on the SED may help to determine T_ mw and improve the accuracy of measuring ISM masses. Specifically, in Section <ref>, we show that T_ mw is well correlated with L_ IR and that the redshift evolution of the L_ IR vs. T_ mw relation is driven by the evolving dust mass. In fact, L_ IR, M_ dust and T_ mw follow a tight scaling relation (Eq. <ref>) for T_ mw≫T_ CMB. Hence, given S∝ M_ dustT_ mw, it should be possible to simultaneously infer M_ dust and T_ mw from a combined measurement of S and L_ IR.

Recent studies have shown that the broadband rest-frame 8 μ m luminosity, L_8, can be a rough tracer of L_ IR over a range of galaxies <cit.>. One important practical interest for using L_8 is that it will be easily accessible by the upcoming James Webb Space Telescope (JWST) for galaxies up to z∼3. The unprecedented sensitivity of the Mid-Infrared Instrument (MIRI) on board the JWST, covering the wavelength range of 5 to 28 μ m, will significantly enlarge the sample size of distant galaxies with measured MIR broadband spectroscopy. We thus propose to use L_8 to infer L_ IR for the (sub)mm-detected galaxies at z3 that have no constraint on SED shape near the emission peak <cit.>.

Hence, we propose to derive M_ dust (as well as T_ mw) of high-redshift galaxies by combining mid-infrared (e.g., from JWST) and far-infrared/submm (e.g., ALMA) data sets. Specifically, by combining Eq. <ref> and <ref>, we obtain


    log (M_ dust/M_⊙)   =1.23  log (S/ mJy)-0.23  log (L_ IR/L_⊙)+ℱ(z)
     or      M_ dust   ∝ (S/L_ IR)^0.23S,


where ℱ(z)=-0.85+1.23  log (ψ(z)Γ_ RJ) and ψ(z) has the unit of mJy M^-1_⊙  K^-1. Assuming that L_ IR=α L_8 <cit.>, we can rewrite the above equation as


    log (M_ dust/M_⊙)=1.23  log (S/ mJy)-0.23  log (L_8/L_⊙)+𝒢(z)


where 𝒢(z)=-0.23logα + ℱ(z). In general, Γ_ RJ is a function of T and Eq. <ref> &  <ref> need to be solved numerically.

It is important to note that α can depend on the variation in the detailed conditions of the star-forming regions. Recent observational evidence has shown that scatter in α can be driven by certain galaxy properties, such as sSFR, Z_ gas and compactness of IR-emitting regions <cit.>. These dependences on intrinsic galaxy properties can then translate to an apparent dependence of α on redshift and starburstiness (SB, see Section <ref> for definition) of galaxy. Therefore, to improve the accuracy of L_ IR estimates from L_8, one needs to rely on either a direct measurement or an observational proxy of these properties and/or galaxy redshift/SB. Furthermore, one should also caution the contribution of PAH molecules and AGN activity to the 8 μ m features <cit.>. This topic is beyond the scope of the present paper and we leave it to a future study.

According to Eq. <ref>, a factor of 2 uncertainty in L_ IR translates into ∼20% uncertainty in the derived dust mass, i.e., matches the intrinsic level of error of the constant T_ mw=25 K method <cit.>. Therefore, increasing complexity by deriving T_ mw and M_ dust from L_ IR and S will only be beneficial if L_ IR can be constrained to within a factor of 2 or better. 

Finally, we note that our sample does not include the most luminous submm galaxies that can have S_ 870 μ m fluxes of much higher than a few mJy <cit.> and, hence, we cannot rule out that T_ mw significantly exceeds 20-30 K in such objects. While submm-luminous galaxies are typically interpreted as having high SFRs, Eq. <ref> shows that submm fluxes are simply the product of M_ dust and T_ mw (Γ_ RJ is a weak function of T). Hence, as long as T_ mw is not significantly higher in these objects, a straightforward interpretation is that the most submm-luminous galaxies are those with the highest M_ dust. In fact, S_ 870μ m (S_ 1.2mm) is nearly doubled as M_ dust gets doubled (by comparing the dust-rich and fiducial cases in Figure <ref>), while T_ mw decreases by only ∼1 K ( 5%). This example also suggests that caution needs to be taken when directly converting (sub)mm flux densities to L_ IR  (and SFR), without taking into account how dust mass (and optical depth) alters the SED shape of galaxy ( , and see the lower panel of Figure <ref>, where we show how SED is altered by M_ dust without having a different SFR of galaxy).




 §.§ The increase of T_ eqv with redshift and its observational evidence



Adopting T_ eqv and an SED shape is another way to estimate the IR luminosity from submm fluxes, see Section <ref>. Hence, if T_ eqv is known, it is possible to use the approach described in the previous section to infer dust masses and mass-weighted temperatures. This could be a particularly useful approach at z3, where the potential MIR diagnostics redshift out of the wavelengths accessible by JWST.

Using the simplified functional forms of SED (an OT-MBB or GP-MBB), the obtained T_ eqv increases monotonically with redshift (Table <ref>). The typical T_ eqv of a z=6 galaxy is as high as 45-50 K for an OT-MBB function, which is significantly higher than the mean T_ mw (∼25 K) at this redshift. This result is consistent with what has been implied by some recent observational findings, including the unusual relationship between the IR excess (IRX≡ L_ IR/L_ UV) and the UV spectral slope (β) of the Lyman break galaxies (LBGs) at high redshifts.

Empirically, IRX is used as a proxy of total dust mass and β is a measure of dust column density. It has been found that galaxies between 0<z<4 follow a well-defined sequence on the IRX-β diagram, although there exists nontrivial scatter that depends on the stellar populations and the detailed dust and ISM properties <cit.>. However, recent observations of higher-redshift LBGs show evident IRX deficit <cit.>. A significant fraction of the z>4 UV-selected galaxies have no observed dust continuum at ALMA bands. For the ALMA-detected objects, their estimated IRX appears to be significantly lower than the value inferred from their measured β using the canonical IRX-β relations found by the local samples. The IRX deficit of the selected high-redshift LBGs is challenging to explain with the current dust attenuation models <cit.>.

Instead of the high-redshift populations having significantly different dust properties, an alternative solution is that they have a higher T_ eqv, which results in a higher derived L_ IR with a given observed (sub)mm flux density. For example, <cit.> report that among 330 LBGs in the Hubble Ultra Deep Field spanning the redshift range of z=2-10, only 6 were detected at 2σ at 1.2mm (band 6) by the ALMA Spectroscopic Survey (ASPECS). This is significantly lower than the number (35) extrapolated from the z=0 IRX-β relation based on their UV properties and by assuming a constant (equivalent) dust temperature of 35 K. The authors suggest that using a monotonic increase of (equivalent) dust temperature with redshift,  T∝ (1+z)^0.32 (an OT-MBB function is assumed), the number of detected sources can be consistent with the SMC IRX-β relation, which is derived based on the local metal-poor populations <cit.>. Encouragingly, this suggested redshift dependence of temperature well agrees with that of T_ eqv found by the MassiveFIRE sample (see Table <ref>). This indicates that a significant deviation of the dust properties of the high-redshift UV-selected populations is not needed for explaining their observed IRX deficit <cit.>.

Finally, T_ eqv should not be deemed equivalent to the mean intensity of the radiation field, <U>∝L_ IR/M_ dust <cit.>, despite that the latter has also been found to increase with redshift <cit.>. This is because T_ eqv is a parametrisation of SED shape, and it depends on both the assumed functional form of SED as well as the observing wavelength (Section <ref>), while <U> (∝ T^5.4_ mw, Eq. <ref>) represents the physical dust temperature. It is also important to note that since T_ mw not only evolves with redshift, but also shows a clear dependence on L_ IR at fixed redshfit (see the lower left panel of Fig. <ref>). The observed redshift evolution of <U>  therefore can potentially depend on the selection function of L_ IR. The selection bias caused by using a flux-limited sample (galaxies at higher redshift are confined to higher L_ IR) can lead to a steeper increase of <U>  with redshift than is measured at fixed L_ IR.



 §.§ The dependence of T_ eqv on δ_ dzr



Fig. <ref> shows that T_ eqv is anti-correlated with δ_ dzr at fixed redshift. Specifically, an order-of-magnitude decrease in δ_ dzr translates to ∼0.13 dex increase of required T_ eqv, corresponding to a factor of ∼4 increase of L_ IR (Eq. <ref>). As noted before (Section <ref>, <ref>), δ_ dzr directly affects the total dust mass (and optical depth) of a galaxy, thereby altering its SED shape. Hence, it can be one source of uncertainty in estimating L_ IR through T_ eqv.

Observationally, while δ_ dzr has been found to be fairly constant across a wide range of galaxies at low redshifts by different studies <cit.>, there is also evidence of a reduced δ_ dzr in the low-metallicity environments <cit.>. This can imply a lower δ_ dzr in higher-redshift galaxies since they are known to have lower metallicities than the low-redshift galaxies <cit.>. A direct implication of the decrease of δ_ dzr with redshift is that it can further mitigate the IRX deficit problem of the high-redshift LBGs, see Section <ref>. A reduced δ_ dzr leads to an additional increase of T_ eqv and hence a higher upper confidence limit of the IRX of the undetected objects, meaning that the dust properties of these high-redshift LBGs are more probable to be consistent with the canonical dust attenuation laws that are derived based on the low-redshift observations.

However, δ_ dzr of high-redshift galaxies is not yet well understood. Recent studies based on foreground absorbers towards γ-ray burst (GRB) afterglows and quasars (QSOs) as well as distant lensed galaxies have shown different trends of how δ_ dzr depends on redshift <cit.>, which can be due to selection effects and observational uncertainties <cit.>, as well as the different choice of δ_ dzr measures <cit.>. Moreover, because these studies are often limited to only a handful of galaxies, it becomes difficult to distinguish the explicit redshift dependence of δ_ dzr from the intrinsic correlation with the galaxy properties, although some studies based on the QSO damped Lyman-α absorbers (QSO-DLAs) find that the QSO-DLAs over a range of redshifts (z=2-6) follow a similar δ_ dzr-Z_ gas relation <cit.>.

Observationally, gas-phase galaxy metallicities (Z_ gas) can be derived using the ratios between (rest-frame) optical auroral and nebular line fluxes and with calibration on theoretical models <cit.>. However, this method can only be used for galaxies up to z∼3, above which the emission lines redshift out of the wavelengths of the current ground-based NIR spectrographs. An alternative method is to use the equivalent width of (rest-frame) UV absorption features, which has been used for galaxies up to z∼5, but is still limited because of the faintness of the features <cit.>.

To overcome these shortcomings, <cit.> have recently proposed a new method of using the (rest-frame) FIR [O_III] 88 μ m / [N_II] 122 μ m line ratio for probing the gas metallicities of galaxies at z>4, where both characteristic lines shift to the submm range that is accessible with ALMA. Using the previously reported FIR line measurements of a sample of local normal and star-forming galaxies by the Infrared Space Observatory <cit.>, <cit.> find that the derived galaxy mass-metallicity relation is consistent with the result derived using optical emission lines <cit.>. The gas metallicities of three z=2∼3 submm-luminous galaxies derived using Herschel measurements are also in good agreement with the high-redshift relationships previously derived by <cit.> and <cit.>, despite that the  stellar mass estimates of these obscured dusty high-redshift galaxies have large uncertainties. These results suggest that FIR emission lines could be promising tool for estimating Z_ gas of z>4 galaxies.

Therefore, if Z_ gas could be used as a predictor for δ_ dzr (which is currently uncertain), it will further help improve the accuracy of T_ eqv (and hence L_ IR) estimates at z>4, which is important since the MIR diagnostics for L_ IR are inaccessible by JWST at this epoch (Section <ref>). This will in turn help improve our constraints on the total obscured SFR in the early Universe, where currently only UV-based SFR estimates are available <cit.>.




 §.§ The sub-resolution structure of the birth-clouds







Observational evidence has indicated that young star clusters reside in dense dusty birth-clouds <cit.>. To check the uncertainty arising from potentially unresolved small-scale ISM structure, we have repeated the analysis presented in this paper with additional RT analysis by SKIRT as <cit.>, where we include a sub-grid model for birth-clouds embedding the young stars (our `alternative' RT model). We summarise the detail of this sub-grid model and the main results from this model in this subsection.

In brief, all the young star particle of a galaxy that has formed less than 10 Myrs ago is assigned a MappingsIII source SED <cit.>. MappingsIII SED templates are parameterised by the SFR and the metallicity of the star-forming regions, the pressure of the ambient ISM, the HII region compactness (log C), and the covering fraction of the associated PDR (f_ pdr). 

To explore how our results depend on this choice, the upper and lower panels in Figure <ref> show the overall SED of one of our galaxies for different values of log C and f_ pdr, respectively. As log C increases, the birth-clouds become more compact and the dust associated with the clouds attain higher mean temperature because of the stronger incident radiation onto dust grains. The source SED of this dust component (shown with dashed lines) shifts to shorter wavelength, and so does the overall SED of the galaxy. f_ pdr is a measure of the survival timescale of birth-clouds <cit.>. Increasing f_ pdr results in a larger fraction of the stellar emission being absorbed by dust in the birth-clouds, which results in more energy being re-emitted as IR light. The mean dust temperature of the birth-clouds, however, decreases. For the total emission of galaxy, a higher f_ pdr typically leads to higher L_ IR. Whether the emission peak of the overall SED shifts to shorter or longer wavelengths with f_ pdr, however, depends on the value of log C. 

The sub-grid model has minor impact on T_ mw of galaxies. T_ mw increases with log C at fixed f_ pdr, and decreases with f_ pdr at fixed log C. The reason is that the photons emitted from the birth-clouds are more energetic if the birth-clouds are more compact (higher log C) and less dust-obscured (low f_ pdr). But the resulting difference of T_ mw is typically no more than ±1 K (±5%) by exploring the parameter space of the MappingsIII model.

T_ peak, however, is more sensitive to the uncertainty of the small-scale ISM structure. In some cases, especially for strongly star-forming galaxies, T_ peak can differ by much as 10 K when the MappingsIII parameters are varied.  T_ peak is typically higher with increasing log C, and for low/intermediate (∼4.0-5.5) value of log C, decreases with  f_ pdr. For the z=2 MassiveFIRE sample, log C=6.5 (max) leads to a median T_ peak higher than log C=4.0 (min) by about 4 K with fixed f_ pdr=0.2, and f_ pdr=1.0 (max) yields a median T_ peak lower than f_ pdr=0 (min) by about 2.5 K with fixed log C=5.5. Uncertainty of the small-scale ISM conditions could introduce scatter in the observed T_ peak vs. L_ IR relation in addition to galaxy-by-galaxy variations of δ_ dzr (Figure <ref>).

Including the sub-grid birth-cloud model strengthens the correlation between T_ peak (and T_ eqv) and sSFRs of galaxies. By pre-processing starlight in birth-clouds, the range of the physical conditions surrounding star-forming regions is reduced. We note, however, that none of the trends reported in this paper change on a qualitative level by including or excluding the MappingsIII birth-cloud model. 



§ SUMMARY AND CONCLUSION




In this paper, we study dust temperatures of high-redshift galaxies and their scaling relationships with the help of cosmological zoom-in simulations and dust RT modelling. Our sample consists of massive (M_ star>10^10 M_⊙) z=2-6 galaxies extracted from the MassiveFIRE suite <cit.>, a set of cosmological hydrodynamic zoom-in simulations from the FIRE project <cit.>. The sample encompasses 18 central galaxies at z=2 and their most massive progenitors up to z=6, together with a disjoint set of 11 central galaxies at z=6. We generate FIR-to-mm broadband fluxes and spectra for our galaxy sample with SKIRT. 

We explicitly define and discuss four different dust temperatures that are commonly used in the literature, T_ mw, T_ peak, T_ eff and T_ eqv. T_ mw is the physical, mass-weighted temperature that can be extracted from RT analysis, but is often not easily accessible to observations. T_ eff and T_ peak are derived from SED fitting: T_ eff is the T parameter in the best-fit modified blackbody function and T_ peak is the inverse of emission peak wavelength. These two are the temperatures that are often adopted for analysing large statistical sample of galaxies by observational studies. And finally, T_ eqv is the temperature one needs to convert single (sub)mm data to total IR luminosity based on an assumed SED shape. 


The main findings of this paper are:


 
  * FIRE simulations together with RT processing successfully reproduce T_ peak of z=2, L_ IR10^11 L_⊙ galaxies, in good agreement with recent observations (Figure <ref>). The observational data shows large scatter, which may be driven by galaxy-to-galaxy variations of δ_ dzr as well as local variations in the physical conditions of unresolved birth-clouds embedding young star clusters (Section <ref>, <ref>).
 
 
  * T_ mw is only weakly correlated with T_ peak over z=2-6 (Figure <ref>). The former sets the slope of the RJ tail (Figure <ref>), and is the temperature needed for estimating dust and gas mass of distant galaxies (Figure <ref>). Using T_ peak, or T_ eff ( derived from full SED fitting), which is strongly correlated with T_ peak, can lead to a systematic bias/error of the derived dust/gas mass, and may lead to an inaccurate interpretation of the star-forming conditions in high-redshift galaxies (Section <ref>). 
 
 
  * T_ peak is well correlated with sSFR (ρ∼0.7) (Figure <ref>). Recently formed stars efficiently heat the dense, warm dust in the close vicinity of star-forming regions. The emission from this warm dust component boosts the overall dust SED at MIR, and helps to shift the emission peak to shorter wavelength (Figure <ref>). T_ mw is less well correlated with sSFR (ρ∼0.55) and the scaling relation shows a flatter slope (Δ T_ mw∝ SB^4.0 vs. Δ T_ peak∝ SB^8.7). The bulk of the cold diffuse dust is not as effectively heated as the warm dust component (Section <ref>).
 
 
  * T_ peak scales as (1+z)^0.25 at fixed L_ IR between z=2-6 driven by the increasing sSFR at higher redshift, which is consistent with recent observations (Section <ref>). T_ mw evolves only weakly with redshift at fixed L_ IR at z=2-6 (Figure <ref>).  
 
 
  * Of the galaxies in our sample, 68% have mass-weighted dust temperatures T_ mw = 25±5 K (Figure <ref>). This temperature range corresponds to an uncertainty of 20% in estimating M_ dust from a single submm band since the mass estimates scale only linearly 
with T_ mw. Furthermore, 90% of our sample lies within T_ mw = 25±8 K. Our findings support the empirical approach of adopting a constant T_ mw = 25 K to estimate the ISM mass of high redshift galaxies <cit.>.
 
 
  * T_ mw is well correlated with L_ IR at T_ mw≫T_ CMB at a given redshift (Figure <ref>). The normalisation of this relation evolves weakly with redshift but the slope does not change. At higher redshift, galaxies of the same L_ IR have higher T_ mw but lower M_ dust. Using the z=2-6 sample, we derive the scaling relation L_ IR∝ M^1.0_ dustT^5.4_ mw, which appears to be shallower than the classical L_ IR∝ M_ dustT^4+β (β≈ 2.0 for our adopted dust model) relation expected from the optically-thin assumption (Section <ref>).
 
 
  * We propose to use this scaling relation to derive M_ dust (and T_ mw) of high-redshift (sub)mm-detected galaxies, assuming that their L_ IR can be constrained, for example, via the mid-IR luminosity probed by the Spitzer telescope and the upcoming JWST. We showed that this method improves over the T_ mw=25 K approach if L_ IR can be constrained to within a factor of 2 or better (Section <ref>).
 
 
  * T_ eqv increases with redshift, meaning that a higher temperature is needed to convert observed (sub)mm broadband fluxes to L_ IR (and hence SFRs) of galaxies at higher redshift. T_ eqv is tightly correlated (ρ0.95) with T_ peak, a much stronger correlation than with T_ mw (Figure <ref>). In particular, two galaxies at different redshifts can have very different T_ eqv (Δ T_ eqv>10 K) but similar T_ mw (Section <ref>).
 
 
  * We find an anti-correlation between T_ eqv and the dust-to-gas ratio, δ_ dzr. Hence, at a given redshift, dust-poorer galaxies need, on the average, a higher T_ eqv for the (sub)mm-flux-to-IR-luminosity conversion. We express T_ eqv as a power-law function of δ_ dzr and (1+z), and perform linear regression analysis using the MassiveFIRE sample at z=2-6. The best-fit parameters of the scaling relation are provided in Table <ref>. We present the result for both ALMA band 6 and 7. We propose to apply the scaling relation to more accurately convert between (sub)mm flux and IR luminosity (and SFR) of high-redshift galaxies (Section <ref>).
 


To summarise our results, we find that the observationally-derived temperatures, in particular, T_ peak, generally differ from T_ mw. T_ peak shows a steeper slope and a stronger correlation with sSFR, and evolves more quickly with redshift compared with T_ mw. We also find that T_ eqv is more strongly correlated with T_ peak than with T_ mw.

The difference between T_ peak and T_ mw may be understood by a `two-phase' picture of ISM dust. T_ mw is set by the diffuse, cold dust component which dominates the total dust mass, while T_ peak is also influenced by the dense, warm dust component in the close vicinity of young star clusters. The former component is typically heated less effectively by young stars than the latter so that T_ peak and T_ mw are not well correlated with each other.

The increase of T_ eqv with redshift is consistent with recent observational evidence, including low number counts of (sub)mm sources in ALMA blind surveys <cit.> and the unusual IRX-β relation of high-redshift galaxies <cit.>. However, as we argue in this paper, the rise of T_ eqv with redshift is not simply a sign of dust being hotter at higher redshift, but it reflects the change in SED shape. In particular, higher T_ eqv is often a consequence of a more prominent MIR emission of galaxies at higher redshift, resulting from more active star formation. However, as T_ mw evolves only weakly between z=2 and z=6, the temperature of the majority of the dust component (∼T_ mw) does not significantly change despite the change in T_ eqv. In this sense, dust in galaxies with higher T_ eqv is not necessarily physically hotter.

In conclusion, dust temperature is important for estimating and probing key physical properties ( dust/gas mass, IR luminosity) and ISM conditions of high-redshift galaxies. A proper interpretation of dust temperatures and their scaling relationships requires taking into account the differences between temperatures derived from the SED shape and the physical, mass-weighted dust temperature. Upcoming facilities, such as JWST, SPICA and CCAT-prime, will significantly improve our capability of constraining key dust properties of galaxies in the distant Universe.



§ ACKNOWLEDGEMENTS


We are grateful to the referee for the careful review and insightful suggestions that have helped improve the quality of this manuscript. We acknowledge the comments from Andreas Faisst, Rob Ivison, Allison Kirkpatrick, Georgios Magdis and Ian Smail to an early version of this manuscript. We thank Caitlin Casey for a discussion on  before they were published on arXiv, which had helped the writing of this paper.  We are grateful to the technical guidance by Peter Camps and Maarten Baes. LL would like to acknowledge the stimulating atmosphere during the Munich Institute for Astro- and Particle Physics (MIAPP) 2018 programme the Interstellar Medium of High Redshift Galaxies, and especially thank the conversation with the programme coordinators. This research was supported by the MIAPP of the Deutsche Forschungsgemeinschaft (DFG) cluster of excellence “Origin and Structure of the Universe". RF acknowledges financial support from the Swiss National Science Foundation (grant no 157591). Simulations were run with resources provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research centre, proposal SMD-14-5492. Additional computing support was provided by HEC allocations SMD-14-5189, SMD-15-5950, by NSF XSEDE allocations AST120025, AST150045, by allocations s697, s698 at the Swiss National Supercomputing Centre (CSCS), and by S3IT resources at the University of Zurich. DK acknowledges support from the NSF grant AST-1715101 and the Cottrell Scholar Award from the Research Corporation for Science Advancement. CAFG was supported by NSF through grants AST-1517491, AST-1715216, and CAREER award AST-1652522; by NASA through grant 17-ATP17-0067; by STScI through grant HST-AR-14562.001; and by a Cottrell Scholar Award from the Research Corporation for Science Advancement. PFH was supported by an Alfred P. Sloan Research Fellowship,
NASA ATP Grant NNX14AH35G, and NSF Collaborative Research Grant #1411920 and CAREER grant #1455342. EQ was supported in part by a Simons Investigator Award from the Simons Foundation and by NSF grant AST-1715070. The Flatiron Institute is supported by the Simons Foundation.

mnras





