Cosmic Sands. II. Challenges in Predicting and Measuring High-z Dust Temperatures

In the current era of high-z galaxy discovery with JWST and the Atacama Large Millimeter/submillimeter Array, our ability to study the stellar populations and interstellar medium conditions in a diverse range of galaxies at Cosmic Dawn has rapidly improved. At the same time, the need to understand the current limitations in modeling galaxy formation processes and physical properties in order to interpret these observations is critical. Here, we study the challenges in modeling galaxy dust temperatures, both in the context of forward modeling galaxy spectral properties from a hydrodynamical simulation and via backwards modeling galaxy physical properties from mock observations of far-infrared dust emission. Using the simba model for galaxy formation combined with powderday radiative transfer, we can accurately predict the evolution of dust at high redshift, though several aspects of the model are essentially free parameters (dust composition, subresolution dust in star-forming regions) that dull the predictive power of the model dust temperature distributions. We also highlight the uncertainties in the backwards modeling methods, where we find the commonly used models and assumptions to fit far-infrared spectral energy distributions and infer dust temperatures (e.g., single temperature, optically thin modified blackbody) largely fail to capture the complexity of high-z dusty galaxies. We caution that conclusions inferred from both simulations—limited by resolution and post-processing techniques—and observations—limited by sparse data and simplistic model parameterizations—are susceptible to unique and nuanced uncertainties that can limit the usefulness of current high-z dust measurements.


Introduction
Like those of the gaseous interstellar medium (ISM), dust temperatures vary spatially within a galaxy, dependent on the proximity to star-forming regions, active galactic nuclei (AGNs), and the galactic disk.Within a single galaxy there exists a dust temperature distribution, with the mass-weighted average temperature typically lower than the luminosityweighted as the cold dust component dominates the mass while the warm dust dominates the luminosity.Characterizing the global "average" dust temperature is nontrivial; if a single dust temperature is used to describe a galaxy, the massweighted temperature is perhaps the most physical if concerned about characterizing the dust mass.By contrast, the luminosityweighted dust temperature is important for characterizing the intensity of star formation or AGN activity.
In the context of galaxy observations, the global dust temperature is most commonly measured from galaxy spectral energy distributions (SEDs) and is significantly dependent on the assumptions made regarding the shape of the far-infrared (FIR) SED (Yang & Phillips 2007).The most common model assumption is that the global dust emission can be well described by a single temperature modified blackbody; though in theory a multitemperature model can more accurately describe the dust temperature distribution in galaxies and yield a better fit to the FIR SED, especially for IR-luminous galaxies (Clements et al. 2010), the range and quality of photometry necessary to constrain such a model is challenging even for low-redshift galaxies (Casey 2012;Utomo et al. 2019) in addition to the inherent degeneracies between the multitemperature, optically thick, and variable spectral index solutions (Dunne & Eales 2001;Klaas et al. 2001;Rangwala et al. 2011).
The relationship between the SED-derived dust temperature and the real dust temperature distribution depends on both the model assumptions and the underlying dust properties (opacity, geometry, etc.; Liang et al. 2019).In some cases, e.g., the standard assumption of optically thin, single temperature sources, correlations between model parameters can be induced, adding to the aforementioned uncertainties and degeneracies.These model dependencies and degeneracies between model parameters (e.g., dust temperature, opacity, and spectral index) are well established (Blain et al. 2003;Shetty et al. 2009aShetty et al. , 2009b;;Kelly et al. 2012;Juvela et al. 2013) and can in principle be constrained with robust data spanning the range of FIR dust emission.
Yet there is a scarcity of FIR constraints at high redshift: dust continuum measurements are difficult to obtain except for the rarest, brightest galaxies, and are typically tied to observations of FIR emission lines such as [C II] or [O III] (e.g., Pavesi et al. 2016;Bakx et al. 2020Bakx et al. , 2021;;Harikane et al. 2020;Sommovigo et al. 2022a;Yoon et al. 2022) and will inevitably compound the issues described above.As cautioned in Spilker et al. (2016) and Drew & Casey (2022), the lack of constraints on the dust column density make constraining the FIR optical depth properties of these galaxies nearly impossible.This can lead to considerable biases in the dust and gas masses measured from broadband SEDs (Cochrane et al. 2022), as well as the potential evolution in dust temperatures across time (Drew & Casey 2022;Sommovigo et al. 2022a).Moreover, the applicability of certain assumptions that may be valid for some low-z galaxies (e.g., optically thin FIR emission) is relatively unknown for high-z galaxies; though these assumptions are typical in the literature, their uncritical use can clearly lead to further uncertainties.
In this paper, we aim to understand these uncertainties with the use of cosmological simulations.We employ the Cosmic Sands suite of early, dusty galaxies (Lower et al. 2023) to probe the efficacy of existing methods to infer dust properties from observational methods.The Cosmic Sands simulations offer the advantage of a self-consistent dust evolution model, which we couple with POWDERDAY dust radiative transfer (Narayanan et al. 2021), allowing us to forward model dust temperatures from the particle data and to backwards model dust temperatures as an observer would from the mock SEDs.Our goal is not to provide predictions for dust temperatures for high-z galaxies, but instead to use the Cosmic Sands galaxies as a testing ground for the methodologies used to infer dust properties in the literature.
This paper is organized as follows: in Section 2, we describe the Cosmic Sands suite of simulations based on the SIMBA galaxy formation model.In Section 3, we present an analysis of the radiative transfer derived dust temperatures of the Cosmic Sands galaxies and discuss the impact that the SIMBA dust model has on these results.In Section 4, we describe and compare the techniques for measuring dust temperatures with a modified blackbody model and the relationship the SEDderived temperature measurement has with the physical dust temperatures.Finally, in Section 5, we discuss the implications of our results for surveys of high-z dusty galaxies as well as the caveats to our galaxy formation model and post-processing radiative transfer.

Numerical Methods
In Lower et al. (2023), we presented the Cosmic Sands suite of zoom-in cosmological simulations to study the buildup of massive, dusty galaxies at Cosmic Dawn, focusing on the processes that impact efficient stellar mass buildup.Here, we use the Cosmic Sands galaxies to study the various methods of measuring dust temperatures and what drives the diversity in galaxy-averaged mass-weighted dust temperatures.Below, we briefly describe the SIMBA galaxy formation physics that Cosmic Sands is based on and the radiative transfer postprocessing done on the simulations snapshots to calculate dust temperatures.

Cosmic Sands and the SIMBA Galaxy Formation Model
Our galaxy sample is generated from the SIMBA galaxy formation model.SIMBA  (Davé et al. 2019) is based on the GIZMO gravity and hydrodynamics code (Hopkins 2015) and includes models describing heating and cooling, star formation, chemical enrichment, feedback from stellar winds, dust production and growth, and black hole (BH) accretion and feedback.We briefly summarize the key aspects of each model.
Physically, star formation occurs in dense molecular clouds; numerically, the rate of star formation is dictated by the density of H 2 divided by the local dynamical timescale, with gas particles below a density of n H = 0.13 cm −3 blocked from star formation.The H 2 fraction is modeled with the subresolution prescription of Krumholz & Gnedin (2011) depending on the gas-phase metallicity and local gas column density.
SIMBA uses the GRACKLE-3 library (Smith et al. 2017) to model radiative cooling and photoionization heating including a model for self-shielding based on Rahmati et al. (2013).The chemical enrichment model tracks 11 elements from Type Ia and II supernovae (SNe) and asymptotic giant branch (AGB) stars with yields following Nomoto et al. (2006), Iwamoto et al. (1999), and Oppenheimer & Davé (2006), respectively.
The subresolution model for stellar feedback includes contributions from Type II SNe, radiation pressure, and stellar winds.The two-component stellar winds adopt the massloading factor scaling from FIRE (Anglés-Alcázar et al. 2017) with wind velocities given by Muratov et al. (2015).Metalenriched winds extract metals from nearby particles to represent the local enrichment by the SNe driving the wind.The subresolution model for feedback via AGNs is implemented as a two-phase jet (low accretion rate) and radiative (high accretion rate) model.Thermal energy is injected into the surrounding ISM at high accretion rates, while BH-driven winds are produced at low accretion rates.
The dust model in SIMBA (Li et al. 2019) includes prescriptions for the production of dust by stellar sources, both in the remnants of Type II SNe explosions and in the metalloaded winds of evolved AGB stars.Dust can grow via metal accretion in the ISM and be destroyed via thermal sputtering, supernova shocks, and astration in star-forming regions.The dust is assumed to be at a fixed grain size of a = 0.1 μm, though we refer the reader to Li et al. (2021) for an evolving grain size distribution model implemented into the SIMBA galaxy formation framework.
The galaxies we use for this study are from the Cosmic Sands suite of cosmological zoom-in simulations presented in Lower et al. (2023).These simulations were selected from dark matter only runs and rerun with baryons at higher particle counts within the selected halos, thereby creating a highresolution "zoom-in" of the halo of interest.The Cosmic Sands simulations were preferentially selected as the most massive halos in 32 distinct boxes with volumes of 25 Mpc each.The baryonic mass resolution within the high-resolution region is 2.8 × 10 5 M e .At redshift z = 6, the halos span a range of masses (gas + stellar + dark matter) from 6 × 10 11 to 8 × 10 12 M e .We use CAESAR (Thompson 2014) to identify the initial dark matter halos and later the galaxies associated with each halo with a 6D (position-velocity phase space) friends-of-friends finder (6D FOF).The spatial linking length for gas and star particles is 0.0056 times the mean interparticle spacing, while the velocity linking length is 1 times the local velocity dispersion; gas and star particles that are bound in this way are associated with galaxies, from which the galaxy properties are calculated, which in turn are associated with their parent halos.

