The early Universe was dust-rich and extremely hot

We investigate the dust properties and star-formation signature of galaxies in the early universe by stacking 111227 objects in the recently released COSMOS catalogue on maps at wavelengths bracketing the peak of warmed dust emission. We find an elevated far-infrared luminosity density to redshift 8, indicating abundant dust in the early universe. We further find an increase of dust temperature with redshift, reaching 100 +- 12 K at z ~ 7, suggesting either the presence of silicate rich dust originating from Population II stars, or sources of heating beyond simply young hot stars. Lastly, we try to understand how these objects have been missed in previous surveys, and how to design observations to target them. All code, links to the data, and instructions to reproduce this research in full are located at https://github.com/marcoviero/simstack3/.


INTRODUCTION
The history of star formation at high redshift has always had a dust problem. Star formation is well traced by ultraviolet emission, but a large fraction of ultraviolet light is obscured by dust, thus proving it to be an incomplete estimator. But dust also emits light, re-radiating the ultraviolet energy it absorbed as roughly a blackbody at farinfrared (FIR) wavelengths. To measure the blackbody is to measure star formation; to measure it across cosmic time is to characterise cosmic star-formation history, and maybe learn something about the evolution of the different galaxy populations.
Fundamental questions about the early Universe still remain, such as: was it dust rich or dust poor? How much does obscured star formation contribute to the star-formation rate density? Do evolutionary trends of galaxies with redshift seen locally (e.g. temperature, specific star-formation rate) continue indefinitely, or do they plateau, or even turn over? Just how active are AGNs in this period? And how much ionizing radiation is escaping the dark curtain of dust exactly?
The catch is that measuring the far-infrared spectrum of a galaxy at high redshift is challenging. A modified blackbody spectrum (or spectral-energy distribution; SED) with a temperature of 18 K peaks in the infrared at a rest-frame of roughly 160 m, or between 400-600 m in the observer frame if emitted at the peak of star formation, ∼ 2-3. To fully bracket that peak requires combining observations from space telescopes like Spitzer and Herschel, and ground based observatories like the JCMT, LMT, and ALMA. FIR space telescopes are wonderful except for one drawback: source confusion. The primary beams are large enough to contain tens (or hundreds) of galaxies, complicating optical counterpart identification. Groundbased FIR telescopes are wonderful except for a different drawback: the atmosphere, which is to infrared light what dust is to optical light.
Until we are able to put an actively-cooled, 50 m mirror in space, ★ E-mail: marco.viero@caltech.edu the solution is stacking. Galaxies with similar physical properties estimated from optical photometry (e.g. redshift, stellar-mass, starforming activity) have, on average, similar infrared properties (e.g. infrared luminosity, dust temperature; Schreiber et al. 2018). Viero et al. (2013) developed an algorithm to stack galaxies, grouped by their similarity, simultaneously (see also Kurczynski & Gawiser 2010), thus providing a way to overcome the bias inherent to stacking in highly confused images.
The techniques in this paper are not new. The reason we are only now able to address these questions is because of newly available catalogues released to the public from COSMOS (Scoville et al. 2007). We combine the COSMOS2020 catalogue (Weaver et al. 2022) with infrared/submillimetre maps publicly available through the Herschel Extragalactic Legacy Project (HELP) and SCUBA-2 Cosmology Legacy Survey (S2CLS), and a newly released, 3-compatible version of , to measure the far-infrared SEDs from redshift 0 to 10 and stellar-masses of log( /M ) = 9.5 to 12.
We find that the early Universe was dust-rich, and unusually hot. Effective dust temperatures increase with redshift, following established relationships to = 4, but then rise more rapidly to ∼ 10. Similarly, we find that the far-infrared luminosity density ( LIR ) tracks the established star-formation rate density (SFRD or SFR ) to = 2, but then decouples from it, remaining flat to = 8. Given the unusual nature of the hot dust at > 4, we explore the mechanisms that could be responsible for this heating, how this could have been missed previously, and strategies for follow-up observations. We assume a Chabrier (2003) initial-mass function (IMF) and Planck18 cosmology (Planck Collaboration et al. 2020), with Ω M = 0.315, Ω Λ = 0.685, 0 = 67.4 km s −1 Mpc −1 , 8 = 0.811, and = 0.054. All code and instructions to reproduce these results can be found at https://github.com/marcoviero/simstack3/ and at the DOI: https://doi.org/10.5281/zenodo.6792395.

