key: cord-0292257-rz6ddky4
authors: Vijayan, Aswin P.; Wilkins, Stephen M.; Lovell, Christopher C.; Thomas, Peter A.; Camps, Peter; Baes, Maarten; Trayford, James; Kuusisto, Jussi; Roper, William J.
title: First Light And Reionisation Epoch Simulations (FLARES) III: The properties of massive dusty galaxies at cosmic dawn
date: 2021-08-02
journal: nan
DOI: 10.1093/mnras/stac338
sha: 3bc6054832933f9e20dfe3b8f09409b4f9729ea0
doc_id: 292257
cord_uid: rz6ddky4

Using the First Light And Reionisation Epoch Simulations (textsc{Flares}) we explore the dust driven properties of massive high-redshift galaxies at $zin[5,10]$. By post-processing the galaxy sample using the radiative transfer code textsc{skirt} we obtain the full spectral energy distribution. We explore the resultant luminosity functions, IRX-$beta$ relations as well as the luminosity-weighted dust temperatures in the Epoch of Reionisation (EoR). We find that most of our results are in agreement with the current set of observations, but under-predict the number densities of bright IR galaxies, which are extremely biased towards the most overdense regions. We see that the textsc{Flares} IRX-$beta$ relation (for $5le zle8$) predominantly follows the local starburst relation. The IRX shows an increase with stellar mass, plateauing at the high-mass end ($sim10^{10}$M$_{odot}$) and shows no evolution in the median normalisation with redshift. We also look at the dependence of the peak dust temperature ($T_{mathrm{peak}}$) on various galaxy properties including the stellar mass, IR luminosity and sSFR, finding the correlation to be strongest with sSFR. The luminosity-weighted dust temperatures increase towards higher redshifts, with the slope of the $T_{mathrm{peak}}$ - redshift relation showing a higher slope than the lower redshift relations obtained from previous observational and theoretical works. The results from textsc{Flares}, which is able to provide a better statistical sample of high-redshift galaxies compared to other simulations, provides a distinct vantage point for the high-redshift Universe.

The Hubble Space Telescope (HST) has been instrumental in the last decade observing the rest-frame UV of high-redshift galaxies (e.g. Beckwith et al. 2006; Bouwens et al. 2006; Bunker et al. 2010; Grogin et al. 2011; Salmon et al. 2020) , finding more than 1000 galaxies at z > 5. These efforts from HST have been complemented by widearea ground based near-IR surveys (e.g. UltraVISTA, Bowler et al. 2014; Stefanon et al. 2019 ) providing samples of rare bright galaxies. Spitzer Space Telescope observations (e.g. Ashby et al. 2013; Roberts-Borsani et al. 2016) , probing the rest-frame optical at z > 5 has provided further contraints on these high-redshift systems. However, the UV/optical alone cannot unravel the nature as well as dynamical properties of these high-redshift systems, such as reliable estimates of the total star formation rates, since it is not an unbiased tracer due to the presence of dust.

Dust plays a major role in the observation of galaxies, with nearly half of all photons in the Universe reprocessed by dust grains during their lifetime (e.g. Savage & Mathis 1979; Dwek et al. 1998 ; Bern-E-mail: apavi@space.dtu.dk stein et al. 2002; Driver et al. 2016) . Even though the average dust content of galaxies in the EoR is expected to be lower compared to the low-redshift Universe (z < 4, e.g. Li et al. 2019; Magnelli et al. 2020; Pozzi et al. 2021) , it still has a significant impact on shaping observations by attenuating the emitted source radiation (see review by Salim & Narayanan 2020) , particularly on the most massive galaxies. Recent observations have also found dust-obscured galaxies (as early as z ∼ 7, e.g. Fudamoto et al. 2021) , which were undetected in the rest-frame UV. Thus it is important that we probe the nature of dust and its affects to better understand the various observationally derived quantities of galaxies in the EoR. The stellar emission in a galaxy, which is predominantly in the UV-to-NIR, gets re-processed by the intervening dust into the IR regime. Over the years, the observations in this regime using far infrared (FIR), millimetre (mm) and sub-millimetre (sub-mm) observatories have been instrumental in mapping the dust content of galaxies. This has been done with the help of instruments like ALMA (e.g. Knudsen et al. 2017; Hashimoto et al. 2018; Smit et al. 2018; Bouwens et al. 2020) , Herschel (e.g. Gruppioni et al. 2013; Wang et al. 2019) , etc. In many cases there have been detections from deep ALMA and PdBI observations of galaxies at extremely high redshifts (z > 6) with large reservoirs of dust (> 10 8 M ; Mortlock et al. 2011; Venemans et al. 2012; da Cunha et al. 2015) .

The early identification of these dusty star-forming galaxies were from single-dish sub-mm surveys finding massive populations at z > 1 (see Casey et al. 2014) . Even though they are rare, these galaxies contribute significantly to the cosmic star formation density during cosmic noon (z ∼ 2 − 3). The picture at higher redshift (z > 4) is still unclear. ALMA and Herschel have been instrumental in filling this space at high-redshift. Recent survey programmes like ALPINE (Le Fèvre et al. 2020; Béthermin et al. 2020; Faisst et al. 2020a ) and ASPECs (Walter et al. 2016; Decarli et al. 2019; González-López et al. 2019 ) are helping us to understand the dusty nature of highredshift galaxies (also see Hodge & da Cunha 2020, for more high-z surveys) by building a large statistical sample. High-redshift studies like Gruppioni et al. (2013) ; Wang et al. (2019) ; Gruppioni et al. (2020) have constructed IR luminosity functions. The jury is still out on the normalisation of the IR LF due to difficulties in de-blending of IR data and smaller volumes probed in some surveys. Other observational studies like Schreiber et al. (2018) and Bouwens et al. (2020) have explored the evolution of the luminosity-weighted dust temperatures, and have found an increase in the value with increasing redshift (up to z ≤ 5). Another important observational space that has been studied widely is the relationship between the Infrared Excess (IRX) -UV continuum slope (β ). Empirical relationships (e.g. Pettini et al. 1998; Meurer et al. 1999; Takeuchi et al. 2012; Reddy et al. 2015) built in this space using observations of low-redshift (mainly z 2) galaxies have been used to correct for dust attenuation in galaxies. There has been a variety of observational studies exploring this space at these high redshifts (e.g. Koprowski et al. 2018; Fudamoto et al. 2020; Bouwens et al. 2020; Schouws et al. 2021) . They have found varying results that favours empirical relation using Calzetti as well as SMC extinction curve. The obtained relation is also strongly influenced by the adopted SED dust temperature and the functional form used to obtain the IR luminosity. With upcoming surveys on facilities like the James Webb Space Telescope (JWST), Euclid, Roman Space Telescope, and the Atacama Large Aperture Submillimeter Telescope (AtLAST, Klaassen et al. 2019 ) are expected to substantially contribute to these efforts to build a comprehensive picture of galaxy formation and evolution in the highredshift Universe.

In addition to these observational efforts, it is crucial to study the nature of these dusty high-redshift systems using theoretical models of galaxy formation and evolution. Several studies have used techniques built with semi-analytical or analytical (e.g. Lacey et al. 2016; Popping et al. 2017; Lagache et al. 2018; Lagos et al. 2019; Sommovigo et al. 2020 ) and hydrodynamical (e.g. Olsen et al. 2017; Narayanan et al. 2018; Ma et al. 2019; McAlpine et al. 2019; Liang et al. 2019; Baes et al. 2020; Trčka et al. 2020; Lovell et al. 2021b; Shen et al. 2022 ) models to explore and understand the trends and variations in observed properties of galaxies like the dust temperatures, fine-structure transitions, infrared luminosity functions, IRXβ relations, submillimeter number counts, etc. Many of them have been successful in reproducing various observational results, and has also been instrumental in understanding the underlying scaling relations.