Dust Radiative Transfer
We post-process the simulated galaxies with 3D dust radiative transfer (RT) to generate UV-FIR SEDs for each massive galaxy.Following the schematic in Figure 1, we use the radiative transfer code POWDERDAY (Narayanan et al. 2021), which is built on YT (Turk et al. 2011), HYPERION (Robitaille 2011), and FSPS (Conroy et al. 2009;Conroy & Gunn 2010), to construct the synthetic SEDs.For each galaxy, an octree is constructed and gas particles (containing the SIMBA dust mass information) are smoothed onto the grid, which is refined to a maximum of 16 particles in each cell, with a minimum cell size ranging between 0.05 and 0.2 kpc.With FSPS, we generate the intrinsic SEDs for the star particles within each cell using the stellar ages and metallicities as returned from the cosmological simulations.For these, we assume a Kroupa (2002) stellar initial mass function and the MIST stellar isochrones (Choi et al. 2016;Dotter 2016).These FSPS stellar SEDs are then propagated through the dusty ISM.
Polycyclic aromatic hydrocarbon (PAH) emission is included following the Robitaille et al. (2012) model in which PAHs are assumed to occupy a constant fraction of the dust mass (here, modeled as grains with size a < 20 Å) and occupying 5.86% of the dust mass.The diffuse dust distribution is derived from the on-the-fly self-consistent model of Li et al. (2019).As SIMBA contains a "passive" dust model, with a single grain size (0.1 μm), we assume the extinction and emission models of Draine & Li (2007).This model includes the Weingartner & Draine (2001) grain size distribution and the Draine (2003) renormalization relative to hydrogen.Specifically, we assume the temperature-and frequency-dependent opacities and emissivities of the Weingartner & Draine (2001), though the emissivities are parameterized in terms of the mean intensity absorbed by grains, rather than the average interstellar radiation field as in the original Draine & Li model.
The radiative transfer propagates through the dusty ISM, with a resolution of 10 6 photon packets, in a Monte Carlo fashion using HYPERION, which follows the Lucy (1999) algorithm in order to determine the equilibrium dust temperature in each cell.We iterate until the energy absorbed by 99% of the cells has changed by less than 1%.We generate SEDs corresponding to 25 random viewing angles for each galaxy, allowing us to probe the relationship between apparent dust temperatures and the intrinsic dust properties of each galaxy.
We note here an important detail relevant to forward modeling of dust temperatures in the following section: the emissivities described above, which account for the imperfect blackbody (graybody) nature of dust grain emission in the IR, are calculated for the specific grain sizes and compositions in the adopted model.While the emissivities are known for each dust cell, there is not a straightforward way to predict the emissivity of the composite FIR SED as the emission also depends on the local temperature.Each grid cell containing dust will emit a blackbody with nonzero emissivity at the local temperature; the galaxy FIR SED is effectively the sum of these graybodies with an emissivity that will depend on the temperature distribution and the dust column density along the line of sight.

Forward Modeling Dust Temperature
We derive the mass-weighted (T mass-weighted ) and luminosityweighted (T lum-weighted ) dust temperatures for the Cosmic Sands galaxies using POWDERDAY radiative transfer, with the full temperature distributions shown in Figure 2. Because Cosmic Sands uses the SIMBA galaxy formation model that includes dust evolution, we do not have to make assumptions about the distribution of dust or its scaling with the gas or metal mass.The mass-weighted temperatures are calculated in each POWDERDAY grid cell i occupied by dust.We calculate T mass-weighted as where T i and M i are the dust temperature and mass in each cell, respectively.
To compute the luminosity-weighted temperatures we instead weight by the luminosity in each cell i, assuming the dust emits as a modified blackbody at the temperature in the cell with absorption and emissivity values from Weingartner & Draine (2001): where B λ is the blackbody function.We note here that neither the mass-weighted nor the luminosity-weighted temperatures are necessarily measurable quantities, as radiative transfer effects can bias temperatures measured from dust luminosities.To demonstrate this, we also plot in Figure 2 the apparent dust temperatures corresponding to the peak of the FIR SED (T peak ) in the solid black lines (calculated from the spectral flux density F ν ).For galaxies with lower dust mass, the peak temperature is typically similar to the luminosity-weighted temperature but the relationship between the peak temperature and the temperature distributions of dustier galaxies is more complex.As we will discuss in subsequent sections, namely Section 3.1, this is primarily driven by optical depth effects in the FIR becoming more significant in dustier galaxies, including complex star-dust geometries with highly opaque dust column densities.We collapse the Cosmic Sands dust temperature distributions into their luminosity-weighted values in Figure 3, where we plot T lum-weighted as a function of the total stellar mass, dust mass, and star formation rate (SFR) surface density.The SFR surface density is calculated using the CAESAR FOF defined stellar half-mass radius, computed from the radially sorted cumulative mass of the star particles tagged by CAESAR as being part of the galaxy structure.The radii range from 0.08 to 1.4 kpc.The dust temperature-mass relations show that more massive galaxies have warmer dust, but there is increased scatter at the high-mass end, where galaxies with M * 2 × 10 10 M e , M dust 8 × 10 7 M e have dust temperatures that vary by 40 K.The relation with SFR surface density shows a clearer trend, implying that luminosity-weighted dust temperature is more sensitive to density than to the integrated quantities; the high-density starburst galaxies host higher surface density and warmer dust reservoirs which increase the luminosity-weighted temperatures.
In Figure 4, we compare the mass-weighted, luminosityweighted, and peak dust temperatures.We find that luminosityweighted temperatures are warmer than the mass-weighted, indicating that while the most luminous dust is the warm dust, most of the dust mass is colder, as expected.The ratio of T lum-weighted /T mass-weighted ranges from 1.6 to 2.1.The massweighted temperatures tend to increase with increasing dust mass.Also demonstrated in Figure 4 is the (lack of) correlation between the apparent temperature, T peak , and T lum-weighted and T mass-weighted .For galaxies with less dust mass, the luminosityweighted and peak temperatures agree relatively well, with some scatter driven by the temperature distributions unique to each galaxy.But the peak temperatures are much lower compared to T lum-weighted for the more dust rich systems, due to the complex interplay between local dust temperature and emissivity properties as will be discussed in the following section.Thus, while the dust temperature distribution of a single galaxy can be represented by a single characteristic temperature, this temperature is (1) not unique to that temperature distribution and (2) not always equivalent to the other characteristic temperatures as each temperature is sensitive to a specific part of the entire dust temperature distribution.