COSMOS2020 Catalogue
Our analysis is performed on data in the COSMOS 1 field (Scoville et al. 2007), centred at 10 h 00 m 26 s , 2 • 13 00 , and covering 1.606 deg 2 after masking. We use the recently released COSMOS2020 catalogue (Weaver et al. 2022), whose = 25.2 (AB)-selected sample is an improvement of nearly two magnitudes over previous releases used for stacking (Muzzin et al. 2013). This extra depth enables access to objects reaching redshift of 10, and completeness improvements of ∼ 2 dex in stellar mass.
The COSMOS2020 release comes with two options for photometry; we use FARMER/L P , generally understood to be superior for faint sources at high redshift. Following Weaver et al. (2022), we split the population into star-forming and quiescent using the − versus − selection (Ilbert et al. 2013).

Maps
Maps are selected to bracket the rest-frame dust SED at = 0-10. Unless otherwise noted, they were obtained from the Herschel Extragalactic Legacy Project (HELP 2 ). The Spitzer/MIPS map at 24 m is a Spitzer Enhanced Imaging Product (SEIP). It has an rms of 27.9 Jy, point-spread function (PSF) full-width at half maximum (FWHM) of 6.32 arcsec, and aperture/colour corrections equaling 1.24 (Béthermin et al. 2010). Maps are converted from MJy sr −1 to Jy beam −1 by dividing the map by the beam solid angle, 1.55 × 10 −9 sr.
Herschel/PACS maps at 100 and 160 m are from the PACS Evolutionary Probe (PEP; Poglitsch et al. 2010). They have rms values of 0.10 and 2.11 mJy, and FWHM of 7.5 and 11.3 arcsec. PACS maps are also converted from MJy sr −1 to Jy beam −1 , with beam solid angles of 2.03 × 10 −9 and 4.66 × 10 −9 sr. Calibration uncertainty is 4%.
The 850 m map is a product of S2CLS 3 (Geach et al. 2017). Map depth is 1.6 mJy beam −1 (1 ), and FWHM 12.1 arcsec. Note, we specifically select the non-match-filtered (NMF) maps, which retain the unresolved fluctuations on which relies.

Unbiased Stacking with
Our stacking analysis is based on , an algorithm designed to overcome the biases inherent with stacking on images dominated by confusion noise. Briefly, uses the 3D positions in a catalogue to construct and fit a model to a map. The model is made up of layers representing bins that have been defined by the user (e.g. by stellar-mass, redshift, etc.), and whose best-fittings are interpreted as the average flux density of the binned objects. The software for stacking used here is the same one detailed in Viero et al. (2013) with some minor modifications.
The first difference, following Sun et al. (2018), is that instead of stacking in redshift slices, now all redshift bins are also stacked simultaneously. This change is to minimize the potential bias that may arise from objects outside the bin boundaries, or from objects correlated with interloping sources. The next difference is motivated by Duivenvoorden et al. (2020), who showed that stacked flux densities can be biased by holes in the catalogue coverage caused by extremely bright objects like stars. Their solution is the addition of a flat foreground layer, which is now an optional flag (default=True) in . The last change concerns the bootstrap estimation. Previously, the bootstrap consisted of randomly sampling the catalogue (with replacement), which is robust in most scenarios, but in this case actually biases the fluxes low. The modification satisfies the central principle of , that all correlated objects must be stacked at the same time to prevent biasing the estimate, by randomly splitting each bin in two, with a 80:20 ratio, and stacking the entire set simultaneously. The 80% bins are retained as the bootstrapped flux, and the 20% bins are examined for outliers. The drawback with this method is that the memory required to perform the stack is doubled, requiring either a more powerful computer, or splitting the bootstrap into chunks.

Spectral Energy Distribution Fitting
To obtain SEDs from stacked flux densities we fit them with a hybrid blackbody, substituting the mid-infrared Wien side with a power law, , with fixed slope of 2.0, and emissivity index fixed at 1.8 (e.g. Casey 2012). Fits then have two free parameters, amplitude and observed temperature. We sample the parameter space with for (Foreman-Mackey et al. 2013), a user-friendly, affineinvariant ensemble sampler for Markov chain Monte Carlo (MCMC). We adopt a flat prior with generous upper and lower limits for amplitude and temperature. For five bins at = 1.5-2 and = 2.0-2.5, we add a prior on temperature to prevent the solution from preferring local minima at higher temperature. Additionally, we account for PAH emission lines contaminating the 24 m flux density by artificially inflating the MIPS errors over the redshift range 0.5 < < 3. Full covariances are derived from bootstraps to account for the highly correlated nature of bands near to one another in wavelength. Flux density measurements consistent with zero are treated as upper limits, whose contribution to the log-likelihood is via a survival function (Sawicki 2012). We use 15 000 samples, discarding the first 3000.
Infrared luminosities (LIR) are estimated as the integral of from rest-frame 8 to 1000 m (Casey 2012), and naïvely converted to the star-formation rate (SFR) using the Kennicutt (1998) relation, with a 0.23 dex correction to convert to a Chabrier (2003) IMF. Error bars account for additional uncertainties in the luminosity distances arising from the photometric redshifts.
We also estimate the SFR via the rest-frame 850 m luminosity 850 . As a proxy for the interstellar medium gas mass, it is insensitive to dust temperature d with a relatively constant luminosity-mass ratio, and is also a reliable SFR estimator (Scoville et al. 2016).

