The early Universe was dust-rich and extremely hot

ABSTRACT

We investigate the dust properties and star-formation signature of galaxies in the early Universe by stacking 111 227 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\pm 12\,\mathrm{ K}$| at |$\mathit{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/.

galaxies: star formation, early Universe, infrared: galaxies, submillimetre: diffuse background

1 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 far-infrared (FIR) wavelengths. To measure the blackbody is to measure star formation; to measure it across cosmic time is to characterize 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 (SFRD)? Do evolutionary trends of galaxies with redshift seen locally [e.g. temperature, specific star-formation rate (SFR)] continue indefinitely, or do they plateau, or even turn over? Just how active are active galactic nuclei (AGNs) in this period? And how much ionizing radiation is escaping the dark curtain of dust exactly?

The catch is that measuring the FIR 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|$\, \rm \mu m$|⁠, or between 400–600|$\, \rm \mu m$| in the observer frame if emitted at the peak of star formation, |$\mathit{z}$| ∼ 2–3. To fully bracket that peak requires combining observations from space telescopes like Spitzer and Herschel, and ground based observatories like the James Clerk Maxwell Telescope (JCMT), the Large Millimeter Telescope (LMT), and the Atacama Large Millimeter/Submillimeter Array (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. Ground-based 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\, \rm m$| mirror in space, the solution is stacking. Galaxies with similar physical properties estimated from optical photometry (e.g. redshift, stellar mass, star-forming 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, python3-compatible version of simstack, to measure the FIR SEDs from redshift 0 to 10 and stellar-masses of log (M/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 |$\mathit{z}$| = 4, but then rise more rapidly to |$\mathit{z}$| ∼ 10. Similarly, we find that the FIR luminosity density (ρ_LIR) tracks the established SFRD or ρ_SFR to |$\mathit{z}$| = 2, but then decouples from it, remaining flat to |$\mathit{z}$| = 8. Given the unusual nature of the hot dust at |$\mathit{z}$| > 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 VI 2020), with Ω_M = 0.315, Ω_Λ = 0.685, |$H_0=67.4\, \rm km\, s^{-1}\, Mpc^{-1}$|⁠, σ₈ = 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.

2 DATA

2.1 COSMOS2020 catalogue

Our analysis is performed on data in the COSMOS¹ field (Scoville et al. 2007), centred at |$\rm 10^h00^m26^s, 2^{\circ }13^{\prime }00^{\prime \prime }$|⁠, and covering 1.606 deg² after masking. We use the recently released COSMOS2020 catalogue (Weaver et al. 2022), whose |$K_s = 25.2\, \rm (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 |$\sim 2\, \rm dex$| in stellar mass.

The COSMOS2020 release comes with two options for photometry; we use farmer/LePhare, 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 NUV − r versus R − J selection (Ilbert et al. 2013).

2.2 Maps

Maps are selected to bracket the rest-frame dust SED at |$\mathit{z}$| = 0–10. Unless otherwise noted, they were obtained from the HELP.²

The Spitzer/MIPS map at 24|$\, \rm \mu m$| is a Spitzer Enhanced Imaging Product (SEIP). It has an rms of |$27.9\, \rm \mu 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⁻¹ to Jy beam⁻¹ by dividing the map by the beam solid angle, |$1.55\times 10^{-9}\, \rm sr$|⁠.

Herschel/PACS maps at 100 and 160|$\, \rm \mu m$| are from the PACS Evolutionary Probe (PEP; Poglitsch et al. 2010). They have rms values of 0.10 and |$2.11\, \rm mJy$|⁠, and FWHM of 7.5 and 11.3 arcsec. PACS maps are also converted from MJy sr⁻¹ to Jy beam⁻¹, with beam solid angles of 2.03 × 10⁻⁹ and |$4.66\times 10^{-9}\, \rm sr$|⁠. Calibration uncertainty is 4 per cent.

Herschel/SPIRE maps at 250, 350, and |$500\, \rm \mu m$| are from the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012), with instrumental depths of 15.9, 13.3, and 19.1 mJy (5 σ); FWHM are 17.6, 23.9, and 35.2 arcsec. Calibration uncertainty is 5 per cent.

The |$850\, \rm \mu m$| map is a product of the S2CLS³ (Geach et al. 2017). Map depth is 1.6 mJy beam⁻¹ (⁠|$1\, \sigma$|⁠), and FWHM 12.1 arcsec. Note, we specifically select the non-match-filtered (NMF) maps, which retain the unresolved fluctuations on which simstack relies.

3 METHOD

3.1 Unbiased stacking with simstack

Our stacking analysis is based on simstack, an algorithm designed to overcome the biases inherent with stacking on images dominated by confusion noise. Briefly, simstack uses the 3 D 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 simstack.

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 simstack, 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 per cent bins are retained as the bootstrapped flux, and the 20 per cent 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.

3.2 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 emcee for python (Foreman-Mackey et al. 2013), a user-friendly, affine-invariant 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 |$\mathit{z}$| = 1.5–2 and |$\mathit{z}$| = 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 polycyclic aromatic hydrocarbon (PAH) emission lines contaminating the 24|$\, \rm \mu m$| flux density by artificially inflating their errors over the redshift range 0.5 < |$\mathit{z}$| < 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 νI_ν from rest-frame 8 to 1000 |$\, \rm \mu m$| (Casey 2012), and naïvely converted to the 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 |$\, \rm \mu m$| luminosity L₈₅₀. As a proxy for the interstellar medium gas mass, it is insensitive to dust temperature T_d with a relatively constant luminosity–mass ratio, and is also a reliable SFR estimator (Scoville et al. 2016).

3.3 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 per cent group of the bootstrap split, finding their relative variance is proportional to their relative size.

Next, we carefully examine the impact that low-|$\mathit{z}$| interlopers have on the stacked SED shape at |$\mathit{z}$| > 4 using simulated catalogs/maps generated with the pysides simulation (Béthermin et al. 2010). We use the redshift probability distributions P(⁠|$\mathit{z}$|⁠) provided by the COSMOS team to i) draw a redshift for each object, ii) look for the closest counterpart in the sides catalog, iii) bin and average their flux densities, iv) and compare the best-fitting 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 P(⁠|$\mathit{z}$|⁠), 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 |$\mathit{z}$| > 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).