The Impact of FIR Optical Depth
As shown in Figure 4, the relationship between the characteristic dust temperatures (T lum-weighted and T mass-weighted ) and the apparent temperature (T peak ) are complex for galaxies with large dust reservoirs, primarily due to radiative transfer effects and dust optical depth.Dust thermal emission in the FIR is typically assumed to be optically thin, or at least transitioning from optically thick to thin around 100-200 μm; this assumption is adopted extensively in the literature to model FIR SEDs of both UV-selected and submillimeter-selected galaxy samples (Swinbank et al. 2014;Scoville et al. 2017;Harikane et al. 2020;Sommovigo et al. 2022a).Here, we argue that this assumption may be invalid for dusty galaxies at high redshift, with significant implications for the observed dust properties.
First, in the left panel of Figure 5, we show the infrared luminosity of the Cosmic Sands galaxies as a function of luminosity-weighted temperature with galaxies color coded by dust mass.The luminosities are calculated from the SEDs generated at a single random viewing angle.For galaxies with dust masses <2 × 10 7 M e , IR luminosities increase in step with increasing temperatures as is expected for optically thin emission, albeit with some scatter driven by the distribution of dust temperatures within each galaxy. 6But this correlation plateaus at larger temperatures and dust masses, implying that the increase in luminosity with increasing temperature does not occur as efficiently in these starbursty, dust-rich systems and could be driven by the fact that these particular galaxies are optically thick at FIR wavelengths.
We demonstrate this with two example galaxies in Figure 6, where we show the dust optical depths for a face-on and edgeon orientation at 100 μm over scales of 5 pc.The optical depths, τ i = ∫κ ρ i ds, are calculated in each POWDERDAY cell i with dust density ρ i along (the viewing angle dependent) line of Points are color coded by galaxy dust mass. 6If the dust in a galaxy were at uniform temperature, we would expect the global IR luminosity to be highly correlated with the apparent temperature (and in the optically thin limit, we should expect the apparent temperature to be very close to the intrinsic luminosity-weighted temperature, as evidenced by Figure 5).But because in all galaxies there is an intrinsic distribution of dust temperatures, there will be some scatter in the T luminosity -LIR relationship, since the FIR SED will essentially be the sum of all of the modified blackbody spectra from the dust at different temperatures.sight ds from −5 to 5 kpc with respect to the center of mass of the galaxy, assuming absorption cross sections of κ = 40.95cm 2 g −1 at 100 μm (Li & Draine 2001;Weingartner & Draine 2001;Draine 2003).In the top panels, which shows a galaxy with M dust = 1.2 × 10 9 M e at z = 4.5, we find that when viewing the galaxy face-on, the central 0.1-0.5 kpc of the galaxy is opaque at 100 μm, with τ ∼ 25 in the densest region.When viewed edge-on, the thin disk has τ ∼ 30, with the thick disk marginally opaque spanning τ ∼ 1-5.In the bottom panel a galaxy at z = 6.5 (left) and z = 5.5 (right) with a dust mass of M dust ∼ 6 × 10 8 M e has dust opacities of 50−100 in the central 10-100 pc.
Similarly high dust opacities are found in compact-obscured nuclei (CONs; Aalto et al. 2019;Falstad et al. 2021) and ULIRG systems with CONs at low-z including Arp 220, whose FIR/submillimeter emission does not become optically thin until ∼2.5 mm (Scoville et al. 2017).The spatial extent of these high dust column densities vary depending on the source: the compact nucleus of the nearby IR-luminous galaxy IC 860 is opaque at millimeter wavelengths on scales of <10 pc (Aalto et al. 2019), while the disks of Arp 220 are opaque at 100 μm (2.6 mm) on scales of 240 pc (100 pc; Rangwala et al. 2011;Scoville et al. 2017), comparable to the highest optical depths of the Comic Sands galaxies at 500 pc.
At high redshift, interferometric imaging of IR-luminous galaxies have revealed opaque structures on kiloparsec scales.For instance, SMA and LABOCA imaging of the Cosmic Eyelash (SMMJ2135-0102; Swinbank et al. 2010) resolved optically thick clumps embedded in the more diffuse ISM, with sizes estimated at ∼500 pc; the average ISM column density of the system on the whole is estimated at ∼10 24 cm −2 over a 1 kpc radius (Danielson et al. 2011).Similarly, Atacama Large Millimeter/submillimeter Array (ALMA) detections of [C II] in the SPT SMG 0311-58 resolved clumpy structure on scales of 0.6-1.3kpc with [C II] emission suggesting high column densities of gas (0.8-12 × 10 4 M e kpc −2 ; Spilker et al. 2022).Lastly, dust properties inferred from NOEMA and ALMA detections of CO in several compact high-z starburst galaxies favored optically thick 100 μm emission (τ ∼ 1.5-4) over kiloparsec scales (Jin et al. 2022).
Explicit predictions for large FIR dust opacities from hydrodynamical simulations are sparse, but Yajima et al. (2014) report UV extinction toward individual stars in z = 3 Lymanbreak galaxies modeled with GADGET-3.For some stars in the most massive halos, they calculate optical depths in excess of 100 −1000 at 1700 Å, assuming the graphite and silicate dust model of Draine & Lee (1984).Though the pencil beam estimates do not compare spatially to the values calculated for the Cosmic Sands galaxies, some sightlines toward the Yajima et al. (2014) halos could be opaque in the FIR if the UV extinction is >1000.Similarly, Lovell et al. (2022) estimate the UV extinction along individual sightlines toward galaxies at z = 2 from the 100 Mpc SIMBA simulation.They find (using the column density of dust as a proxy for the dust emission along that line of sight) that some galaxies have UV opacities in excess of 1000, again implying that these systems could also be opaque in the FIR.All together, these findings suggest that while the column densities of the Cosmic Sands galaxies may be extreme, they are comparable to other dusty, star-forming systems at both low and high redshift.
The consequence of these high FIR dust opacities is that, observationally, only the dust at the τ  1 surface is detected, which can bias the measured dust luminosity, temperature, and mass of these systems if the optical depth is not accounted for (e.g., Cortzen et al. 2020;Fanciullo et al. 2020;Jin et al. 2022).Additionally, these large dust opacities in the Cosmic Sands galaxies explain the slowing trend between the luminosity and temperature for the most dust rich systems in left panel of Figure 5: the hottest dust in these galaxies is obscured by cooler dust in the foreground from which the IR photons we measure originate.Returning to Figure 5, we show this explicitly in the right panel, where we plot IR luminosity versus T lum-weighted but this time color coded by the maximum dust optical depth at 100 μm calculated in cells of 0.1 kpc. 7The galaxies whose luminosities are not as high as their luminosity-weighted temperatures would suggest have dust opacities >1, implying that we cannot see emission from the densest, hottest dust.In combination with Figure 3, higher-density starburst galaxies have more intense UV/optical radiation fields and have warmer, higher-density dust with higher optical depths.As we will explore in the following sections, one of the major simplifying assumptions made when fitting a modified blackbody function to galaxy SEDs is that the emission is optically thin at long wavelengths, which could severely limit the ability to accurately infer dust temperatures in the most dust rich systems.

Orientation Angle Bias
For both observed galaxies (e.g., Nagaraj et al. 2022) and galaxies from hydrodynamical simulations (e.g., Lovell et al. 2022), it has been shown that variations in measured galaxy properties (e.g., luminosity, UV-NIR dust attenuation) are significantly impacted by the galaxy's morphology and orientation with respect to the observer.We briefly make note of this in the context of biases in measuring dust properties from FIR SEDs.
To demonstrate this, we generated SEDs at 25 different viewing angles for each Cosmic Sands galaxies.Figure 7 now shows the dependence of the IR luminosities and peak temperatures on the orientation angle of each galaxy.The maximum spread in luminosity (0.33 dex) and temperature (8 K) are shown in each panel.What these plots demonstrate is the anisotropic nature of the IR emission, and the assumption of isotropic emission neglects the effects of star-dust geometry and galaxy morphology.The consequences of this for observations of real galaxies is that the IR luminosity and dust temperatures measured from galaxy SEDs are subject to an uncertainty that is not trivial to quantify.

Backwards Modeling Dust Temperature
We now approach modeling dust temperatures from the opposite side: inferring dust properties from galaxy SEDs.In practice, this is a relatively challenging task, given the degeneracies between FIR luminosity, dust temperature, and dust mass that drive variations in the shape of the FIR SED that can be difficult to constrain with limited data.As such, several approaches to determine dust temperatures exist in this regime, all of which depend on the assumed shape of the FIR SED and dust optical properties from the UV to FIR.We focus our analysis on the most common method for fitting FIR SEDs of galaxies at both low and high redshift: fitting a modified blackbody function to broadband FIR data. 88 Of course, there are other models that have been successful at modeling thermal dust emission on both molecular cloud scales and galactic scales, including the variable interstellar radiation field models of Dale & Helou (2002) and Draine & Li (2007), which effectively allow for a distribution of dust temperatures due to the differential heating of dust grains.But, like the modified blackbody model, the dust properties inferred from these models are dependent on the model assumptions, namely, the nature of the interstellar radiation fields and the chosen dust grain properties.