Null Tests and Other Consistency Checks
To check that the signal we measure is not noise, we redo the stack but randomly shuffle source positions. We find that the measured flux densities are centred around zero, with uncertainties consistent with the rms in the maps. We further check that the measured signal is not dominated by a few bright sources by inspecting the distribution of the 20% group of the bootstrap split, finding their relative variance is proportional to their relative size.   Next, we carefully examine the impact that low-interlopers have on the stacked SED shape at > 4 using simulated catalogs/maps generated with the simulation (Béthermin et al. 2010). We use the redshift probability distributions ( ) provided by the COS-MOS team to i) draw a redshift for each object; ii) look for the closest counterpart in the catalog; iii) bin and average their flux densities; iv) and compare the best-fit SEDs. We repeat this 100 times (notebooks and downloads of the simulation results can be found in the same GitHub repository as for the main results).
Objects at higher redshift or with lower stellar masses have broader ( ), sometimes even double-peaked; provided they are reliable probabilistic estimators of the true redshifts, this method should give a reasonable representation of the outlier bias on the mean SED. We find that the overall bias in rest-frame dust temperature is small, but that the scatter is large, particularly at > 5. We apply these excess uncertainties to the errors in the dust temperatures.
Finally, we check that the cumulative flux density is consistent with the cosmic infrared background (e.g. Dole et al. 2006), and that the contributions from different redshift and stellar-mass bins agree with Viero et al. (2013).

RESULTS
The full set of SEDs for star-forming galaxies are shown in Figure 1. The best fits are in blue, with blue shaded regions outlining the 25th and 75th quartiles of the full chain of MCMC samples. Overlaid in magenta are SEDs drawn from the star-forming main sequence of Schreiber et al. (2015), as described in § 3.3, showing a clear departure from existing trends at = 0-4.
(1)  Figure 3. Star-formation rate density vs. redshift, binned by stellar mass (coloured lines), and summed (black line with 1 shaded region). Starformation rates are converted from far-infrared luminosities -derived from the integral of intensities from (rest-frame) 8 to 1000 m-by scaling with the Kennicutt (1998) relation. We find we are in broad agreement with a selection of the latest SFRD from models (Béthermin et al. 2017;Pillepich et al. 2018) and measurements (Zavala et al. 2021). Also plotted is the restframe 850 m-derived star-formation rate density as a light blue line and shaded region (1 ). The LIR-derived SFRD tracks 850 m-derived relation to = 4, but then they diverge, similar to the dust temperature behaviour, suggesting that the FIR luminosity may not be a good tracer of star-formation at high redshift.
Naïve explanations for these findings (after considering heating by the cosmic microwave background; da Cunha et al. 2013) are inadequate. Ambient heating by the elevated cosmic microwave background is subdominant in these objects (da Cunha et al. 2013). Further, at ∼ 7, assuming typical, steady-state dust grains and an interstellar radiation field with total energy density 1500 eV cm −3 , to supply enough energy by star formation would require an SFR of greater than 1000 M yr −1 . Even for a halo as massive (and rare) as 10 13 M , this entails an high, order-of-unity star formation efficiency (Sun & Furlanetto 2016).
However, the dust composition of galaxies at these redshifts is very different from local (e.g. Behrens et al. 2018). Multiple dust components (e.g. cold vs. warm) may be present (Strandet et al. 2017;Liang et al. 2019). As pointed out by De Rossi et al. (2018), dust in primeval galaxies (Pop II) is typically younger than 400 Myr; silicate rich, being produced by low metallicity AGB stars (Ventura et al. 2012), with poor emission efficiency at low and high wavelengths (Koike et al. 2003); and located in environments near to the heating source (Spilker et al. 2016). AGNs are another likely significant source of heating (Lyu et al. 2016), whose host galaxies are very compact (Decarli et al. 2018), with SFRs of 1000 − 5000 M yr −1 (and SFR surface densities reaching 300 − 2000 M yr −1 kpc −2 ; Lyu et al. 2016). A fuller treatment must specifically consider the nature of galaxies in the early Universe.