In comparison to semi-analytical models, hydrodynamical simulations model in greater detail the evolution of dark matter, gas, stars, and black holes, allowing for a more detailed exploration of galaxy structure and observed properties. A drawback of some of the current state-of-the-art cosmological simulation periodic boxes [e.g. MASSIVEBLACK (Matteo et al. 2012) , EAGLE Crain et al. 2015) , BLUETIDES (Feng et al. 2016 ), ILLUSTRIS-TNG (Naiman et al. 2018; Nelson et al. 2018; Marinacci et al. 2018; Springel et al. 2018; Pillepich et al. 2018) , SIMBA , COSMIC DAWN II (Ocvirk et al. 2020) , etc] is that they have fewer massive galaxies in the EoR, which are thought to be biased towards the most overdense regions (see Chiang et al. 2013; Lovell et al. 2018; Ito et al. 2020 ). This is mainly due to the unfeasible amount of computational time required to run much larger periodic volumes with the needed resolution to resolve the relevant scales at these high redshifts. Hence, they lack the statistical power to investigate the bright galaxies that will be discovered and investigated with the current or future generation of telescopes.

To overcome this dearth of a representative statistical sample of massive galaxies when studying the EoR, we have run a suite of zoom simulations, the First Light And Reionisation Epoch Simulations, FLARES; introduced in Lovell et al. (2021a); Vijayan et al. (2021) to re-simulate a wide range of overdensities. The simulations were run using the well studied EAGLE Crain et al. 2015) model, which has been shown to be in good agreement with a large number observables at low and high redshift regime not used in the calibration (e.g. Furlong et al. 2015; Lagos et al. 2015; Trayford et al. 2017) . In this paper, using the radiative transfer code SKIRT (Camps & Baes 2015) we post-process the simulation to study the spectral energy distributions (SEDs) of massive galaxies in the EoR. We aim to understand how well the EAGLE physics model is able to reproduce the high-redshift Universe, mostly in comparison to observations in the infrared part of the spectrum like the luminosity functions, IRX-β and luminosity-weighted dust temperatures. This work complements other theoretical studies in the high-redshift Universe and provide insights into how the intrinsic galaxy properties are connected to their observed and derived properties. This paper is structured as follows, in section §2 we introduce the suite of simulations, our galaxy sample and the method for SED generation. In section §3 we show our results, including the UV and IR luminosity in §3.1, the IRX-β relation in §3.2, and the variation and evolution of dust temperatures in §3.3. We finally summarise our findings and present our conclusions in section §4. Throughout we assume a Planck year 1 cosmology (Ω m = 0.307, Ω Λ = 0.693, h = 0.6777; Planck Collaboration et al. 2014).

We use the First Light And Reionisation Epoch Simulations (FLARES, Lovell et al. 2021a; Vijayan et al. 2021 ), a suite of 40 zoom simulations chosen from a 3.2 cGpc a side dark matter only box resimulated with full hydrodynamics using the EAGLE Crain et al. 2015) simulation physics, down to z = 4.7. The simulations were run with a heavily modified version of P-GADGET-3, which was last described in Springel et al. (2005a) , an N-Body Tree-PM smoothed particle hydrodynamics (SPH) code. The model uses the hydrodynamic solver collectively known as ANARCHY (described in Schaye et al. 2015; Schaller et al. 2015) , that adopts the pressure-entropy formulation described by Hopkins (2013) , an artificial viscosity switch (Cullen & Dehnen 2010) , and an artificial conduction switch (e.g. Price 2008) . The model includes radiative cooling and photo-heating (Wiersma et al. 2009a) , star formation (Schaye & Dalla Vecchia 2008) , stellar evolution and mass loss (Wiersma et al. 2009b) , black hole growth (Springel et al. 2005b ) and feedback from star formation (Dalla Vecchia & Schaye 2012) and AGN (Springel et al. 2005b; Booth & Schaye 2009; Rosas-Guevara et al. 2015) . We use the AGNdT9 configuration of the EAGLE simulation physics, which produces similar mass functions to the Reference model but better reproduces the hot gas properties of groups and clusters (Barnes et al. 2017 ). This configuration gives less frequent, more energetic AGN outbursts.

FLARES has an identical resolution to the EAGLE reference run, with a dark matter and an initial gas particle mass of m dm = 9.7 × 10 6 M and m g = 1.8 × 10 6 M respectively, and has a gravitational softening length of 2.66 ckpc at z ≥ 2.8. The regions are selected at z = 4.7 and have a radius of 14 cMpc/h, spanning a wide range of overdensities, from δ = −0.479 → 0.970 (see Table A1 in Lovell et al. 2021a) . We have deliberately selected a large number of extreme overdense regions (16) to obtain a large sample of massive galaxies. To obtain a representative sample of the Universe, the galaxies in various overdensities can be combined using an appropriate weighting scheme. For a more detailed description of the simulation and weighting method we refer the readers to Lovell et al. (2021a) . Since the galaxies analysed in this study are sampled from the different regions, the mean or median results quoted are weighted values based on that scheme.

Galaxies in FLARES, similar to the standard EAGLE, are identified with the SUBFIND algorithm (Springel et al. 2001; Dolag et al. 2009 ), which runs on bound groups found from via the Friends-Of-Friends algorithm (FOF, Davis et al. 1985) . The galaxy stellar masses are defined using star particles within a 30 pkpc aperture centred on the most bound particle of the self-bound substructures. For the purpose of this study we concentrate only on the most well resolved galaxy systems that have more than 1000 star particles. This coincides with galaxies more massive than ∼10 9 M in stellar mass (see Fig. 1 ). This selection also overlaps very well with the observationally inferred mass ranges (for e.g. the ALPINE survey; Le Fèvre et al. 2020; Béthermin et al. 2020; Faisst et al. 2020a ) of the galaxies detected in the infrared with ALMA and other instruments.

In Fig 1 we plot the stellar mass of the selected galaxies against their star formation rate (SFR, quoted values are averaged for stars formed in the last 10 Myr) for z ∈ [5, 10]. We also show histograms of the galaxy stellar masses and SFR distributions. Our selection samples ∼ 7000 galaxies in this redshift and mass regime. At z = 10, our sample of galaxies is only 44, and thus any inferences drawn can be subject to large scatter. The SFR seen in our selection has a maximum value just below 10 3 M /yr.

There are various methods to obtain the SED of a galaxy in simulations. A comprehensive method to get them, so that the properties of the dusty medium is captured, is to perform radiative transfer. There are numerous codes [e.g. SUNRISE (Jonsson 2006) , RADMC-3D (Dullemond et al. 2012 ), SKIRT (Camps & Baes 2015) , POW-DERDAY , etc] available, most relying on sophisticated Monte Carlo methods. For this study we use the publicly available code SKIRT, version 9 (Camps & Baes 2020) .

FLARES does not inherently model dust formation and destruction, and thus cannot reliably estimate the amount, nature and distribution of dust in the different galaxies. For the purpose of obtaining the amount and distribution of dust we assume a constant dust-tometal ratio (DTM=M dust /(M metal + M dust )) per galaxy, in SPH gas particles below temperatures of 10 6 K or in star-forming gas particles. This temperature is higher than what was adopted in previous Figure 1 . Shows the relationship between the galaxy stellar mass and the star formation rate (SFR) averaged over the star particles that were formed in the last 10 Myr for z ∈ [5, 10]. Also shown is the histogram of the distribution of stellar mass and SFR in different bins for these redshifts. The total number of galaxies at these redshifts are indicated within brackets alongside the legend. EAGLE-SKIRT work (e.g. Camps et al. 2016; Trayford et al. 2017) , and ensures that dust is only destroyed in the very hot gas phase in the galaxies due to thermal sputtering (see Draine & Salpeter 1979; Tielens et al. 1994 , for more details). Changing the threshold to lower temperatures has negligible impact on the results presented in this work. The DTM ratio is calculated using the DTM fitting function in Vijayan et al. (2019, Equation 15 in that work), obtained from the dust model implemented in the L-Galaxies semi-analytical model. In that work, the DTM ratio is parameterised as a function of the mass-weighted stellar age and the gas-phase metallicity. Fig. 2 shows the evolution and spread in the DTM ratio used in this work as a function of the galaxy stellar mass for z ∈ [5, 10]. It can be seen that there is an increase from a value of ∼ 0.01 at z = 10 to ∼ 0.2 by z = 5. The spread in the value also increases with decreasing redshift. More details on the evolution of the DTM ratio with redshift, and its dependence on various other galaxy properties, can be found in Vijayan et al. (2019) . The use of a varying DTM ratio dependent on galaxy properties, as opposed to a constant value of 0.3, is another difference from previous EAGLE-SKIRT work. The evolution of the median DTM ratio with redshift seen here is similar to the one observed in Vogelsberger et al. (2020) for ILLUSTRIS-TNG galaxies.