Model Definitions
If we assume modified blackbody emission from a source with a single temperature (Hildebrand 1983), T MBB is the dust temperature that fits the SED in the following, neglecting any distance/redshift terms: where B ν is the Planck function, τ ν is the frequency-dependent optical depth, and ò ν is the emissivity, which accounts for the "imperfect" blackbody (graybody) nature of dust grain emission in the FIR.We can simplify this equation since τ is proportional to ò and assume the typical frequency-dependent power-law parameterization for τ to get9 : where the optical depth is defined in such a way to be unity at ν 0 and the power-law slope is set by β, the spectral emissivity index.For an individual grain of known size and composition the emissivity can be measured, but in practice the emissivity of the entire dust population is difficult to determine due to dependencies on the dust temperature distribution and line-ofsight column density.The circular dependence of β on the local dust properties (temperature, and grain size and composition) and of T MBB on β when fitting SEDs alludes to one of the primary difficulties in constraining dust temperatures.Because a modified blackbody fit alone does not typically describe the mid-IR excess found in the Wien's side of the thermal emission peak due to the presence of warm dust, the short wavelength FIR data (λ < 50 μm) are typically fit with a power law: where α is the slope of the MIR power-law component and λ c is the wavelength where the SED transitions from power-law dominated to modified blackbody dominated.In both cases (Equation ( 4) MBB or Equation (5) MBB + power law), T MBB will correspond to the temperature of the dust at the τ ∼ 1 surface.
The spatial variation of the underlying distribution of dust temperatures in a galaxy can necessitate the need for additional flexibility that a single temperature modified blackbody model does not allow; as shown in Figure 2, the temperature distributions could be most accurately described by multiple temperature components (e.g., cold, warm, hot) that dominate the emission at different wavelengths.For galaxies with wellconstrained FIR SEDs, the best route would be to fit the data with a multitemperature component model (Casey 2012) instead of the more commonly used single temperature modified blackbody model above.But as demonstrated by Rangwala et al. (2011), while the single temperature model neglects the details of the spatial variation of dust temperatures within a galaxy, attempts to use a two-temperature optically thin model to fit SEDs of an extremely optically thick galaxy can result in erroneous estimates of other dust properties (e.g., dust masses that do not align with estimates from molecular gas constraints even though the inferred dust temperatures are physically reasonable).Thus, it is also important to understand the degeneracies between multiple temperature components and optical depth, as one can mimic the other's SED fingerprints (Dunne & Eales 2001).
We note here that Equation (4) is typically used in the literature to constrain the dust temperature, which can then be used to measure the dust mass.However, as noted in Casey (2012) and demonstrated in Cochrane et al. (2022), the dust temperature has a profound effect on the inferred dust mass when not in the Rayleigh-Jeans regime of the FIR SED: since the Planck function is now proportional to - h kT 1 MBB , a 5 K difference in measured T MBB (from 30 to 35 K) results in a ∼125% increase in the inferred dust mass.

Characterizing the Cosmic Sands FIR SEDs
In Figure 8, we show the process of inferring the Cosmic Sands dust temperatures from modified blackbody fits to photometry sampled from the POWDERDAY SEDs.The basic idea is to treat the Cosmic Sands galaxies and their SEDs as though we are "observing" them, and test the efficacy of modified blackbody fitting on high-z dusty galaxies.First, however, to understand how well the Cosmic Sands galaxies are characterized by a modified blackbody function, we fit Equation (5) to the full POWDERDAY SEDs from 30 to 1000 μm. 10 We use the fitting routine MCIRSED introduced in Drew & Casey (2022) and allow all parameters (T MBB , β, ν 0 , and α) to vary.MCIRSED uses the PYMC3 package (Salvatier et al. 2016) to sample parameter posteriors, which are constrained by uniform priors spanning T CMB (z = 6.5) < T MBB < 250 K, 0.5 < β < 3.0, 0.5 < α < 4.0, and 30 < λ 0 < 400 μm (where CMB is cosmic microwave background).
In Figure 9, we plot T MBB as a function of the true dust temperatures shown in Figure 2: luminosity-weighted, massweighted, and the temperature corresponding to the peak wavelength of the FIR SED.The points are color coded by the median of the β posterior distribution, and the sizes of the points correspond to the dust mass of each galaxy, with larger points indicating larger dust masses.First, we note that the Cosmic Sands SEDs span a range of β values but are generally concentrated around β = 1.5 and β = 1.9, with the dust-rich  systems favoring larger β values.Recall that, as part of the forward modeling process, the dust emission from each particle is a function of the particle temperature, assuming modified blackbody emission with the emissivity in turn dependent on the Weingartner & Draine (2001) grain properties.Each particle containing dust can theoretically emit a blackbody with a unique emissivity, resulting in an effective FIR SED on galactic scales made up of various blackbody spectra at different temperatures and emissivities.The observed effective emissivity is then dependent on the emissivity of each particle and the column density of the dust, which will impact the optical depth toward the observer.What this means is that from the observed SED, we can really only probe this effective emissivity and optical depth, which are both dependent on the underlying temperature distribution.A consequence of this is that an apparently cold system with a measured emissivity of β = 2 can either be due to optically thin emission from cold dust or optically thick cold dust obscuring a temperature gradient toward warm dust (Jin et al. 2022).
Unsurprisingly, the modified blackbody temperatures are generally a poor estimate for both the luminosity-weighted and the mass-weighted dust temperatures.For dust-rich systems, the blackbody fit underestimates both temperatures because the dust column density obscures the rest of the warmer dust distribution from view, thereby appearing colder.By contrast, the galaxies with a lower dust content have blackbody temperatures that are too warm compared to their physically weighted temperatures.Because the optical depth in these systems is much smaller, we can "see" into the entire temperature distribution, but the SED will dominated by the emission from the warmest dust.Thus, while still biased warm, the blackbody temperatures are a bit closer to the luminosity-weighted temperatures than the mass-weighted.For both high and low dust mass galaxies, the blackbody temperatures are closest to the temperatures corresponding to the peak of the FIR SED as they are both measures of apparent temperature.However, we note that there is an average bias of 6.25 K, which is smaller than the offsets from the physical temperatures but relatively large compared to the dynamic range of the effective dust temperatures and the uncertainties associated with the fit.

Fitting FIR Photometry with a Modified Blackbody Model
Fitting a modified blackbody+power-law function to the Cosmic Sands SEDs generally produces reasonable SED models, but even in a "best-case scenario" in which the true shape of the FIR SED is known, the accuracy of the measured dust temperatures depends on the dust content of the galaxy, specifically the dust optical depth in the FIR.The primary challenges in inferring dust properties from SEDs are the typically sparse FIR SED coverage and the degeneracies between both the intrinsic dust properties and the chosen model parameters that are only marginally broken with shorterwavelength FIR constraints.To explore how these factors contribute to modeling dust temperatures, we test two different scenarios that differ in SED coverage using the z = 6.5 Cosmic Sands galaxy sample.

Scenario 1: Well-sampled FIR SED
First, we assume optimal SED coverage and fit the modified blackbody + MIR power-law model in Equation (5) to five FIR bands (rest-frame λ rest ∼ 46, 65, 116, 185, 267 μm, which correspond to observations at 350, 500, 870, 1400, and 2000 μm for a z = 6.5 galaxy).We again use MCIRSED to sample parameter posteriors.While the effective β of the Cosmic Sands galaxies is a bimodal distribution, the intent of this exercise is to show the impact of balancing the number of free parameters with available SED constraints.Thus, we fix β = 1.5 and α = 2.0 but allow ν 0 = c/λ 0 to vary, enabling optically thick solutions toward λ 0 ∼ 400 μm in the rest frame.The parameter posteriors are constrained by uniform priors spanning T CMB (z = 6.5) < T MBB < 250 K and 30 < λ 0 < 400 μm.We note that the model and parameter choices are similar to analysis of high-z dusty galaxies with reasonable FIR coverage that have stellar and dust properties comparable to the Cosmic Sands galaxies (e.g., Greve et al. 2012;Spilker et al. 2016;Lower et al. 2023) and can therefore serve as illustrative comparisons.
We plot the median and 16th-84th percentile width of the dust temperature posteriors for each galaxy in Figure 10, comparing to the mass-weighted temperature (left panel) and the luminosity-weighted temperature (right panel).The results from these fits are nearly identical to the baseline fits to the entire FIR SED: for galaxies with low dust opacities, the dust temperatures inferred from modified blackbody fits to a wellsampled SED overestimate the mass-weighted temperatures but are reasonable estimates for the luminosity-weighted temperatures, with an average bias of 7.3 K for galaxies with τ 100 μm < 0.5.In these galaxies, the luminosity of the FIR SED is dominated by warm dust that is probed by the modified blackbody fits, given the ability to properly model the transition from optically thick to thin emission.Also similar to Figure 9, the dust temperature estimates for the most optically thick galaxies tend to underestimate both the mass-weighted and the luminosity-weighted temperatures, as the bulk of the warm-hot dust mass is still highly obscured by colder dust.From this, we note the importance of allowing for general opacity solutions but also caution the utility of single temperature fits to determine the underlying dust properties, as even with a relatively well-sampled FIR SED the apparent temperatures do not necessarily align with physically meaningful temperatures.