4 RESULTS

The full set of SEDs for star-forming galaxies are shown in Fig. 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 Section 3.3, showing a clear departure from existing trends at |$\mathit{z}$| = 0–4.

Figure 1.

Best-fitting SEDs of star-forming galaxies, binned by redshift (increasing left to right) and stellar mass (increasing top to bottom). Data at 24|$\, \rm \mu m$| (blue circles) are from Spitzer/MIPS; at 100 and 160 |$\, \rm \mu m$| (green circles) from Herschel/PACS; at 250, 350, and 500|$\, \rm \mu m$| (red circles) from Herschel/SPIRE; and at 850 |$\, \rm \mu m$| (yellow circles) from SCUBA2. Error bars represent |$1\, \sigma$| Gaussian uncertainties, derived with the bootstrap estimator, on 150 iterations. Non-detections are shown at |$3\, \sigma$| upper limits. The SED model is a hybrid blackbody with fixed emissivity index β = 1.8 on the Rayleigh–Jeans side, and fixed power–law approximation slope α = 2 on the Wien side. MCMC fits (cyan) show the 50th percentile line bracketed by a shaded region marking the 25th and 75th percentiles. For comparison we also show hypothetical main-sequence galaxy SEDs (thin magenta lines), extrapolated from the Schreiber et al. (2015) model at |$\mathit{z}$| < 4.

4.1 Dust temperature

We find dust temperatures increasing with redshift (Fig. 2), in agreement with previous relationships found to redshift 4 (Viero et al. 2013; Schreiber et al. 2018; Bouwens et al. 2020). At |$\mathit{z}$| > 4, we find temperatures diverge from those relationships. Temperatures are incompatible with other estimates that follow (extrapolated) existing relationships (e.g. Hashimoto et al. 2019; Béthermin et al. 2020; Faisst et al. 2020; Ferrara et al. 2022; Sommovigo et al. 2022). They are in better agreement with the lower limit of Bakx et al. (2020), and estimate of Behrens et al. (2018) based on data from Laporte et al. (2017); these findings had been considered outliers. We model the T_d–|$\mathit{z}$| relationship with a polynomial:

$$\begin{eqnarray*} T_\mathrm{d}(\mathit{z}) = (20.9\pm 3.5) + (5.9\pm 1.9)\, \mathit{z} + (0.5\pm 0.2)\, \mathit{z}^2. \end{eqnarray*}$$

(1)

Figure 2.

Rest-frame dust temperatures from best-fittings to SEDs plotted as open circles. A polynomial fit for T(⁠|$\mathit{z}$|⁠) is shown as a dashed red line. Closed circles are the five bins fit with priors, and are not included in the fit. Positions on the redshift axis correspond to median redshifts of the objects in the bins. Also plotted are T_d–|$\mathit{z}$| relationships (Viero et al. 2013; Schreiber et al. 2018; Bouwens et al. 2020), extrapolated to |$\mathit{z}$| = 9. Temperatures of individual objects include Hashimoto et al. (2019), Faisst et al. (2020), Sommovigo et al. (2022), Béthermin et al. (2020), and Laporte et al. (2017)/Behrens et al. (2018), and lower limit from Bakx et al. (2020).

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 |$\mathit{z}$| ∼ 7, assuming typical, steady-state dust grains and an interstellar radiation field with total energy density |$U\gtrsim 1500\, \rm {eV\, cm^{-3}}$|⁠, to supply enough energy by star formation would require an SFR of greater than |$1000\, \rm M_{\odot }\, \mathrm{yr^{-1}}$|⁠. Even for a halo as massive (and rare) as |$\rm 10^{13}\, M_{\odot }$|⁠, 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 versus 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\, \rm 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, Rieke & Alberts 2016), whose host galaxies are very compact (Decarli et al. 2018), with SFRs of |$1000-5000\, \rm M_{\odot }\, yr^{-1}$| (and SFR surface densities reaching |$300-2000\, \rm M_{\odot } \, yr^{-1} \, kpc^{-2}$|⁠; Lyu et al. 2016). A fuller treatment must specifically consider the nature of galaxies in the early Universe.