In this work we use the Small Magellanic Cloud (SMC) dust grain type and size distribution (Weingartner & Draine 2001 ) built in to SKIRT. Due to its low-metallicity, the SMC is considered to be a good analogue to high-redshift galaxies. We use 8 grain size bins for silicate and graphite dust types to compute the thermal emission. The dust grid for this setup is constructed using the built-in octree grid in SKIRT, using the previously defined dust particle distribution obtained from SPH gas particles. The octree is refined between a minimum refinement level of 6 and maximum of 16, with the cell splitting criterion set to a dust fraction value of 2×10 −6 times the total dust mass in the domain, as well as a maximum V-band optical depth of 10. We use 10 6 photon packets per each radiation field wavelength grid, giving good convergence in observed properties. Our radiation field wavelength grid, as well as the dust emission grid, is spanned by a logarithmic grid between 0.08 − 1500 µm, with 200 points. We include dust self-absorption and re-emission in the set-up, with this procedure iterated such that the change in the absorbed dust luminosity is less than 3 per cent. We place our detector to record the SED at a distance of 1 pMpc enclosing a 60 pkpc a side square region. In this work we record multiple orientation sightlines, but the fiducial orientation is along the z-axis.

Similar to previous EAGLE-SKIRT work, we apply differing amount of dust attenuation to old and young stellar populations. Young stars, with stellar ages less than 10 7 yr, are still embedded in their birth clouds, and as such experience higher dust attenuation (e.g. Charlot & Fall 2000) . We perform the same resampling technique that was employed in Camps et al. (2016) ; Trayford et al. (2017) to designate young and old stellar population from star particles and star-forming gas particles in the simulation. A difference from those works is that we do not subtract the contribution of dust from young stars which were part of the star-forming gas particles when we perform the resampling. This is because these particles already have gas/dust intrinsic to them (see section 2.4.4 in Camps et al. 2016 , about introducing 'ghost' gas particles) unlike the resampled star particles which were converted to young stars. The emission from old stellar populations is modelled using the BPASS (Stanway & Eldridge 2018) SPS library and the young stars with MAPPINGS III (Groves et al. 2008) templates. The former uses the Chabrier (Chabrier 2003 ) IMF while the latter uses Kroupa IMF (Kroupa 2002) . We do not expect this difference to have a big effect since they are very similar. The BPASS model is characterised by the age and metallicity of the stellar particle while the MAPPINGS III template uses the SFR, metallicity, the pressure of the ambient ISM, the compactness of the HII region (log 10 (C)), and the covering fraction of the associated photo-dissociation region ( f PDR ) of the star-forming particles. We use the same prescription for deriving the SFR, pressure and log 10 (C) of the star particle as in Camps et al. (2016) . However, we set the PDR covering fraction, f PDR to 0.2, higher than 0.1 which was used in Camps et al. (2016) . Our adopted value is same as the fiducial value used in Groves et al. (2008) ; Jonsson et al. (2010) . A higher f PDR results in more of the stellar emission to be absorbed by the dust present within the birth-clouds, implying that more of the light is re-processed to the IR. Thus a higher value of f PDR implies a higher value of the IR luminosity. However, it should be noted that some features of the galaxy SED (e.g. change in the position where the IR flux density peaks) also depends on the value of log 10 (C) (also see §5.4 in Liang et al. 2019) . For more details about the SKIRT set-up we have used, we refer the interested reader to Camps et al. (2016) and Trayford et al. (2017) .

We use the local thermal equilibrium set-up in SKIRT which means that the dust grains are in local equilibrium with the radiation field. This condition (as opposed to being in non-thermal equilibrium) will pre-dominantly affect the fluxes in the rest-frame mid-IR, but have very negligible effect on our predictions in this work (for more details see Appendix A, where we have run SKIRT with the non-thermal equilibrium setup). We also include dust heating from CMB radiation, which at high-redshifts (since, T CMB (z) =T CMB (z = 0) × (1 + z); also see da Cunha et al. 2015) can be non-negligible. We do not include the effect of AGN on the SEDs (SKIRT has the capability to model AGN emission, see Stalevski 2012; Stalevski et al. 2016 ); we will show in Appendix B how the predictions are affected when adding the AGN bolometric luminosity to the infrared emission (the effect is negligible and only seen at the bright IR luminosity end). In a future work (Kuusisto et al. in prep) we will explore in more detail the effect of AGN on the UV emission from FLARES galaxies.

We had previously modelled the UV to near-IR SED of the FLARES galaxies in Vijayan et al. (2021) using a line-of-sight (LOS) dust extinction model. That work calibrated the dust attenuation based on matching to the UV luminosity function and the UV luminosity-UV-continuum slope relation at z = 5, as well as the [OIII]λ 4959, 5007+Hβ equivalent width relation at z = 8. Here we do not perform any calibration, and only adopt the dust-to-metal ratio from the L-Galaxies SAM which was successful in reproducing many of the seen observational trends. This will enable us to better understand many of the successes and shortcomings of the EAGLE model when applied at high-redshift. We compare the UV luminosity of the galaxies from this model to the LOS model in Appendix E.

In this section we will look at what we can learn about the dust properties of massive high-redshift galaxies from the FLARE simulations, focussing on z ∈ [5, 10]. In §3.1 we will look at the infrared (IR) luminosity function, while exploring the IRX-β space in §3.2. In §3.3 we will look at the dust temperatures of these galaxies, exploring both the SED-inferred as well as the peak dust temperatures. All these observables are also compared to the current observations. It should be noted that we do not model any observational effects (such as modelling the PSF or associated noise) that are inherent to the observed datasets that we compare to; this could impact derived properties and the associated systematic errors.

In Fig. 3 we show the observed UV luminosity (measured at 1500Å, plotted in magnitudes) function for the FLARES galaxies for z ∈ As can be seen, our model reproduces the UV LF reasonably well within the scatter seen in the observational data for M 1500 −21.

The turnover at the faint-end is mainly due to our selection of wellresolved massive galaxies, whose contribution are at the bright end. It should be noted that there is a hint of galaxy number densities being slightly lower at z = 8 compared to the observations. There is also a similar trend at z = 10, but is harder to draw conclusions from, since the Oesch et al. (2018) data contains only 5 galaxies. Some of this tension can be attributed to the slightly lower normalisation (∼ 0.3 dex) of the SFR function of the EAGLE reference volume or FLARES at intermediate SFR (1 < SFR < 10 M yr −1 ) as noted for high-redshift galaxies in Katsianis et al. (2017) ; Lovell et al. (2021a , also see Furlong et al. 2015 ) when compared to observed values. Vijayan et al. (2021) also showed that the unobscured SFR density of FLARES galaxies at z ∈ [5, 7] showed slightly lower normalisation (∼ 0.2 dex) in comparison with the unobscured value from Bouwens et al. (2020) , even though the dust model was explicitly calibrated to match the UV LF, indicating that either the star formation rates are generally lower in the simulation, or the chemical enrichment rate (and thus the derived dust content) is higher, giving rise to higher attenuation than expected in these model galaxies. In the future with JWST we will be able to put tighter constraints on galaxy metallicities in the high-redshift regime. There is really good agreement at the high UV luminosity end at all the redshifts. Since the UV LF is predicted considerably well against observations (with the caveats noted) we will now try to draw meaningful conclusions from comparing against other observational spaces.