Scenario 2: Sparsely Sampled FIR SED
Imitating large FIR surveys of high-z galaxies that typically have sparse wavelength converge, we fit the modified blackbody model in Equation (3) to two FIR bands (rest-frame λ = 100, 160 μm) that allow us to neglect the power-law component that would affect data blueward of ∼50 μm.As in the first scenario, we fix β = 1.5 and assume that emission approaches optically thick at λ 0 = 100 μm.We use the EMCEE package (Foreman-Mackey et al. 2013) to sample the dust temperature posteriors, which are constrained by a uniform prior spanning T CMB (z = 6.5) < T MBB < 250 K. Here, our mock data and model choices are similar to analysis from Algera et al. (2023), and can again serve as illustrative comparisons for dust temperature measurements from sparse data.
We show the modified blackbody dust temperatures as a function of mass-weighted and luminosity-weighted temperatures in the left and right panels of Figure 11, where we plot the median and 16th-84th percentile width of the T MBB posterior for each galaxy.Though this model accounts for FIR dust optical depth, the wavelength at which dust transitions from optically thick to thin is fixed at 100 μm.First, the uncertainties on the dust temperatures inferred from these fits are much larger than in the previous model, with an average posterior posterior width of 35.5 K, similar to the uncertainties in dust temperature measurements in Algera et al. (2023) across their three model variants.This is mainly driven by the fact that only two photometric bands are being used to constrain two model parameters (amplitude and temperature), and that, while the shortest wavelength constraint at 100 μm is nearer to the peak of emission, it is still not a strong enough lever arm to break the degeneracies between the wide parameter space that fits the Rayleigh-Jeans tail of these SEDs.
If we focus on the median inferred temperatures, these are marginal overestimates of the mass-weighted temperatures of less dusty galaxies and marginal underestimates of the luminosity-weighted temperatures.The dust-rich galaxies, by contrast, all have similar temperature estimates (40-50 K).As is clear from both exercises, the use of single temperature models with fixed opacity and spectral emissivity parameters severely hampers the measurement of dust temperatures by imposing strong biases and implicit uncertainties.We note that such a methodology is common in the literature, especially at high redshift (e.g., Strandet et al. 2016;Faisst et al. 2020;Algera et al. 2023), where the dearth of FIR data necessitates simplifying assumptions and can bias any conclusions inferred from these models.

Discussion
In this analysis, we have demonstrated challenges in measuring dust temperatures for a sample of z ∼ 6 galaxies generated from a high-resolution cosmological simulation.These galaxies are an ideal testing ground for dust temperature techniques, as they span a wide range of stellar, dust, and radiative properties at the earliest epochs.Below, we discuss our results in the context of previous dust temperature predictions from hydrodynamical simulations, dust detections of high-z galaxies, and the potential caveats owing to assumptions made when forward modeling our simulated galaxies.

Comparisons to High-z Dust Temperatures from
Hydrodynamical Models in the Literature The dust temperature modeling we present above demonstrates the array of challenges in characterizing the dust properties of high-z galaxies from both a forward modeling perspective and the backwards modeling perspective.For instance, while the Cosmic Sands galaxies benefit from having a self-consistent dust evolution model within SIMBA, we do need to make assumptions related to the optical properties of the dust grains within POWDERDAY, the implications of which are described in further detail in Section 5.4.And while our intent is to highlight the difficulties in modeling dust properties rather than predict the dust temperatures of high-z galaxies in comparison to either observational studies or predictions from other hydrodynamical models, it would be useful to discuss the modeling choices made here in comparison to other simulations.
In recent years, several cosmological hydrodynamical galaxy formation models have studied the impact of dust on the evolution of galaxies, including galaxy chemical enrichment and dust mass evolution (e.g., McKinnon et al. 2016;Graziani et al. 2020;Choban et al. 2022;Di Cesare et al. 2023;Lewis et al. 2023) and the effect of dust on the radiative properties of galaxies (e.g., Liang et al. 2019;Mushtaq & Puttasiddappa 2022;Shen et al. 2022;Vijayan et al. 2022).The former analysis necessitates the implementation of a dust evolution model within the galaxy formation framework (akin to SIMBA), while the latter can be studied using post-process radiative transfer models such as SKIRT (Camps & Baes 2015, 2020) and POWDERDAY (Narayanan et al. 2021).Below we focus specifically on studies that presented dust temperature predictions for high-z galaxies.
In Liang et al. (2019), the authors present dust temperature predictions from the MASSIVEFIRE suite of simulations across a range of redshifts.The dust emission properties were modeled with SKIRT with an adopted fiducial dust-to-metals (DTM) ratio of 0.4 and includes the impact of both diffuse dust and dust in stellar birth clouds, implemented with the MAPPINGSIII model.This is in contrast to Cosmic Sands in which we neglect subresolution radiative transfer processes.Focusing on the galaxies at z = 6 (with stellar masses ranging from 10 9 < M star < 10 11 M e and dust masses ranging from 10 6.5 < M dust < 10 8 M e ), they predict mass-weighted dust temperatures ranging from ∼20 to 40 K in contrast to the temperatures from Cosmic Sands that range from 30 to 80 K with one galaxy at ∼120 K. Notably, however, the modified blackbody dust emission models the authors fit to mock photometry sampled from the SKIRT SEDs struggled to reproduce the shape of the true SEDs, resulting in a large scatter in dust temperatures inferred from these fits, similar to the results for the Cosmic Sands galaxies shown in Section 4. And though the predicted range of z > 6 dust temperatures are different between the two simulations, they do agree that temperatures derived from modified blackbody SED fits tend to overestimate the temperature corresponding to the peak of the FIR SED.
While the IR luminosities and dust masses of the FIRE galaxies are comparable to those of the Cosmic Sands galaxies, we note that the predicted (and directly modeled) global DTMs ratio in the Cosmic Sands galaxies range from ∼0.05 to 0.6 with significant variations within each individual galaxy, in contrast to the value of 0.4 Liang et al. (2019) assumed.Furthermore, the use of the MAPPINGSIII model to account for dust emission in subresolution star-forming regions can have significant impacts on the predicted dust temperatures, as we describe in Section 5.4.Therefore, we conclude that differences in the post-processed treatment of dust emission and the treatment of the star-dust geometry/dust abundance drive the differences in the mass-weighted temperature distributions between the two models.Shen et al. (2022) presented high-z dust temperatures of massive (10 10 < M star < 10 11 M e ), IR luminous galaxies from the IllustrisTNG simulations.Dust absorption and emission are modeled with SKIRT, assuming the redshift-dependent DTM mass ratio presented in Vogelsberger et al. (2020), which is applied to cold, star-forming gas cells.Dust temperatures inferred from the peak of the FIR SKIRT SEDs range between ∼45 and 60 K for galaxies at z = 6-8, increasing to ∼70-90 K if measured from modified blackbody fits to the full FIR SED.While the authors describe both the peak temperature and the dust temperature inferred from modified blackbody fits as luminosity-weighted temperatures, it is important to note the caveats related to optical depths discussed in Section 3.1.Compared to the Cosmic Sands galaxies, the TNG predicted dust temperatures derived from SED fits are comparable to the peak and luminosity-weighted temperatures we present in Section 3, perhaps owing to the more nuanced treatment of the dust abundance distribution in relation to star-forming regions.
Finally, in Vijayan et al. (2022), the authors present dust properties including temperatures for galaxies from the FLARES simulation across a stellar mass range of 10 8 < M star < 10 11 M e .As above, dust absorption/emission is modeled with SKIRT, with dust abundances scaled to the metal mass in cool (T < 10 6 K) gas cells following the DTM mass relation of Vijayan et al. (2019), which relates the gas-phase ISM metallicity and the mass-weighted stellar population age to the dust abundance in each galaxy.For galaxies z = 6-10, the DTMs ratio ranges from 0.01 to 0.2, which is much lower than both the assumed DTM of Liang et al. (2019) and the DTM values found in the Cosmic Sands galaxies but is closer to the redshift-dependent DTM value used in Shen et al. (2022; DTM ∼ 0.1 at z = 6.5).Also similar to Shen et al. (2022), the dust temperatures corresponding to the peak of the FIR SED are ∼10-25 K lower than those derived from modified blackbody fits to the entire SED.The latter dust temperatures are comparable to the luminosity-weighted temperatures of the Cosmic Sands galaxies and to those inferred from modified blackbody fits to photometry sampled from our POWDER- DAY SEDs.