4.2 The LIR and SFRD

We estimate the FIR luminosity density (ρ_LIR) as a number-weighted sum of the LIR in different stellar-mass bins derived above,

$$\begin{eqnarray*} \rho _\mathrm{LIR} =\frac{1}{V} \displaystyle \sum ^m _{i=0} L_i N_i, ~~ \delta \rho _\mathrm{LIR}^2 = \frac{1}{V^2} \displaystyle \sum ^m _{i=0} (N_i \delta L_i)^2 + (L_i \delta N_i)^2, \end{eqnarray*}$$

(2)

where |$\delta N^2_i = N^2\epsilon ^2_{CV}$|⁠, and ϵ_CV is the cosmic variance term (Moster et al. 2011). Similarly, we estimate the SFR density (ρ_SFR) from SFRs derived from rest-frame 850|$\, \rm \mu m$| flux densities.

In Fig. 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 |$\mathit{z}$| = 4, but then remains elevated, while the models decline. Given our finding of hot dust at |$\mathit{z}$| > 4 this excess is not unexpected, as we know that hot dust SEDs are poor estimators of SFRs (Bouwens et al. 2016). On the other hand, the 850 |$\, \rm \mu m$|-derived ρ_SFR, shown as a blue line and shaded region, exhibits the expected decline after peaking at |$\mathit{z}$| = 1 − 2. Note, we exclude the |$\mathit{z}$| = 8 − 10 bins, whose galaxy abundances suggest a large contribution from interlopers.

Figure 3.

SFRD versus redshift, binned by stellar mass (coloured lines), and summed (black line with 1 σ shaded region). SFRs are converted from FIR luminosities – derived from the integral of intensities from (rest-frame) 8 to 1000|$\, \rm \mu m$|– by scaling with the Kennicutt (1998) relation. We find we are in broad agreement with a selection of the latest |$\rm SFRD$| from models (Béthermin et al. 2017; Pillepich et al. 2018) and measurements (Zavala et al. 2021). Also plotted is the rest-frame 850 |$\, \rm \mu m$|-derived SFRD as a light blue line and shaded region (⁠|$1\, \sigma$|⁠). The LIR-derived SFRD tracks 850 |$\, \rm \mu m$|-derived relation to |$\mathit{z}$| = 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.

It is reasonable to ask how so much of the |$\rm SFRD$| can be recovered when the incompleteness levels of the catalogue at |$\mathit{z}$| > 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-|$\mathit{z}$| 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.

5 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 high-redshift 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 |$\mathit{z}$| > 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. Faisst et al. 2017; Hashimoto et al. 2019; Sommovigo et al. 2020; Ferrara et al. 2022), so that the behaviour of the T_d − |$\mathit{z}$| 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 |$\, \rm \mu 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 |$\mathit{z}$| > 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.

6 CONCLUSION

We estimate FIR luminosities and dust temperatures from galaxies at |$\mathit{z}$| = 0–10, by stacking 111 227 objects from the COSMOS2020/Farmer catalogue on maps in the FIR/submillimetre. We find agreement with previous measurements to |$\mathit{z}$| = 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 ii] 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 ii] 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 |$\mathit{z}$| > 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/near-infrared observations. On the horizon are up-to-50 m class ground-based single-dish observatories like CCAT-prime/FYST (CCAT-Prime collaboration 2021), AtLAST (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).

ACKNOWLEDGEMENTS

The authors are grateful to the referee for their careful reading and constructive suggestions, thank you. DTC holds a CITA/Dunlap Institute postdoctoral fellowship. The Dunlap Institute is funded through an endowment established by the David Dunlap family and the University of Toronto. The University of Toronto operates on the traditional land of the Huron-Wendat, the Seneca, and most recently, the Mississaugas of the Credit River; DTC is grateful to have the opportunity to work on this land.

DATA AVAILABILITY

The simstack package needed to reproduce these results, and Jupyter Notebooks guiding the user through each step, are available at https://github.com/marcoviero/simstack3/tree/main/viero2022/ and at the DOI: https://doi.org/10.5281/zenodo.6792395.

Footnotes

1

https://cosmos.astro.caltech.edu

2

https://herschel.sussex.ac.uk

3

https://zenodo.org/record/57792