In Fig. 4 we show the IR luminosity functions for z ∈ [5, 10]. The IR luminosity of the FLARES galaxies are obtained by integrating the observed SED between rest-frame wavelength of 8 − 1000 µm. We also plot alongside observational data from Gruppioni et al. (2013, using Herschel data); Wang et al. (2019, using the Herschel catalogue generated by the Bayesian source extraction tool XID+ in the COSMOS field); Gruppioni et al. (2020, using the the ALPINE-ALMA data) as well as theoretical results from FIRE-2 (Ma et al. 2019, their IR LF fit obtained from running SKIRT) for similar redshifts. The FIRE-2 IR LF fit is only plotted till L IR = 10 12 L due to the bulk of their sample being mostly lower mass galaxies compared to FLARES. It can be seen that FLARES is in agreement with the observational data for luminosities 10 12 L for z = 5. There is a sharp decline in extremely bright IR galaxies in our simulation at z = 5. The very bright end of the function is under-estimated by ∼ 1 dex (similar to the FIRE-2 IR LF fit) compared to Gruppioni et al. (2013 Gruppioni et al. ( , 2020 , which are collated measurements within broad redshift ranges. Due to this broader redshift range, the normalisation can be higher, since lower redshifts are expected to have higher number densities. However, a difference of ∼ 1 dex is in tension with our predictions. Zavala et al. (2021) have also described the IR LF measurements in Gruppioni et al. (2020) to be representative of an overdense patch in the high-redshift Universe. This inference . The rest-frame 250 µm luminosity function of the FLARES galaxies for z ∈ [5, 7]. The errorbars show the Poisson 1σ uncertainties for the different bins. Bins with fewer than 5 galaxies are represented by dashed lines. The data is incomplete at the faint-end due to our galaxy selection. We plot alongside observational data from Koprowski et al. (2017) ; Gruppioni et al. (2020) at similar redshift range.

comes from the observational targets being highly clustered massive galaxies (log 10 (M/M ) 10.5). This is in good agreement with the plotted IR LF of highly overdense regions in FLARES shown in Fig. 6 . In case of the Wang et al. (2019) data at z = 5 there is a similar case of underprediction of bright IR luminous galaxies. Thus we are inconsistent in a regime where two independent measurements (Gruppioni et al. 2013; Wang et al. 2019 ) agree, though it should also be noted that they are both obtained from the Herschel catalogue, and are subject to uncertainties associated with the deblending techniques employed. Thus they can be ideally treated as upper limits on the IR luminosity function. A reason for this sudden decrease is that our extreme IR-bright galaxies are biased towards the most overdense regions, having much lower contribution to the IR LF (see §3.1.1). Following from our argument in the UV LF section, the lower normalisation of the star formation rate function in our model at these redshifts also contributes to this lower number density (also discussed in McAlpine et al. 2019; Baes et al. 2020 ). This has also been investigated at lower redshifts and has been similarly attributed to the lower star formation rate as well as the lack of 'bursty' star formation in the EAGLE model (see McAlpine et al. 2019) . However, our result is not an isolated case and has been a feature of many other cosmological and zoom simulations like ILLUSTRIS-TNG (Shen et al. 2022 ) and FIRE-2 (Ma et al. 2019 , also plotted in Figure 4 ) at these redshifts. The SIMBA suite of simulations shows a higher normalisation of the SFR function than EAGLE at high-redshift. In Lovell et al. (2021b) , where they post-processed SIMBA galaxies using POW-DERDAY, ) they find reasonable agreement with observationally inferred 850 µm number counts, which they partly attribute to the higher SFRs. So the lower star formation rate at high-redshift is a likely cause of the deficit in IR luminosities in our model.

There are caveats that come along with physics recipes to produce higher star formation. For example, there is a lack of quiescent galaxies at high-redshifts in SIMBA compared to observations and EAGLE, as explored in Merlin et al. (2019) . A fine interplay of feedback and star formation is fundamental to match the various observational results and thus provide test beds for improving the model recipes. There has also been suggestions of changes to the initial mass function to a top heavy one in the most luminous galaxies to produce the seen higher number density of IR luminous galaxies (see for e.g. Motte et al. 2018; Schneider et al. 2018; Zhang et al. 2018) . Any of these two scenarios would imply higher dust content from increased star-formation, thus reconciling the increase in the intrinsic emission with higher attenuation.

The FLARES IR LF at z = 6 is in very good agreement with the Wang et al. (2019, also similar to what was seen for the IR LF of EAGLE galaxies in that study) and FIRE-2 results, while still being more than ∼ 0.5 dex lower compared to the Gruppioni et al. (2020) data for 4.5 ≤ z ≤ 6. At z = 7 the FLARES IR LF shows a higher normalisation (∼ 0.3 dex) compared to the fit from FIRE-2 in the range 10 11 − 10 12 L . There are no observational constraints on the IR LF for < 10 11 L at z > 4 yet, hence it is hard to draw conclusions from the FIRE-2 data.

In Fig. 5 we plot the rest-frame 250 µm luminosity function of the FLARES galaxies in z ∈ [5, 7] . We also compare to observational data from Koprowski et al. (2017, using sub-mm/mm imaging from SCUBA-2 and ALMA); Gruppioni et al. (2020) . In general, we see an underprediction of the rest-frame 250 µm luminosity function. We are in agreement with the 3.5 < z < 4.5 Gruppioni et al. (2020) data within the uncertainties, while at z = 6, our results are ∼ 1 dex lower at the extreme bright end. It can also be seen that the Gruppioni et al.

(2020) data in the two redshift range show no clear decline in the number density galaxies, while in our case there is a reduction in number density by ∼ 0.5 at the bright end of the function. But this is inconsistent with Koprowski et al. (2017) where they found a lower number density in the range 3.5 < z < 4.5 compared to the Gruppioni et al. (2020) values. This discrepancy in the two data sets could be due to incompleteness in the sample selection associated with the Koprowski et al. (2017) data. Nevertheless, our data does not extend to the extreme luminosities that the Koprowski et al. (2017) sample covers.

In Appendix B we add the AGN bolometric luminosity to the IR luminosity for comparison to observations to gauge the effect AGN has on the IR LF. We see very small changes, not enough to reconcile an order of magnitude difference at the bright end of the IR LF. Also, in Appendix C, we look at how our results are consistent with the current setup when excluding emission from birth clouds of young stars for z ≥ 8.

In Fig. 6 , we show the IR LF in different matter overdensity bins for z ∈ [5, 7]. The composite function sits within the boundary of the positive and negative overdensity bins as expected. The plot shows that there is an increase in the number densities of IR luminous galaxies with increasing overdensity. The smallest overdensity bin is not visible in the plot since it is below the plotted IR luminosity range. At z = 5, there is an increase of ∼ 1.5 dex in number density of galaxies with IR luminosity of ∼ 10 11 L . This is very similar across the rest of the plotted redshift range. At z = 5, only the most overdense regions contribute to the very bright end of the IR LF, which results in the rapid fall in the number densities seen in the composite function.

In this section we will look at the infrared excess (IRX) of the galaxies in the FLARE simulation. The infrared excess is defined as the ratio of the total infrared luminosity (L IR ) over the UV luminosity (L UV ), and can be derived as follows

where L 1500 = L 1500 × 1500 Å, with L 1500 being the far-UV luminosity calculated at 1500 Å. The UV-continuum slope, β is defined such that f λ ∝ λ β or alternatively f ν ∝ λ β +2 , for λ in the rest-frame UV range. We measure β using the following prescription, β = log 10 (L 1500 /L 1500 ) log 10 (1500/2500) − 2,

where L 1500 and L 2500 are the far-UV and near-UV luminosity, respectively.

The IRX-β relation has been explored in numerous theoretical (e.g. Popping et al. 2017; Safarzadeh et al. 2017; Ma et al. 2019; Trčka et al. 2020; Schulz et al. 2020; Liang et al. 2021 ) and observational studies (e.g. Reddy et al. 2015; Wang et al. 2018; Fudamoto et al. 2020; Bouwens et al. 2020) in low-and high-redshift galaxies. Some studies have provided empirical relations for this plane assuming a dust screen model. These relations are expected to arise from the simple assumption that with increasing dust attenuation the UV-continuum slope becomes redder with the L IR to L UV ratio increasing. The assumption of different dust attenuation curves will Bouwens et al. (2020) . Also shown is the empirical relation for dust screen models for SMC (Pettini et al. 1998 ) and Calzetti (Meurer et al. 1999) , and from Reddy et al. (2015) .

determine the trajectory of this relation. However, in galaxies one would expect different attenuation for young and old stars, as well as different dust distribution across the galaxy that can provide large variation to the relationship (e.g. Narayanan et al. 2018; Schulz et al. 2020) . The empirical relation between IRX and β is widely used to correct for the amount of dust-obscured star formation within galaxies at high-redshift. This is based on the assumption that high-redshift galaxies follow the same relation as their local analogues. However, some high-redshift galaxies seem to have smaller IRX than their local analogues (e.g. Capak et al. 2015; Fudamoto et al. 2020 ). This has been largely attributed to the low dust temperatures adopted in modelling the observational data (e.g. Sommovigo et al. 2020; Bouwens et al. 2020 , by changing the usual assumption of dust tem- peratures of 35 − 45 K to 45 − 50 K). Another cause of concern in using this relation at high redshift is the observed spatial offset between the UV and IR emission (e.g. Bowler et al. 2018) , possibly due to the birth cloud dispersal time in galaxies (Sommovigo et al. 2020) .

We look at the variation of the IRX-β relation across z ∈ [5, 10] in Fig. 7 . The plane is represented by hexbins which are coloured by their median sSFR values. Due to our simulations containing a large selection of extreme overdensities, there is an overabundance of IR luminous dusty galaxies in our data. We also plot observational data at similar redshifts from Capak et al. . Also shown is the empirical IRX-β relations from Pettini et al. (1998, SMC) , Meurer et al. (1999, Calzetti) and Reddy et al. (2015) .

From the figure we can see that IRX-β relation of the FLARES galaxies agree well with plotted observational values. There are a few exceptions in case of data with high-β and low-IRX values (low dust content and older stellar populations) found for a few galaxies in the Hashimoto et al. (2019) and Harikane et al. (2020) sample. This could be a drawback of our implemented model that does not fully capture the diverse star-dust geometries of galaxies in these observations. These galaxies having high-β and low-IRX imply that they are moving away from the main-sequence relation towards quiescence. The FLARES sample could be inefficient in producing such galaxies. The quiescent galaxy population in FLARES will be probed in a future work where their number densities will also be explored. A similar dearth of high-β , low-IRX galaxies is seen in Ma et al. (2019) . However, their probed galaxy stellar mass range is lower than ours and thus it could also be due to the fact that there are not any low-sSFR high-mass galaxies, that are moving into the quiescent regime.

At high-redshifts (z > 4), no clear picture has emerged on what kind of attenuation relation galaxies follow, whether or not they are consistent with the local relation, following a starburst or Calzettilike attenuation law, or a steeper one like the SMC. As can be seen in Fig. 7 , the massive galaxies in FLARES predominantly follow the Calzetti or Reddy et al. (2015) like relation, with a small proportion of galaxies following the SMC relation at z ∈ [5, 8], even though we adopted the SMC grain distribution. We also note that the galaxies which drop below the canonical relations are the ones that exhibit very low sSFR values, or with older stellar populations, in agreement with theoretical studies like Narayanan et al. (2018) . At higher redshifts (z ≥ 8) there is a hint of some transition away from the Calzetti relation towards the SMC one in FLARES, which we fail to properly capture due to the limited mass resolution of our simulation. This leads to the lack of well resolved low-mass galaxies in our sample to populate this space.

The reason for the majority of FLARES galaxies following the local starburst relation is mostly because these high-redshift galaxies are intrinsically bluer with higher sSFR compared to their lowredshift counterparts. The inhomogeneous nature of dust distribution at high redshift will also have an effect, with the β values being dominated by unobscured young stars while the IRX is dominated by dust emission near the highly obscured dust patches. With our selection of the most massive galaxies, this is a very likely outcome due to their higher dust content. Narayanan et al. (2018) , using cosmological zoom simulations run using GIZMO (Hopkins 2015) , also attributed the major drivers of the observed deviations from the canonical relation to older stellar populations, complex star-dust geometries and variations in dust extinction curves. A similar result was also seen at high-redshift in Schulz et al. (2020) , studying the IRX-β relation in ILLUSTRIS-TNG galaxies at z = 0 − 4. They concluded that the seen deviations could be best described in terms of sSFRs driving the shift in β with the star formation efficiency possibly being a good indicator of the variations in star-dust geometry. Liang et al. (2021) explored in detail the secondary dependencies of the IRX-β relation, concluding that the main driver of the scatter is the variations in the intrinsic UV spectral slope and thus the age of the underlying stellar population. Thus, due to large degeneracies among sSFRs or ages, star-dust geometry as well as dust compositions, it would be hard to pin-point a global track in the IRX-β relation for galaxies at these redshifts for different stellar masses.

In Fig. 8 we look at the relationship between the galaxy stellar mass and IRX. We also plot observational results from the publicly available dataset of the ALPINE collaboration (Le Fèvre et al. 2020; Béthermin et al. 2020; Faisst et al. 2020a , values obtained from SED fitting) as well as the stacked and median results from Fudamoto et al. (2020) and Bouwens et al. (2020) respectively. Also shown is the weighted median and the 16-84th percentile variation for the sample of galaxies. This plane provides further insights into what we saw in Fig. 7 . We can see that the galaxies following the SMC relation have stellar masses of ∼ 10 9 M . These galaxies are the ones in the process of transitioning to the higher attenuation relation from rapid dust enrichment (see figs 9 and 10 in Vijayan et al. 2021 where we also see a rapid rise in the UV attenuation ∼ 10 9 M in stellar mass). We also see that there is a general trend towards high median IRX values with higher stellar masses, plateauing or slightly dropping at the highest masses. If we extrapolate the median to lower masses, it would lead to lower IRX values. This would indicate that the region near the SMC relation would be occupied by the lower mass galaxies. This is in agreement with what was seen in Ma et al. (2019) using the FIRE-2 simulation, where the majority of the galaxies below a stellar mass of 10 9 M follow the SMC relation.

The trend at the massive end (∼ 10 10 M ), which is clear at 5 ≤ z ≤ 8, points towards an increase in the UV luminosity not being reflected to the same extent in the IR luminosity. This points towards decreasing dust attenuation in the most massive galaxies. This was also seen in the LOS dust attenuation model applied on the FLARES galaxies in Vijayan et al. (2021, see fig. 10 ). The observed UV LF being better fit by double power-law at these high-redshift (e.g. Bowler et al. 2014 Bowler et al. , 2020 Shibuya et al. 2022 ) also points towards decreasing dust attenuation at the bright and massive end, that can contribute to this decline. However, studies such as Ferrara et al. (2017) posit that some of the IRX deficit galaxies could have dust embedded in large gas reservoirs, thus not contributing to an increase in the IR luminosity.

On comparing to the observational data, the median values from Bouwens et al. (2020, note that the highest mass bin has just one galaxy) are a good match. In case of the ALPINE data as well as the stacked results from Fudamoto et al. (2020, which also uses ALPINE data), our median relation is higher than their dataset. However, this can be explained by the UV selection of the galaxies observed in the survey, which can miss the dustier systems.

There is no noticeable evolution of IRX-stellar mass relation with redshift. The normalisation of the median value does not show any redshift evolution.

In this section we will look at the dust temperatures of the galaxies in FLARES. There are different definitions of the dust temperature, both observational and theoretical. Here we will look at the SED-derived and peak wavelength temperature, which are both measures of the light-weighted dust temperatures, but can be considerably different depending on the functional forms being used (see e.g. Casey 2012). These measures are very different to the mass-weighted temperature, a measure of mostly the cold dust content of galaxies, which is expected to be largely independent of redshift and galaxy properties as well as significantly lower than the light-weighted temperatures (Liang et al. 2019; Sommovigo et al. 2020 ).

The peak dust temperature (T peak ) can be obtained from the Wein displacement law from the rest-frame wavelength at which the infrared flux density peaks (λ peak ). This is defined as

which follows the relation for a true blackbody (that has a dust emissivity index β of 2). This measure has been used in many observational studies to understand the evolution of luminosity-weighted dust temperature across redshifts and galaxy properties (e.g. Casey et al. 2018a; Schreiber et al. 2018; Burnham et al. 2021) . It suffers less from model dependent biases compared to other dust temperature values. It should also be noted that T peak is only a proxy for λ peak , and the choice of the normalisation is to compare with other theoretical and observational studies. In Fig. 9 , we show how the weighted median as well as the 16-84th percentile (shaded region) peak of the IR emission varies with various galaxy properties like the stellar mass (left panel), IR luminosity (middle panel) and the specific SFR (sSFR, right panel) for z ∈ [5, 10]. In case of the variation in λ peak with galaxy stellar mass, we do not find any significant trend in the ranges we are considering. There is a general decrease (increase) in the λ peak (T peak ) value with increasing redshift. The observational data from Schreiber et al. (2018) , which shows the median data in 3.5 < z < 5 in the stellar mass range 10 11 − 10 11.5 M , is in agreement with our data at z = 5.

In a similar vein, the middle panel of Fig. 9 shows the variation of λ peak with the IR luminosity for different redshifts. We also plot is well above our values as well as the high-redshift observations. The FLARES galaxies agree well with the z = 5 fit from FIRE-2, but the slope of the relation at higher redshifts (z ≥ 7) is steeper in FLARES. There is a trend of lower (higher) values of λ peak (T peak ) with increasing IR luminosity similar to the what is seen in Ma et al. (2019) and Shen et al. (2022, with the SKIRT, postprocessed ILLUSTRIS-TNG galaxies). However, by ∼ 10 11 L we see a flattening in this relation similar to what was found in Shen et al. (2022) at z = 4, 6, with a higher normalisation than the one here. As posited there as well as in other studies (Jin et al. 2019), this trend could be due to the increasing optical depth in the most luminous galaxies, hiding the warm dust associated with star-forming compact regions, making the contribution to the dust temperature minimal. At the high IR luminosity end, there is a strong evolution towards lower (higher) λ peak (T peak ) values with redshift, showing that FLARES also prefers an evolving relation similar to results in Ma et al. (2019) ; Shen et al. (2022) .

The right panel of Fig. 9 shows the variation of λ peak with the sSFR (SFR calculated using stars born in the last 10 Myr) for different redshifts. We have also over-plotted two galaxies at z ∼ 4.5 from Burnham et al. (2021) which are in agreement with our z = 5 relation within the scatter. We see a very tight relation for the FLARES galaxies at the high sSFR (sSFR/Gyr −1 0) end. This has been observed in studies like Magnelli et al. (2014) ; Ma et al. (2019) . This strong correlation can be understood by looking at the following relation for an isothermal modified blackbody (see Hayward et al. 2011) ,

SEDs can be qualitatively described by such a form. For main sequence galaxies, L IR /M dust has been found to be proportional to the sSFR (e.g. Magdis et al. 2012; Magnelli et al. 2014; Ma et al. 2019) , implying an inverse correlation with λ peak . We show on the figure that this matches our results well at high sSFR (the T peak ∝ sSFR 1/6 dashed line). This relation can also be used to understand the increasing dust temperatures with redshift (explicit redshift evolution is shown in Fig. 10 ). Lovell et al. (2021a) has already shown that there is a systematic increase in the normalisation of the sSFR of the FLARES galaxies at constant stellar mass. This would imply that the redshift dependence of the dust temperature can be attributed to the increasing sSFR. We explore the evolution of λ peak with different stages of galaxy star-formation activity in Appendix D.

To obtain the SED dust temperature we follow Casey (2012) by parameterising our galaxy SEDs using the sum of a single modifiedblackbody and a mid-infrared powerlaw. The addition of the powerlaw to the functional form provides a better fit to the mid-infrared which is dominated by warm dust. Using this prescription, the luminosity at a rest-frame frequency ν can be written as

Here,

where h and k b are the Planck's constant and the Boltzmann constant, respectively. The free parameters in the form are N pl (normalisation factor), β (emissivity index), T SED (SED dust temperature), ν 1 (frequency where optical depth is unity, usually taken as ∼ 100 µm or 3 THz in high-redshift studies) and α (mid-IR power-law slope). We adopt the same parameterisation of ν c as given in Casey (2012) (or λ c there). We use LMFIT (Newville et al. 2014 ), a non-Linear least-squares minimization and curve-fitting package in PYTHON, to fit this parametric form to the SKIRT SEDs and obtain T SED . In the fitting, we impose the criteria that the dust temperature is higher than the CMB temperature at that redshift. It should also be noted that there is degeneracy between the dust temperature and the emissivity; a higher temperature can compensate for a lower value of the emissivity, and vice-versa. Thus when comparing to observations there is considerable maneuverability when choosing the values of T SED and β , and thus deviations or agreement with the values can also be achieved based on the ranges being probed. In our case we see that in general the SED temperature increases while the emissivity decreases with increasing redshift when they are both kept as free parameters.

The Rayleigh-Jeans (RJ) part of the SED can also be approximated by a generalised modified-blackbody function in the optically thin case (equation 2 in Casey 2012) by

where A is a normalisation constant and the rest of the terms are defined as before. In this case the free parameters are A, β (parameter search is confined to the range 1.5 − 2.5 to be consistent with high-redshift observational works referenced here) and T SED,RJ . We refer to the SED dust temperature obtained from this functional form as T SED,RJ . Similar to the fitting function in equation (5), a similar degeneracy exists here as well. This form is usually used by observational studies to derive the total infrared luminosity of galaxies (e.g. Knudsen et al. 2017; Hashimoto et al. 2019; Bouwens et al. 2020) . In many cases where there is only a single detection in the dust-continuum, β and T SED,RJ are kept constant and a fit for the normalisation is obtained. It should also be noted that the SED dust temperature obtained from this form closely matches with the galaxy peak dust temperatures, while T SED is typically higher than the peak dust temperature (see fig. 2 in Casey 2012) . From our analysis we also see that this form can lead to an overprediction of the obtained total infrared luminosity (median deviation of ∼ +17 per cent at z = 5 and lowering to 1 per cent by z = 10) compared to results using equation (5) or the true SED. This is seen when we use the full range of the dust SED. We have also tried to constrain our fits by only using wavelength ranges in the RJ tail. In this case, most fits underpredicted the total IR luminosity. Thus it is important to have some constraints at wavelengths short of the RJ tail to produce reliable SED temperature estimates that can retrieve the total IR luminosity.

In Fig. 10 we show the weighted median evolution of the different dust temperatures for the FLARES galaxies (red circles for T peak , grey squares for T SED and brown squares for T SED,RJ ). We also show the 16-84 th percentile spread of the values as well as the error on the median. Overplotted are several dust temperatures from observations of high-redshift galaxies (circle for T peak and square for T SED or T SED,RJ values, which is explicitly stated in the figure caption and the text) as well as fit functions for T peak and T SED,RJ from Liang Jin+2019 Faisst+2020 (z > 6) Hashimoto+2019 Harikane+2020 Bakx+2020 Bethermin+2020 (Stacked) Figure 10 . We show the evolution of the peak dust temperature (T peak , red circles) and the SED dust temperature (brown squares for T SED,RJ , grey diamonds for T SED ) from the FLARE simulation. The markers indicate the weighted median, the 16-84 percentile spread and the error on the median (the errorbars, negligible due to the high number counts) at z ∈ [5, 10]. Observational data from studies at high redshift (circle for T peak , diamond for T SED and square for T SED,RJ values); included are data from Strandet et al. (2016, both T peak The figure clearly indicates that the median values of the dust temperatures consistently increase towards higher redshift, as expected since higher redshift galaxies are intensely star-forming (higher sSFR), which leads to higher UV emission resulting in warmer dust. There is also a spread of ∼ 10 K in all the temperature values across the redshift range. T peak increases from a median value of ∼ 40 K at z = 5 to ∼ 70 K at z = 10, which is higher than the increase in the CMB temperature across this redshift. Using our sample of massive galaxies, we fit a linear relation to the median redshift evolution of T peak and obtain T peak /K = (40.03 ± 0.15) + (5.77 ± 0.14)(z − 5).

The relation has a higher slope compared to the one in Schreiber et al. (2018) (slope of 4.60 ± 0.35) obtained for an observational sample using stacked SEDs of main-sequence galaxies at z ≤ 4. This indicates that towards higher redshift in the EoR, the evolution of T peak is stronger for the massive galaxies in FLARES than their low redshift relation. In Fig. 10 , we also compare our values to other theoretical and observational results at similar redshifts. We compare to observa-tional values from Strandet et al. (2016); Faisst et al. (2020b) , with most values from FLARES in good agreement or otherwise within the constraints. There are a few galaxies in Strandet et al. (2016) which show slightly colder T peak values, in tension with our predictions. FLARES fails in this case to be fairly representative of such cold dusty galaxies. It can be seen that our T peak values are offset from the relation obtained from the MASSIVEFIRE simulations in Liang et al. (2019) for 2 ≤ z ≤ 6. The results are in agreement in the region for which the fit was obtained, but would be an overprediction on extrapolation. This difference could be due to the smaller sample size (29 massive galaxies in total) as well as the use of a higher dust-to-metal (DTM=0.4, increasing the optical depth) ratio in that study. Our result is also very similar to the recent values from the ILLUSTRIS-TNG (Shen et al. 2022) suite of simulations using SKIRT, with their reported median T peak being slightly higher at z = 4, 6, 8.

We can compare our T SED values to the ones in Faisst et al. (2020b) , since they use the same fitting function as in equation (5) (in their work α is fixed at 2.0, since the SED is not constrained blueward of rest-frame ∼ 110 µm). The values obtained in that study provide a very reasonable match to our constraints within the median spread, while the other observational results are lower by 10 K. The measurement from Jin et al. (2019) uses an optically thick MBB that gives dust SED temperatures that are very similar to the ones obtained from equation (5) (see Figure 2 in Casey 2012, which shows that an optically-thick MBB gives similar dust temperatures to T SED ). One of the values is in very good agreement with our measures for T SED , while the other galaxy at z ∼ 5 in the work has a very cold dust SED temperature.

For all the other measurements it would be fairer to compare the values to T SED,RJ , as equation (8) is used in many of the high-redshift studies to fit for the observed fluxes/luminosities. As such our predicted values provide a reasonable match to the observational data from Strandet et al. (2016) , Hashimoto et al. (2019) and Béthermin et al. (2020, measures the temperature using stacked galaxies with SFR ≥ 10 M /yr). There are a few exceptions in the observational data that deviate strongly from our predictions. For example, Harikane et al. (2020) fit an optically thick MBB to their galaxies at z ∼ 5, and find very cold dust SED temperatures. However, in the case of Bakx et al. (2020) the lower limit they provide is very high compared to our predictions, which could be reconciled if using a much high emmisivity index. Some of the large dispersions seen in the dust temperature measures in observations point towards either a wide range of values existing in the diverse populations of galaxies in the early Universe, or an indication of more varied dust grain properties such as their size, shape, and composition that is not captured in the models being employed to study them. The Bouwens et al. (2020) fit to the dust SED temperatures (with some of the Faisst et al. 2020b peak dust temperature values also being used), which were used in their modified blackbody fitting, does have a few values at z ≥ 5, however the majority of the constraints are from lower redshift (see their fig. 1 ). The FLARES dust SED temperatures (referring to T SED,RJ ) have slightly higher normalisation compared to their fit, but are in reasonable agreement within the 16-84th percentile spread. Similar to our results, the T SED values in Shen et al. (2022) are also consistently higher than the observationally quoted SED temperatures, with their median being slightly higher than our values.

It should also be noted that the plotted values and the spread are weighted based on which overdensity region they are from, and thus contribution from the extreme overdensities will be down-weighted. Even though they are rare, all the observational values outside the 1σ scatter are well within the maximum and minimum values of the seen dust temperatures across redshifts in FLARES, except for the lower limit measurement from Bakx et al. (2020) .

We have presented the dust SED properties of galaxies in FLARES, a suite of zoom simulations that uses the EAGLE Crain et al. 2015) physics to probe a range of overdensities in the EoR. We select massive galaxies ( 10 9 M ) in the simulation to make a comprehensive statistical study of galaxies that are accessible to current telescopes. These galaxies are post-processed with the radiative transfer code SKIRT (Camps & Baes 2015 to generate their full SEDs. The dust-to-metal ratios were derived from the fitting function from the dust model implemented in the L-Galaxies SAM (Vijayan et al. 2019 ). We do not calibrate any of the parameters in SKIRT to produce the SEDs.

Our main findings are as follows:

(i) The predicted UV LF is in agreement with available observational data. We also compare the IR LF to the observations, finding good agreement at z = 5, 6 for luminosities < 10 12 L . We underestimate the number densities of the most IR luminous galaxies at z = 5. We attribute this mainly to the lack of high star formation rates in the massive galaxies in the simulation, similar to what previous FLARES or EAGLE studies have shown (e.g. Katsianis et al. 2017; Baes et al. 2020; Lovell et al. 2021a) . We also underpredict the number of luminous rest-frame 250 µm galaxies. However, the observations Koprowski et al. (2017) ; Gruppioni et al. (2020) show discrepancies among each other. The extreme IR objects are biased towards the highest matter overdensities.

(ii) The FLARES IRX-β relation for 5 ≤ z ≤ 8 is consistent with the starburst relation (e.g. Meurer et al. 1999; Reddy et al. 2015) from local redshifts. We see a shift towards the SMC relation (Pettini et al. 1998 ) for z > 8. We see that it is predominantly the lower-mass ( 10 9 M ) galaxies that deviate towards this relation. Also galaxies with low-sSFR lie further away from these empirical relations. We see a good match with current available observations, missing a few low-IRX high-β galaxies.

(iii) The IRX shows a gradual increase with stellar mass, showing a flattening at high-stellar masses (∼ 10 10 M ). We do not see any evolution in the normalisation of the median relation with redshift.

(iv) We look at the evolution of the peak of the IR emission (λ peak ) with redshift on properties like the galaxy stellar mass, total IR luminosity and the sSFR. We see flattening of λ peak at high IR luminosity (∼ 10 11 L ). The λ peak (T peak ) -L IR relation is offset from the observed local relation (Casey et al. 2018a ) to lower (higher) values. λ peak strongly correlates with the galaxy sSFR.

(v) Luminosity-weighted dust temperatures (peak dusttemperature: T peak , SED temperature fit from mid-IR power-law+MBB: T SED and SED temperature fit from optically-thin MBB: T SED,RJ ) increase with increasing redshift. We find that, for the massive galaxies in FLARES, the evolution of T peak with redshift is stronger than the low-redshift relation obtained from observational (Schreiber et al. 2018 ) and other theoretical (Liang et al. 2019) studies.

(vi) The SED temperatures (T SED and T SED,RJ ) are mostly in agreement with the observational values. However we find a lack of extremely cold temperatures seen in some observations (Strandet et al. 2016; Jin et al. 2019 ).

Future observations from many of the planned surveys and obser-vations on ALMA, JWST, Roman Space Telescope, Euclid, as well as future IR missions sampling more of the SED will be able to put better constraints on these dust driven properties. In a future work we will explore the dust-continuum sizes of these galaxies.

Through this study as well as previous other works referenced here, there is some evidence in favour of the EAGLE physics model requiring higher star-formation rates to match some of the observations at high-redshift like the UV LF, IR LF or the sub-mm number counts. However, reconciliation of such limitations must be achieved without losing some of the remarkable successes of the model across the low-redshift Universe. The strive to succeed in this extremely non-trivial challenge has been the goal of all theoretical studies of galaxy formation and evolution. The high-redshift Universe is a regime where EAGLE as well as other periodic boxes have not been well studied due to its lack of massive galaxies. Studies with FLARES allows for a statistical exploration of this regime due its novel re-simulation strategy targeting massive overdensities. This will inevitably help to improve theoretical models of galaxy formation and evolution in terms of providing insights into the different feedback mechanisms as well as star formation recipes implemented. In this section we will test the convergence of our SEDs in relation to some of our SKIRT parameter choices. We will first explore how the choice of an increase in the photon number, increasing resolution of our radiation and dust wavelength grid, and adding Stochastic (Non Local Thermal Equilibrium, NLTE) heating to the set-up can change our results. In Fig. A1 , we plot the SED of a galaxy at z = 8 for these different configurations. Our higher photon number test runs have 10 7 (10 times higher compared to the default run) photons per radiation field wavelength grid, while for the increased wavelength grid set-up we have doubled the number of bins from our default configuration. The total IR luminosity only changes by ∼ 0.01 dex for the galaxy between these choices, while there is even smaller effect in the UV and optical part of the SED. We have checked this for a few other galaxies and find that the changes are similar. There are no noticeable dramatic changes.

We also change our dust grain distribution choice from SMC like to Milky Way (MW) like in Fig. A2 . The MW configuration has Poly-Aromatic Hydrocarbons (PAHs) included in the dust grain distribution, which will have an effect on the mid-IR range of the SED. There is a clear indication of an absorption feature near 2175 Å, due to the bump in the extinction curve for the MW dust distribution. This ultimately leads to β values that are always negative. It can be seen from the figure that there is also stronger extinction at wavelengths short of the mid-IR for the MW type, with higher emission in the mid-IR compared to SMC type. Next generation instruments that can scan the mid-IR SED at high-redshift are needed to put constraints in this regime, and to aid our understanding of emission from PAHs. The change to MW type grain has negligible effect on where the peak of the IR emission is. The total IR luminosity also sees negligible change towards higher values. However, the value of β is affected and thus can drive changes in the IRX-β plane by making the β values more negative.

In order to understand the effect AGN have on the observed galaxy SEDs, we perform a simple analysis. For this purpose, we obtain the intrinsic bolometric luminosity of the SMBH in our galaxy sample, and add it to the total IR luminosity of our galaxies. This would represent an upper limit on the total IR luminosity that the galaxy can have from AGN contribution. The SMBH bolometric luminosity is calculated using

where dM • /dt is the accretion rate and η is the efficiency, assumed to be 0.1. In Figure B1 , we plot the IR LF with the AGN bolometric luminosity added to the FLARES galaxies for z ∈ [5, 10]. For comparison we also show the IR LF (which only includes stellar reprocessed dust emission) plotted in Figure 4 We do not plot any of the observations that were shown in Figure 4 . We can see that there is a small change at the very bright end of the function. However this increase is not enough to reconcile the relation with the observational data. Thus, as explained in §3.1, the main driver of the difference can be attributed to the lack of more intense star formation activity in massive/bright galaxies in the model or the observations probing a biased region in the Universe as suggested in Zavala et al. (2021) .

To see the effect of birth cloud emission in the model, we ran SKIRT for galaxies only at z = 8 (in order to drastically reduce the computational costs, since there are only fewer galaxies at z = 8 compared to lower redshifts), by treating the young star-forming particles as regular star particles, without the added dust extinction implemented in the MAPPINGS III template due to the birth cloud. We now model both these radiation sources using BPASS, ignoring this extra extinction. We plot the UV and IR LF of the results in Figure C1 , as well as compare to our fiducial set-up. It can be seen that the UV LF is within the scatter at each bin. Similarly, in case of the IR LF at z = 8, there is only very negligible change. The small impact on both these functions is due to the low metallicity of these systems, owing to them being at extremely high-redshift. It is expected that with the increase in metallicity of systems at lower redshifts, birth cloud attenuation will have more of an influence.

To better understand the relation between λ peak and the sSFR, we separate the FLARES galaxies into 3 groups, based on an evolving piecewise fit to the stellar mass-SFR relation presented in Lovell et al. (2021a, see §3.4 in the work, equations 11 and 12). The groups have been classified based on their deviation from the piecewise fit. The 3 groups have been plotted in Figure D1 , with the median λ peak as a function of the IR luminosity. These groups are galaxies

• below 1σ (left panel in Figure D1 ) from the fit, which can be classified to include the green valley and the passive galaxies,

• within 1σ (middle panel in Figure D1 ) of the fit relation on either sides, termed the main-sequence, and

• above 1σ (right panel in Figure D1 ) of the fit, which can be termed as starbursts.

We only show in Figure D1 the relation for z ∈ [5, 9] , since the starforming sequence fit could not be constrained properly at z = 10 (see Lovell et al. 2021a) .

It can be seen from Figure D1 that the shape of the relation between λ peak and L IR is different for the 3 groups. Galaxies in the green valley/passive regime (left panel) show a consistent decrease (increase) in λ peak (T peak ) with L IR . This is mainly due to the smaller dust content within these galaxies and thus there is a direct correlation between the increase in dust-temperature and L IR . The other two groups exhibit a flat relationship with L IR , due to their high dust content and thus the hot dust being optically thick, similar to that seen in some observations at high-redshift (e.g. Cortzen et al. 2020 ). The starburst galaxies have a lower median λ peak due to their higher sS-FRs. At z = 5, towards high L IR there is also a hint of increasing . We show the variation of λ peak (corresponding T peak values is shown on the right y-axis) with the IR luminosity for z ∈ [5, 9] . The panels represent galaxies with suppresed star-formation (left), galaxies on the main-sequence (middle) and starburst galaxies (right).

λ peak values, indicating the rapid build up of dust in these extreme objects.

In this section we will compare the UV luminosity obtained from our SKIRT modelling here to the line-of-sight (LOS) dust model we implemented on the same galaxies in Vijayan et al. (2021) . In that work we assumed a dust attenuation curve and modelled the dust attenuation parameters for the old stars and the young stars to find a good match to the z = 5 UV LF and the UV-β relation, as well as observations of the [OIII] +Hβ EW relations at z = 8. This method is computationally much faster, and allows for more flexibility in the modelling, allowing you to explore changes and their effects more easily. However, it can not treat certain phenomena, such as the scat-tering of light away from the LOS or dust self-absorption, as these processes are dependent on the chosen extinction curve. We compare the UV luminosity of the galaxies selected in this work, using the two dust models in Fig. E1 , for z ∈ [5, 10]. We can see that the values obtained in this work are systematically lower than the ones from the LOS model by∼ 0.4dex. This is mainly due to the lower dust optical depth (parameterised by the κ parameters) along the LOS adopted in that study compared to this work. This by construction matches the observations that the model was calibrated for. We have already explained one of the reasons for the slightly lower number densities in §3.1 for the non-calibrated model presented here. This paper has been typeset from a T E X/L A T E X file prepared by the author.

Commercon B., Flock M., 2012, RADMC-3D: A multi-purpose radiative transfer tool

LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python

Webplotdigitizer: Version 4

Same as Fig. 3, but here we also include the UV LF obtained from the LOS dust extinction model (thicker lines) implemented in Vijayan

We thank the anonymous referee for their helpful comments. Thanks to Mark Sargent, Rebecca Bowler, Romeel Davè and Seb Oliver for helpful discussions. We also thank the ALPINE-ALMA collaboration for making their data products public. We also thank Jorge Zavala and Yueying Ni for helpful correspondences. Thanks also to Caitlin Casey for providing their updated IRX data. This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This manuscript was written during the global Covid-19 pandemic. We would like to acknowledge our privileged positions that provided us the opportunity to pursue research in these tough times, while many members of our global society cannot afford to do that. We would also like to thank everyone behind the development of various Covid-19 vaccines. APV acknowledges the support of his PhD studentship from UK STFC DISCnet. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant No. 140. CCL acknowledges support from the Royal Society under grant RGF/EA/181016.

The data associated with the figures in the paper can be found at https://flaresimulations.github.io/data.html.