Implications for High-z Dust Temperature Measurements
As we enter the joint ALMA-JWST era, evidence of greaterthan-expected dust enrichment at the earliest epochs is continuing to push our understanding of galaxy formation (e.g., Sommovigo et al. 2022a;Akins et al. 2023;Labbé et al. 2023).The challenges in modeling dust temperatures for dusty, high-z galaxies outlined in this analysis serve to underscore the findings from previous studies spanning decades for both the Milky Way and nearby galaxies (Hildebrand 1983;Dale & Helou 2002;Draine 2006;Casey 2012;Kirkpatrick et al. 2014;Lee et al. 2015) to galaxies at z = 2 and beyond (Cortzen et al. 2020;Jin et al. 2022;Sommovigo et al. 2022a) that sought to characterize the dust content in various systems despite the onslaught of systemic model degeneracies.These degeneracies, also explored by previous studies primarily focused on low-redshift systems (e.g., Klaas et al. 2001;Shetty et al. 2009aShetty et al. , 2009b;;Rangwala et al. 2011;Kelly et al. 2012), are amplified at high redshift, where the additional dearth of FIR data makes the characterization of dust problematic.The assumptions that we make when fitting FIR SEDs at low redshifts may not hold for high redshifts as shown in Section 4, adding implicit uncertainties to already uncertain measurements.
Galaxies at these epochs that have relatively well sampled FIR SEDs have primarily been FIR-or submillimeter-selected dusty star-forming galaxies.For instance, the majority of z > 5 South Pole Telescope (SPT) submillimeter galaxies (SMGs) have detections from Herschel/SPIRE, APEX/LABOCA, SPT, and ALMA, effectively covering rest-frame 40-500 μm.
In Reuter et al. (2020), the authors constrain the dust properties for these galaxies by fitting a modified blackbody function to 4 −6 flux bands, depending on the galaxy, and account for optically thick emission in the FIR, similar to our Scenario 1.Like the estimates for the Cosmic Sands galaxies, the inferred blackbody dust temperatures for the SPT SMG galaxies range between 30 and 70 K, with a characteristic uncertainty of ∼16 K. Similar studies suggest the typical dust temperatures found in SMGs ∼30-40 K (e.g., Greve et al. 2012;Swinbank et al. 2014;Strandet et al. 2016;Danielson et al. 2017;Sun et al. 2022).For UV-selected galaxies targeted by ALMA, studies tend to find higher dust temperatures ranging from ∼40 to 100 K (e.g., Bakx et al. 2020;Faisst et al. 2020;Inami et al. 2022;Sommovigo et al. 2022a;Algera et al. 2023).
As shown above, for the relatively dust poor Cosmic Sands galaxies, the inferred temperatures generally align with the luminosity-weighted dust temperatures, while for the dust-rich galaxies the inferred temperatures align with the massweighted.The lower apparent temperatures in high-z SMGs may also be the result of high dust opacities that obscure the bulk of the dust temperature distribution (Jin et al. 2022).We stress the importance of allowing general opacity solutions when inferring dust properties of high-z dusty galaxies.The relation between inferred and "true" dust temperatures depends on both the shape of the SED (which folds in the dust temperature distribution, the star-dust geometry, and the dust optical properties), and the aforementioned model assumptions.As noted in Section 3, there is typically not a 1:1 nor even a linear relation between the apparent dust temperature and the mass-weighted or luminosity-weighted dust temperature of a galaxy.This is supported by the analyses of Liang et al. (2019) and Vijayan et al. (2022) from the FIRE and FLARES simulations, respectively, when relating apparent dust temperatures to the mass-weighted temperatures.
Inferring the dust properties of UV-or optically selected galaxies at z > 6 is potentially more challenging, as significant observation time must be dedicated to detecting faint dust emission.Large ALMA programs such as ALPINE (Le Fèvre et al. 2020) and REBELs (Bouwens et al. 2022) may only have a single dust continuum detection in the FIR (ALMA bands 6 or 7), which can lead to significant uncertainties when inferring dust temperatures as modified blackbody functions are not well constrained by the Rayleigh-Jeans tail.To constrain degenerate solutions, Sommovigo et al. (2021) proposed using [C II] detections to independently constrain the galaxy dust masses, thereby breaking the degeneracy between dust mass and temperature.In Sommovigo et al. (2022aSommovigo et al. ( , 2022b)), the authors applied this methodology to high-z galaxies from the ALPINE and REBELs surveys, inferring dust temperatures by combining the rest-frame 158 μm [C II] and dust continuum detections, assuming a modified blackbody function and optically thin emission.The dust temperatures inferred from this method range from 30 to 60 K with a median temperature for the REBELs galaxies of 47 K (Sommovigo et al. 2022b).Faisst et al. (2020) fit the ALMA-detected FIR SEDs of four Keck/DEIMOS-selected z = 5-6 galaxies in the COSMOS field, and are careful to point out that the dust temperatures measured from the SEDs are more similar to luminosityweighted temperatures and may miss the bulk of dust mass in colder regions.Nevertheless, the galaxy SEDs are assumed to be optically thin past 100 μm, but constraining both the emissivity index and the dust temperature with rest-frame photometry at ∼100-200 μm is challenging and significant uncertainties persist.The authors predict dust temperatures of 30−40 K for these galaxies, matching the T d (z) relation predicted from hydrodynamical simulations presented in Liang et al. (2019).These measurements also agree with our result from Scenario 2, in which we infer dust temperatures between 30 and 60 K, albeit with relatively large uncertainties.Importantly, however, those temperatures do not correlate significantly with either the true mass-weighted or luminosityweighted temperatures and so we caution the use of singleband dust continuum measurements to characterize the dust properties in observed galaxies.Fudamoto et al. (2023) presented new estimates for high-z galaxies derived from ALMA observations, including A1689-zD1 (Watson et al. 2015;Knudsen et al. 2017;Bakx et al. 2021), B14-65666 (Hashimoto et al. 2019(Hashimoto et al. , 2023)), and MACS0416-Y1 (Tamura et al. 2019;Bakx et al. 2020).Their dust emission model, based on work by Inoue et al. (2020), accounts for clumpiness in the ISM and enables the derivation of dust temperatures and luminosities from single-band dust continuum measurements.With this method, the authors found their temperature measurements generally agreed with existing literature estimates of dust temperatures (e.g., Bakx et al. 2020;Sommovigo et al. 2021) but with relatively significant uncertainties from the dust/ISM geometry.Witstok et al. (2023) compiled SEDs and measurements of dust emission areas of 17 high-redshift galaxies, including SMGs, "normal" star-forming galaxies, and QSOs, to selfconsistently model the dust properties with modified blackbody fitting assuming a general opacity model, linking independently derived dust mass surface density measurements to the FIR emission.For galaxies at z = 6-8 with stellar masses around 10 10 M e , they measure dust temperatures and IR luminosities similar to ULIRG systems at z = 0, ranging from 35 to 80 K, and find most of the galaxies in their sample are fit with an emissivity index of β ∼ 2, similar to the dustier systems in the Cosmic Sands sample fit in Section 4.2.The authors also discuss in detail the nontrivial nature of deriving dust properties of high-z galaxies, highlighting the complex relationship between measures of apparent temperature (T peak , T MBB ) and physical temperatures (T mass-weighted , T luminosity-weighted ).
Finally, Algera et al. (2023) presented new dust temperature measurements for three previously studied REBELs galaxies for which additional ALMA band 8 detections were obtained at rest-frame ∼90 μm.With the leverage of a shorter-wavelength detection, the authors fit modified blackbody functions with three techniques: β = 1.5, β = 2.0, and a variable β for both optically thin and thick assumptions, in which the dust is assumed to transition from optically thick to thin at λ c = 100 μm (ν c = 3000 GHz).The temperatures inferred from these fits (assuming optically thin emission with β = 1.5) are similar to the temperatures inferred in Sommovigo et al. (2022b) but with larger uncertainties.The temperatures inferred from the optically thick model are 10−15 K larger than those from Sommovigo et al. (2022b), implying that either the assumption of optically thin emission at λ > 40 μm is not valid for the REBELs galaxies or the [C II]-derived dust mass constraint does not hold when an additional photometry constraints are available, or both (see, e.g., Vizgan et al. 2022).
What is clear is that the technique of inferring dust temperatures via modified blackbody fits to FIR SEDs is nontrivial for even well-sampled SEDs, driven by degeneracies between dust properties and the fact that FIR emission may be optically thick out to several hundred microns, rendering a bulk of the dust mass invisible.The application of such models to high-z galaxies, where dust abundances and densities, star-dust geometries, and ISM conditions may be significantly different than in low-z galaxies, while common, can introduce biases in measurements of dust properties as demonstrated in Section 4. The relation between luminosity-weighted temperatures and mass-weighted temperatures depends on the star formation and dust properties of galaxies, making it challenging to apply a "one-temperature-fits-all" model to a large sample of galaxies.The various techniques tested in this analysis demonstrate the large uncertainties on individual measurements, even when the FIR SED is relatively well constrained.