The LIR and Star-Formation Rate Density
We estimate the far-infrared luminosity density ( LIR ) as a numberweighted sum of the infrared luminosities in different stellar-mass bins derived above, where 2 = 2 2 , and is the cosmic variance term (Moster et al. 2011). Similarly, we estimate the SFR density ( SFR ) from SFRs derived from rest-frame 850 m flux densities.
In Figure 3 we show the LIR as a grey line and shaded region, which when converted to SFR is in good agreement with previous measurements (e.g. Bouwens et al. 2020;Zavala et al. 2021) and models (Béthermin et al. 2017;Pillepich et al. 2018) to = 4, but then remains elevated, while the models decline. Given our finding of hot dust at > 4 this excess is not unexpected, as we know that hot dust SEDs are poor estimators of star-formation rates (Bouwens et al. 2016). On the other hand, the 850 m-derived , shown as a blue line and shaded region, exhibits the expected decline after peaking at = 1 − 2. Note, we exclude the = 8 − 10 bins, whose galaxy abundances suggest a large contribution from interlopers.
It is reasonable to ask how so much of the SFRD can be recovered when the incompleteness levels of the catalogue at > 4 must be high. The reason is because we are stacking on images with large PSFs, which would result in a bias from galaxies that were not in, but were correlated with, those in the catalogue.
The high-SFRD has significant consequences for the evolution of the early Universe (Robertson 2021); leveraging statistical techniques like stacking to access the full unresolved population (see Viero et al. 2015;Cheng et al. 2021) will complement traditional methods to provide a complete picture of star-formation activity.

DISCUSSION
How could so much hot dust in the early Universe have been missed? In hindsight, all the clues were there: there was a dearth of highredshift dusty galaxies (Casey et al. 2018) that hot dust could explain (Sommovigo et al. 2020). There was known tension in low IRX values observed for given at high redshift (Capak et al. 2015) that could be alleviated with warmer dust (Bouwens et al. 2016). There were unusually high dust-mass estimates that challenged dust production time-scales (Leśniewska & Michałowski 2019), which is resolved with LIR coming from hotter dust (Bakx et al. 2020). And then there were simple models showing that ALMA observations tended to miss galaxies with high-temperature dust (Chen et al. 2022).
So why has the idea of hot dust at high redshift not gained traction? First, a handful of hot dusty galaxies at > 6 had been observed (e.g. Behrens et al. 2018;Bakx et al. 2020), but many others were found to be consistent with extrapolations to previous relationships (e.g. Hashimoto et al. 2019;Faisst et al. 2017;Sommovigo et al. 2020;Ferrara et al. 2022), so that the behaviour of the d − relationship -increasing or plateauing -was still unclear (Faisst et al. 2017;Sommovigo et al. 2020). However, those temperature estimates were derived from measurements with ALMA that did not actually bracket the peak of the SED (largely staying below Band 8, so no shorter than 600 m). If in fact these objects were much hotter, no adequate prior to extrapolate model SEDs even existed.
Alongside hot dust is an elevated infrared-luminosity density at > 4, which appears to decouple from the cold-dust derived SFRD. The source of this excess heating is likely some combination of young, Pop II stars and silicate rich dust, and AGN, but in what combination, and is that enough? Such details impact our understanding of reionization time-scales, or of the presence of quiescent galaxies at high redshift (Glazebrook et al. 2017).
Understanding the nature of the dust, and the extreme source of heating -is it stars, AGNs, PBHs? -ranks high on the list of fundamental questions going forward.

CONCLUSION
We estimate far-infrared luminosities and dust temperatures from galaxies at = 0-10, by stacking 111 227 objects from the COSMOS2020/F catalogue on maps in the farinfrared/submillimetre. We find agreement with previous measurements to = 4, but beyond this find elevated temperatures and luminosities, indicating an early Universe filled with hot dust. The temperatures measured suggest a direct detection of Pop II galaxies -compact, close-packed, and silicate rich -whose SEDs carry the signature of localized hot dust.
If confirmed, this will have far-reaching consequences. It will inform future probes into the epoch of reionization, whether directly with the James Webb Space Telescope (JWST), or via statistics of spectral lines like CO and [C ] with line-intensity mapping experiments like COMAP (Cleary et al. 2022), TIME (Sun et al. 2021), or CONCERTO (Catalano et al. 2022). In particular, a warm dust background could present problems for [C ] experiments by decreasing the line-continuum contrast.
Clearly, follow-up observations targeting the peak of hot dust emission are in order. ALMA band 10 will probe nearer to the peak than band 9 observations; targeting the most massive and reliable subset of the > 7 sample should achieve detections with reasonable signal-to-noise, but still will not fully bracket the peak.
Long term, dedicated observatories targeting the peak of warm dust SEDs in the submillimetre will be necessary to complement deep optical/NIR observations. On the horizon are up-to-50 m class ground-based single-dish observatories like CCAT-prime/FYST (CCAT-Prime collaboration et al. 2021), At-LAST (Klaassen et al. 2020), and the Japanese 30 m Antarctic THz telescope intended for New Dome Fuji (Chen et al. 2022), as well as space missions (e.g. Origins Space Telescope; Wiedner et al. 2020).