Implications for High-z Dust Mass Measurements
While our analysis has focused solely on the techniques used to infer galaxy dust temperatures from FIR continuum measurements, we briefly discuss the implications of these measurements, and their uncertainties, on their application to inferring dust masses.Complications arise when attempting to infer other dust properties from dust temperatures measured from the FIR SED; luminosity-weighted dust temperatures typically do not trace the mass-weighted dust temperature, which is what is needed to estimate other dust properties such as mass.As pointed out by Casey (2012) and Cochrane et al. (2022), relatively small uncertainties on dust temperature result in significant errors on dust mass, especially when measured from photometry on the Rayleigh-Jeans tail that do little to constrain modified blackbody emission at shorter wavelengths.
If we use the dust temperature estimates from Scenario 2 to infer dust masses, again assuming modified blackbody emission, the median offset is −0.25 dex-nearly a factor of 2-from the true dust masses.The uncertainties propagated from the inferred dust temperature posteriors result in an average fractional error of 122%, significantly impacting our ability to conclusively determine the true dust masses of these galaxies.We highlight that this dust mass offset should be interpreted as a lower limit prediction for the actual offset between inferred and true dust masses for observed galaxies, as we are matching the dust optical properties (κ, β) between the backwards and forwards modeling of the dust emission.As explored in Inoue et al. (2020) and Fanciullo et al. (2020), the adopted dust opacity and ISM geometry models have a significant impact on the derived dust masses even if the input dust temperature is correct, leading to inconsistencies in mass measurements between methodologies (Jin et al. 2022).

Caveats to Our Model
A major limiting factor in the predictive power of the Cosmic Sands galaxies is the reliance on subresolution and post-process modeling of dust emission properties, including the treatment of dust in stellar birth clouds and the lack of onthe-fly ionized radiative transfer in H II and photodissociation regions.In this work, we generate UV-FIR SEDs with postprocess radiative transfer but neglect additional attenuation terms toward young stars in birth clouds; i.e., we model only the diffuse dust attenuation and emission.The inclusion of a subresolution dust model (e.g., Charlot & Fall 2000, in which young stars are attenuated by an additional birth cloud dust component), can significantly change the shape of the FIR SED as well as the dust temperature distribution.
For instance, placing a dust screen with τ V = 0.5 mag in front of stars <10 Myr old increases the total IR luminosity by a factor of 1.3-2 on average across the Cosmic Sands galaxies and can shift the dust temperature corresponding to the peak of the FIR SED down by nearly 30 K as less UV light reaches the diffuse dust.We highlight two galaxies in Figure 12 to demonstrate this scenario.The two galaxies have relatively low (top panel) and high (bottom panel) dust optical depths at 100 μm.The black curves show the fiducial SEDs and the blue curves show the impact of adding a uniform dust screen in front of young stars.For the low dust opacity galaxy, the addition of birth clouds significantly changes the shape of the FIR SED, increases the IR luminosity by a factor of 1.6, and decreases the peak temperature by 25 K.By contrast, the high dust opacity galaxy experiences a modest increase in IR luminosity with the peak of the FIR SED shifting slightly to higher wavelengths.In this regard, the dust temperatures presented here may not be directly comparable to dust temperature distributions of real galaxies, or to the trends found in other works analyzing the dust temperatures of high-z galaxies from hydrodynamical simulations (e.g., Liang et al. 2019).However, these model assumptions will not significantly impact the comparisons between the characteristic dust temperatures (luminosityand mass-weighted) and the temperatures inferred from modified blackbody fits in our analysis as both are related to the same post-process radiative transfer assumptions.
A perhaps more fundamental uncertainty is the underlying galaxy formation physics in SIMBA, namely the treatment of stellar sources of feedback on the ISM.A significant obstacle in numerically modeling galaxy evolution is the limitation of resolution such that star formation will not be self-regulating at the smallest scales.To get around this, SIMBA adopts an artificial equation of state that sets a pressure floor that suppresses fragmentation below the scale of the smoothing volume (Davé et al. 2016(Davé et al. , 2019)).This is in contrast to models such as FIRE (Hopkins et al. 2014) and SMUGGLE (Marinacci et al. 2019), in which star formation and feedback from stellar sources is explicitly resolved, allowing the multiphase medium to be more accurately modeled.
In the context of this analysis, the nonlinear coupling of stellar feedback modes on the surrounding medium can play a significant role in the predicted dust evolution, star-dust geometry, and dust temperature distribution of high-z galaxies (Esmerian & Gnedin 2023).For instance, in the SIMBA model, the mass-loading factor of the star formation driven galactic winds is artificially reduced at z > 3 to better match observations of galaxy growth at high redshift (Davé et al. 2019).This could lead to the pileup of gas in the nuclear regions of high-z galaxies, yielding large surface densities and dust abundances that could otherwise be dispelled by stellar winds.The consequences of these fundamental, theoretical uncertainties clearly add to the overall uncertainty budget when predicting dust temperatures from hydrodynamical models.
Lastly, we are unable to predict the impact of variable dust extinction on the FIR emission of high-z galaxies.Specifically, we do not model an evolving grain-size distribution in SIMBA, which could significantly change the thermodynamics of the ISM in these galaxies and drive changes in the ways dust absorbs and emits as a function of wavelength.It is known that dust of differing grain sizes have different optical properties and thus dominate galaxy spectra at different wavelengths; e.g., PAHs with grain sizes ∼1-10 Å dominate galaxy emission in the near-to mid-IR (e.g., Narayanan et al. 2023).While we cannot test the impact of galaxy dust grain size distributions and how this manifests in IR SEDs, we predict that the addition of a self-consistent grain size evolution model, resulting in nonuniversal dust optical properties and perhaps a wider diversity in FIR SED shapes, would result in greater uncertainties in the inferred dust temperatures from modified blackbody fits to FIR SEDs.

Outlook and Paths Forward
The question remains of how can we go forward from this analysis, given the nontrivial nature of measuring dust temperature and the limited relation this temperature has with the real, physical distribution of dust temperatures within a galaxy.If one's goal is simply to estimate some representative dust temperature, the results from Section 4.3 demonstrate the necessity of robust photometric coverage of the FIR SED.Even then, the utility of this temperature measurement is unclear; the relation between it and the physical temperatures describing the underlying distribution are rarely 1:1, rendering any measure of, e.g., dust mass or the thermal state of the ISM, highly uncertain.
If instead the goal of measuring dust properties is to gain a better understanding of the galaxy as a whole, it would be more powerful to place FIR continuum observations in the context of the entire SED, including nebular and molecular emission lines where available, to jointly map UV-FIR photometry to the underlying physical properties.The constraints on dust and ISM properties from UV and optical photometry (e.g., dust composition, grain size distribution, metallicity, and star-dust geometry) are complemented by the constraints in the FIR (gas mass, and dust mass and temperature distribution), allowing one to estimate galaxy physical properties simultaneously and with respect to the correlations that physically exist between galaxy components (e.g., the dust emission in the FIR will be influenced by the model SFH and UV-optical attenuation).Balancing the energy transfer between starlight and dust provides more information than either measurement gives individually.
In other words, the exercises conducted here have limited utility in terms of offering solutions for improving dust temperature measurements.This said, studies such as Haskell et al. (2023), Jones et al. (2022), andJones &Stanway (2023) highlight the benefits of a panchromatic perspective, enforcing symmetry between dust properties modeled in both the UVoptical regime and in the FIR.Jones et al. (2022) fit IRluminous galaxies from COSMOS (0 < z < 0.5), with Jones & Stanway (2023) extending to higher redshifts (z < 6), and find that the inferred dust and stellar properties are mutually dependent on the assumed SFH, stellar population synthesis model, and dust emission model.In testing the utility of the energy balance constraint on high-z galaxies simulated with FIRE, Haskell et al. (2023) demonstrate improvements in inferring dust mass and IR luminosity from panchromatic photometry compared to estimates from FIR photometry alone.
Of course, there are significant caveats to this methodology.For instance, as the model parameter space increases, unless robustly constrained by high-fidelity spectro-photometric observations, uncertainties will inflate to reflect the degree to which degenerate solutions exist.Moreover, the current survey of high-z galaxies is dominated by galaxies selected based on rest-frame UV-optical or FIR emission, and thus many systems are missing complementary photometry in other bands.In this case, relying on bright emission lines in the FIR, such as CO or [C II], can be beneficial when trying to characterize the Figure 12.SEDs with (blue curves) and without (black curves) subresolution dust absorption and emission for two Cosmic Sands galaxies chosen to roughly span the range of dust optical depths at 100 μm.The galaxy in the top panel has a relatively low dust optical depth, and the addition of a birth cloud model significantly changes the shape of the FIR SED, increases the IR luminosity by a factor of 1.6, and decreases the peak temperature by 25 K.The bottom panel shows a galaxy with higher dust optical depth, and the changes to the IR luminosity and peak temperature by adding birth cloud dust are less severe but still meaningful.thermodynamics of the ISM (e.g., Reuter et al. 2023).And though many studies have been devoted to robustly modeling stellar emission with sophisticated implementations used in codes such as PROSPECTOR (Johnson et al. 2021) and BAGPIPES (Carnall et al. 2018), the same is only marginally true for dust emission, which is primarily based on course-grid templates such as Draine & Li (2007) and potentially ineffective for high-z galaxies.
Lastly, hydrodynamical simulations can be used to predict the mapping between physical properties and observables.For instance, Cochrane et al. (2022) model the FIR emission from galaxies from the FIRE simulation suite and find tight relations between certain flux density ratios and the mass-weighted temperatures of both star-forming and quenched galaxies, though this is highly dependent on the underlying physics of the FIRE galaxy formation model and the details of the forward modeling of the dust emission.Similarly, techniques instead using machine-learning algorithms trained on simulations to find the nonlinear mappings between luminosity and physical properties also show promise: indeed, Gilda et al. (2021) demonstrate with MIRKWOOD significant improvements in inferring galaxy physical properties from broadband photometry compared to traditional parametric SED modeling across a range of galaxy types and masses.This said, similar to the Cochrane et al. (2022) results, the Gilda et al. (2021) machinelearning model is dependent on the simulations and forward modeling methods used to train the machine-learning algorithm.

Conclusion
Using galaxies from the Cosmic Sands suite of massive, dusty galaxies in the EoR, we have demonstrated the challenges in modeling dust temperatures at high redshift, both forwards from hydrodynamical simulations post-processed with 3D dust radiative transfer and backwards from modified blackbody fits to broadband FIR photometry.We first examined the "true" dust temperatures of the Cosmic Sands galaxies, calculated from POWDERDAYʼs radiative transfer coupled with the SIMBA dust evolution model.Intuitively, the dust temperature distributions are diverse from galaxy to galaxy, but also within each galaxy as shown in Figure 2, with mass-and luminosity-weighted temperatures producing unique distributions for each galaxy.As shown in Figures 4 and 5, the star-dust geometry and ISM column densities contribute to this diversity by breaking correlations between intrinsic and apparent dust properties for dustier, star-forming galaxies.It also makes comparisons between theoretical predictions and observations of dust properties challenging as the mapping between intrinsic temperature and observable temperature depends on each galaxy's properties such as mass, star-dust geometry, and age.Lastly, in an attempt to characterize the full dust FIR SEDs with a modified blackbody model, we found that the Cosmic Sands FIR dust opacities are not well described by the commonly assumed power-law relation with emission frequency, t n n = b ( ) 0 .Then, starting from the POWDERDAY FIR SEDs, we attempted to infer dust temperatures as an observer would: we ran two cases that differed by SED coverage, spanning the range of photometry available for z = 7 galaxies, and by model assumptions regarding the dust optical depth.Our results, summarized in Figures 10 and 11, demonstrate the difficulties in constraining dust temperatures.The inherent degeneracies between dust properties (temperature, mass, star-dust geometry, dust grain size, composition, and structure) make it challenging to robustly constrain the dust temperature even for well-sampled FIR SEDs (Kovács et al. 2010;Casey 2012).As is clear from this exercise, fitting FIR SEDs with a singletemperature component and a simple parameterization for the dust optical depth does not work for a majority of high-z galaxies.Until high-z galaxy FIR SEDs can be more robustly detected over a wide wavelength range, these uncertainties pertaining to dust properties will persist.
From this analysis, we have shown that both theoretical and observational characterizations of galaxy dust properties are fraught with uncertainties at high-z.We note that predictions for dust temperatures at high redshift are highly dependent on model assumptions regarding dust evolution (e.g., dust abundances and star-dust geometry) and dust radiative transfer (e.g., subresolution dust in star-forming regions as shown in in Figure 12, and dust optical properties) and many of these factors remain relatively unconstrained or difficult to implement in a self-consistent manner.However, results from observational studies of high-z galaxies can be highly uncertain, both implicitly and explicitly, depending on the model assumptions and the ability to constrain FIR dust emission with sparse photometric data.Thus, any prediction for high-z dust properties from theoretical models should be placed in the appropriate context with regard to specific model assumptions, and caution should be taken when comparing theoretical predictions to results from observational studies.

Figure 1 .
Figure1.Schematic detailing the forward dust temperature modeling with POWDERDAY radiative transfer.Starting from Cosmic Sands galaxies evolved with the SIMBA galaxy formation framework, we generate an octree grid onto which the gas particles (containing the dust mass information) are smoothed.We then assume dust grain properties (grain size distribution and composition) corresponding to the Weingartner & Draine (2001) R V = 3.1 extinction curve.This sets the emissivities and opacities of the dust in each cell.From there, the radiative transfer propagates through the dusty ISM, producing UV-FIR SEDs and dust temperature distributions for each galaxy.

Figure 2 .
Figure 2. Dust temperature distributions for a selected sample of the Cosmic Sands galaxies, generated via POWDERDAY radiative transfer.Galaxies increase in dust mass from the top-left to bottom-right subplots.The mass-weighted (luminosity-weighted) dust temperature distributions are shown in blue (red).The whole-galaxy mass-weighted and luminosity-weighted temperatures are shown by the blue and red dotted lines.The black solid lines show the temperatures corresponding to the peak of the FIR SED.

Figure 3 .
Figure 3. Luminosity-weighted dust temperatures as a function of total stellar mass, dust mass, and SFR surface density at z = 6.5.

Figure 4 .
Figure 4. Luminosity-weighted dust temperatures as a function of mass-weighted temperature (left) and temperature corresponding to the peak of the SED (right).Points are color coded by galaxy dust mass.

Figure 5 .
Figure 5. Infrared luminosity, calculated from POWDERDAY SEDs chosen at a random viewing angle for each galaxy, as a function of luminosity-weighted dust temperatures.Points are color coded by galaxy dust mass in the left panel and by dust optical depth at 100 μm in the right panel.

Figure 6 .
Figure 6.Top: dust optical depth maps at 100 μm for a galaxy at z = 4.5 viewed face-on (left) and edge-on (right).The pixel scale in both images is 5 pc.Bottom: dust optical depth maps for a different galaxy at z = 6.5 (left) and z = 5.5 (right).

Figure 7 .
Figure 7. Infrared luminosity (left) and temperature corresponding to the peak for the FIR SED (right) as a function of luminosity-weighted dust temperature.Points are color coded by galaxy dust mass.Each Cosmic Sands galaxy is represented by 25 points, corresponding to 25 SEDs generated at different viewing angles, showing the dependence of IR luminosity and apparent temperature on the orientation angle of the galaxy.The maximum spread in luminosity (0.33 dex) and temperature (8 K) are shown in each panel.

Figure 8 .
Figure8.Schematic detailing the backwards dust temperature modeling with modified blackbody fitting.We first sample photometry from the POWDERDAY SED for each galaxy (five bands for Scenario 1 and two bands for Scenario 2).For Scenario 1, we fit a modified blackbody function (MBB) with a power-law component on the Wien's side of the FIR SED.For Scenario 2, we fit just a modified blackbody function since the two photometric constraints will not constrain the power law on the Wien's side.The parameter posteriors are sampled with PYMC3 or EMCEE from which we infer the dust temperature for each galaxy.

Figure 9 .
Figure9.Dust temperature inferred from modified blackbody+power-law fits to the Cosmic Sands SEDs as a function of luminosity-weighted temperature (left), massweighted temperature (middle), and temperature corresponding to the peak of the FIR SED (right).Points are color coded by the median of the emissivity spectral index (β) posterior from the modified blackbody+power-law fit.The sizes of the points reflect their dust mass, where galaxies with larger dust masses have larger markers.

Figure 10 .
Figure 10.Dust temperatures inferred from modified blackbody fits described in Scenario 1 as a function of mass-weighted (left) and luminosity-weighted (right) temperature.Points are color coded by FIR optical depth at 100 μm as calculated in Section 3.1.

Figure 11 .
Figure 11.Dust temperatures inferred from modified blackbody fits described in Scenario 2 as a function of mass-weighted (left) and luminosity-weighted (right) temperatures.Points are color coded by FIR optical depth at 100 μm as calculated in Section 3.1.